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

时间：2022-09-16 06:45:30

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

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

引子

在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单元：
求解偏微分方程开源有限元软件deal.II学习--Step 2
对于Q2单元：
求解偏微分方程开源有限元软件deal.II学习--Step 2
对于Q3单元：
求解偏微分方程开源有限元软件deal.II学习--Step 2
每种单元上形函数的模样可以参见这里。
通过创建一个有限元对象，然后用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
 
SparsityPattern sparsity_pattern;
sparsity_pattern.copy_from (dynamic_sparsity_pattern);

然后存储到文件：

1
2
 
std::ofstream out ("sparsity_pattern1.svg");
sparsity_pattern.print_svg (out);

结果为：
求解偏微分方程开源有限元软件deal.II学习--Step 2
其中每个小红方块都是矩阵中的一个非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);

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

 #deal.II
标签：

相关文章