rust的学习笔记

电气化、自动化、数字化、智能化、智慧化

0%

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

引子

在step1中创建了网格,下面就是在网格上定义自由度。此例中使用一阶线性有限元,其自由度的个数与网格的顶点数相关。后面的例子将展示更高次的单元,其上面的自由度与顶点、边、面及cell都有关。
自由度可以理解为形函数中的系数个数,因为它们是未知的,所以称之为未知量或自由度。
定义网格上的自由度很简单,因为deal.II已经内置该功能了,唯一要做的是创建有限元空间。

头文件

1
#include <deal.II/dofs/dof_handler.h>

该头文件将自由度与顶点、线、cell联系起来。

1
#include <deal.II/fe/fe_q.h>

该头文件包含双线性有限元的描述,即只在顶点上有自由度,在边上和cell内部无自由度。
1
#include <deal.II/dofs/dof_tools.h>

该头文件包含对自由度的操作工具。
1
2
3
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/dofs/dof_renumbering.h>

产生网格

这里用了step-1中的方法,只不过这里将网格triangulation作为参数返回,同时将manifold object声明为static,防止其过早销毁:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
void make_grid (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,
5 );
static const SphericalManifold<2> manifold_description(center);
triangulation.set_all_manifold_ids(0);
triangulation.set_manifold (0, manifold_description);
for (unsigned int step=0; step<3; ++step)
{
Triangulation<2>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
for (unsigned int v=0;
v < GeometryInfo<2>::vertices_per_cell;
++v)
{
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;
}
}
triangulation.execute_coarsening_and_refinement ();
}
}

创建DoFHandler

目前为止,只创建了一个网格,包含几何信息(顶点的位置)和拓扑信息(顶点怎样连成线,线连成cell,cell之间怎样连接)。为了执行数值运算,还需要一些逻辑信息,比如将自由度赋给顶点,创建矩阵和矢量,用来描述网格上的场量。
首先描述自由度是如何分布的。这里使用类模板FE_Q来创建拉格朗日单元,它的成员函数需要一个参数来描述单元的多项式次数,此处是1,表明是双线性单元,也就意味着自由度只在顶点上。如果参数是3,那么意味着是双三次单元,自由度分布为:每个顶点上一个,每条边上两个,每个cell内有四个。
示意图为:
对于Q1单元:

对于Q2单元:

对于Q3单元:

每种单元上形函数的模样可以参见这里
通过创建一个有限元对象,然后用DoFHandler分配自由度:

1
2
static const FE_Q<2> finite_element(1);
dof_handler.distribute_dofs (finite_element);

将自由度分配到每个顶点上去后,不是很容易直接可视化来看到它们,但这也不重要,因为一般情况下自由度的标号是随机的。
与网格每个顶点对应的还有形函数。注意:形函数仅在它们对应的顶点上为1,在其他顶点上则为0。那么也只相邻顶点形成的矩阵不为0,由于顶点的标号是随机的,那么总矩阵应该是稀碎的。
首先创建一个结构来存储非0元素的位置。这个类是SpasityPattern,但它有一些缺点,因为它需要事先估计每排最多有多少个,这会造成不必要的内存浪费。因此这里换用DynamicSparsityPattern这个类,传入的参数是矩阵的大小:
1
2
DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs(),
dof_handler.n_dofs());

然后根据自由度分布给出非零元素的位置:
1
DoFTools::make_sparsity_pattern (dof_handler, dynamic_sparsity_pattern);

然后将DynamicSparsityPattern的信息传回SpasityPattern:
1
2
SparsityPattern sparsity_pattern;
sparsity_pattern.copy_from (dynamic_sparsity_pattern);

然后存储到文件:
1
2
std::ofstream out ("sparsity_pattern1.svg");
sparsity_pattern.print_svg (out);

结果为:

其中每个小红方块都是矩阵中的一个非0元素。

对自由度重新编号

上面的结果可以看出,非0元素离对角线很远。对于有些算法,如不完全LU分解和Gauss-Seidel预条件子,这样的分布不好,因此需要改进。
注意:对于矩阵中非0的元素(i,j),对应的形函数i和j必须相交,而此时其所在的顶点需要相邻,因此,同一个cell内顶点的编号不能差太多才行。这可以通过一种简单的步进方法实现:首先给定一个顶点标识为0,然后对它的邻居连续标号。这里使用的是Cuthill_Mckee提出的方法:

1
2
3
4
5
6
DoFRenumbering::Cuthill_McKee (dof_handler);
DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs(),
dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dynamic_sparsity_pattern);
SparsityPattern sparsity_pattern;
sparsity_pattern.copy_from (dynamic_sparsity_pattern);

结果为: