求解偏微分方程开源有限元软件deal.II学习--Step 1

求解偏微分方程开源有限元软件deal.II学习--Step 1

Posted on 2016-08-02   |   Incomputational material science   |   暂无评论

简介

deal.II是一款开源的求解偏微分方程的有限元软件，它有如下几个特点：

• 使用C++编写
• 有多种单元类型
• 可以大规模并行
• 可以自适应网格
• 文档和范例齐全
• 与其他库有良好的接口

安装

deal.II最新版本为8.4.1,可从官网上下载源码，解压后进入源文件目录安装：

1 2 3 4 5

mkdir build cd build cmake -DCMAKE_INSTALL_PREFIX=/path/to/install/dir ../deal.II make install make test

如果期间需要其他依赖如cmake、doxygen等，自行安装好即可。

编译运行

deal.II的文档和example特别完备，刚开始接触可以它多达54个例子进行。
具体入口在这里。
编译运行命令为：

1 2 3

cmake . make make run

第一条命令用来创建Makefile文件，指明程序所依赖的文件、怎样编译和运行。此命令应该能找到之前安装deal.II后的库文件，如果不能找到，需要人为指定路径：

1	cmake -DDEAL_II_DIR=/path/to/installed/deal.II .

第二条命令将源文件编译成可执行文件，第三条是运行该可执行文件。其实可以省略第二条命令。后面的example遵循同样的命令。

第一个教学实例——step-1

1	#include <deal.II/grid/tria.h>

此头文件声明了Triangulation类模板，其用途是生成各种单元。如Triangulation<1,1>表示一维线段，Triangulation<1,2>或<2,3>表示二维中的曲线和三维中的面，通常用于边界单元中。

1 2	#include <deal.II/grid/tria_accessor.h> #include <deal.II/grid/tria_iterator.h>

1	#include <deal.II/grid/grid_generator.h>

该头文件声明了GridGenerator命名空间，用于生成标准网格

1	#include <deal.II/grid/manifold_lib.h>

1	#include <deal.II/grid/grid_out.h>

该头文件声明了GridOut类，用于生成多种格式的数据，如dx、gnuplot、msh、eps、svg、vtk等。

1	using namespace dealii;

为了防止命名冲突，该包的函数和类都包含在dealii这个命名空间中。

1 2	GridGenerator::hyper_cube (triangulation); triangulation.refine_global (4);

用超立方体来初始化给定的Triangulation，并将网格加密4次。
第一个函数first_grid生成一个由一个单元网格加密四次后的16×16=256个正方形单元的网格。

1	GridGenerator::hyper_shell (triangulation,center,inner_radius,outer_radius,10);

生成一个三维的壳体或二维的环。注意：这里的边界是弯曲的，而默认情形下边界都是直线。这时需要指定“manifold indicator”来指明边界，告诉在哪个地方细化。
如果不指明manifold indicator，那么就只会得到在周向为10个cell，然后加密3次：

1 2 3 4 5 6 7 8

Triangulation<2> triangulation; const Point<2> center (1,0); const double inner_radius = 0.5, outer_radius = 1.0; GridGenerator::hyper_shell (triangulation, center, inner_radius, outer_radius, 10); triangulation.refine_global (3);

结果如图：

对边界进行特定标识后，效果会好一些：

1 2 3 4 5 6 7 8 9 10

Triangulation<2> triangulation; const Point<2> center (1,0); const double inner_radius = 0.5, outer_radius = 1.0; GridGenerator::hyper_shell (triangulation, center, inner_radius, outer_radius, 10); const HyperShellBoundary<2> boundary_description(center); triangulation.set_boundary (0, boundary_description); triangulation.refine_global (3);

结果如图：

可以看出此时能较好地表现出内外边界的圆形特征，但仍然能从切线上的弯折分辨出最初的10个cell。这可以通过不光对边界上的线，而是对所有的线都指定indicator来优化：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Triangulation<2> triangulation; const Point<2> center (1,0); const double inner_radius = 0.5, outer_radius = 1.0; GridGenerator::hyper_shell (triangulation, center, inner_radius, outer_radius, 10); const SphericalManifold<2> boundary_description(center); triangulation.set_manifold (0, boundary_description); Triangulation<2>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); for (; cell!=endc; ++cell) cell->set_all_manifold_ids (0); triangulation.refine_global (3);

效果如图：

之前设定的hyper_shell初始周向网格为10个cell，如果设置为3,且只对边界上的cell设定indicator：

1 2 3 4 5 6 7 8 9 10

Triangulation<2> triangulation; const Point<2> center (1,0); const double inner_radius = 0.5, outer_radius = 1.0; GridGenerator::hyper_shell (triangulation, center, inner_radius, outer_radius, 3); // four circumferential cells const HyperShellBoundary<2> boundary_description(center); triangulation.set_boundary (0, boundary_description); triangulation.refine_global (3);

结果为：

可以看出细化效果很差，但即使初始为3,也可以通过对全部的cell设定indicator来优化：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Triangulation<2> triangulation; const Point<2> center (1,0); const double inner_radius = 0.5, outer_radius = 1.0; GridGenerator::hyper_shell (triangulation, center, inner_radius, outer_radius, 3); // three circumferential cells const SphericalManifold<2> boundary_description(center); triangulation.set_manifold (0, boundary_description); Triangulation<2>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); for (; cell!=endc; ++cell) cell->set_all_manifold_ids (0); triangulation.refine_global (3);

结果为：

以上分析过程见这里。

在step-1的second_grid函数中，将上述细化更进一步，不再是全局细化，而是局部细化，因此首先是要得到能指向每个cell的指针，可以想象一个triangulation是所有cell的集合，而cell在其中并不是一个序列，因此这里不能直接用指针，而是用迭代器iterator，从第一个cell开始，遍历所有cell。不过这里没有使用遍历所有cell的迭代器，而是使用只遍历active cell的迭代器active_cell_iterator，active cell是没有children的cell，其后来将被细化：

1 2 3 4

Triangulation<2>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); for (; cell!=endc; ++cell)

然后在这个for循环中再遍历每个cell的所有顶点：

1 2 3

for (unsigned int v=0; v < GeometryInfo<2>::vertices_per_cell; ++v)

然后从这些顶点中通过判断该顶点与圆心的距离找到属于内边界的顶点，从而标记该顶点所在的cell，用于后续的细化：

1 2 3 4 5 6 7

const double distance_from_center = center.distance (cell->vertex(v)); if (std::fabs(distance_from_center - inner_radius) < 1e-10) { cell->set_refine_flag (); break; }

然后开始执行细化操作：

1	triangulation.execute_coarsening_and_refinement ();

结果为：

注意：当函数结束，开始销毁所创建的对象时，按相反顺序进行，因为定义的manifold object是在triangulation之后，所以manifold object先销毁，此时将会报错，因此必须先释放它：

1	triangulation.set_manifold (0);

另一种简单的方法是在triangulation之前定义manifold object。

second_grid给了我们提示，可以设定不同的细化条件，比如设定当cell的中心的y坐标大于0时才加密：

1 2 3

for (; cell!=endc; ++cell) if (cell->center()[1] > 0) cell->set_refine_flag ();

结果为：

最后

deal.II的作者建议尽早学会使用一个debugger！

http://blog.sciencenet.cn/home.php?mod=space&uid=441611

http://qixinbo.info/tags/deal-II/

向原作者

亓欣波

致敬！！！！