完整版见：qixinbo.info.
deal.II的程序结构

deal.II采用面向对象编程，其中包含了很多的Modules，各自实现不同的功能，并有机地结合起来。如上图所示。具体为：

1. Triangulation
   Triangulations是单元及其更低维度的边界的集合。triangulation存储了网格的几何和拓扑性质：单元怎样接触，它们的顶点在哪里。triangulation不知道将要在它上面使用的有限元的任何信息，甚至不知道它自己的单元的形状，它只知道在二维情形下有4条线段和4个顶点，三维下有6个面、12条线段和8个顶点。不过其他所有信息都定义在映射类mapping中，由该类将参考单元的坐标映射到真实单元的坐标上，通常采用线性映射。
   当需要访问triangulation的性质和数据时，通过使用iterators迭代器对单元进行循环。

2. Finite Element
   Finite element类用来描述定义在参考单元上的有限元空间的性质，比如在单元的顶点、线段或内部各有多少自由度，此外还给出节点上的形函数的值及其梯度。

3. Quadrature
   跟Finite element相同，Quadrature也定义在单元上，用来描述参考单元上积分点的位置及其权重。

4. DoFHandler
   DoFHandler对象是triangulations和finite elements的汇合点：finite element类描述了triangulation单元的点、线或内部需要多少自由度，而DoFHandler分配了这种空间，从而使得点、线或内部都有正确的数目，同时也给这些自由度统一编号。
   也可以这样理解：triangulation和finite element描述了有限元空间的具体性质(有限元空间就是用来得到离散解的空间)，DoFHandler则枚举了该空间的基本框架，使得我们可以用一系列有序的系数Uj来表示离散解uh(x)=∑jUjϕi(x)。
   正如triangulation对象，DoFHandler也可以通过iterators迭代器对单元进行循环，从而得到具体的信息，比如单元的几何和拓扑信息(这些也可以由triangulation迭代获得，其实两个类是派生关系)，以及当前单元上的自由度数目和数值。需要注意的是，跟triangulation一样，DoFHandler也不知道从参考单元到它上面的真实单元的映射，也不知道对应于它所管理的自由度的形函数。

5. Mapping
   Mapping类就是建立从参考单元到triangulation的单元的映射，包括形函数、积分点和积分权重，同时提供了这种派生的梯度和Jacobian行列式。

6. FEValues
   这一步就是真正地取出某个单元，计算它上面在积分点上的形函数值及其梯度。注意，有限元空间不是连续的，积分都是在特定积分点上计算。

7. Linear Systems
   一旦知道了怎样使用FEValues在单个单元上计算形函数值及其梯度，同时知道怎样使用DoFHandler获得自由度的全局标识，那么就可以使用矩阵类和向量类来存储和管理矩阵和向量的元素，从而形成线性系统。

8. Linear Solvers
   构建好线性系统后，就可以使用求解器来求解该系统。

9. Output
   求解完毕后，就可以输出结果到可视化软件中进行后处理。

源码解析头文件

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

这两个头文件分别处理triangulation和自由度。

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

该文件用于生成网格。

1 2 3

#include<deal.II/grid/tria_accessor.h> #include<deal.II/grid/tria_iterator.h> #include<deal.II/dofs/dof_accessor.h>

这三个文件用来对单元进行循环并获得单元上的信息。

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

该文件包含拉格朗日插值的描述。

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

#include<deal.II/dofs/dof_tools.h> #include<deal.II/fe/fe_values.h> #include<deal.II/base/quadrature_lib.h> #include<deal.II/base/function.h> #include<deal.II/numerics/vector_tools.h> #include<deal.II/numerics/matrix_tools.h> #include<deal.II/lac/vector.h> #include<deal.II/lac/full_matrix.h> #include<deal.II/lac/sparse_matrix.h> #include<deal.II/lac/dynamic_sparsity_pattern.h> #include<deal.II/lac/solver_cg.h> #include<deal.II/lac/precondition.h> #include<deal.II/numerics/data_out.h> #include<fstream> #include<iostream>

以上头文件还需仔细理解。

step3类

跟之前两个例子不同，这次将信息都封装在了一个类里面。

1 2 3 4 5

class Step3 { public:        Step3 (); voidrun();

这是类的public部分，包含一个构造函数和一个run函数，用来说明执行顺序。

1 2 3 4 5 6

private: voidmake_grid(); voidsetup_system(); voidassemble_system(); voidsolve(); voidoutput_results()const;

这是类的private部分的成员函数，函数名说明了其要实现的功能。

1 2 3 4 5 6 7 8

Triangulation<2> triangulation; FE_Q<2> fe; DoFHandler<2> dof_handler; SparsityPattern sparsity_pattern; SparseMatrix<double> system_matrix; Vector<double> solution; Vector<double> system_rhs; };

还有成员变量，用来存储各种信息。
以下是step3类的各个成员函数的详解。

1 2 3 4 5

Step3::Step3 ()    :    fe (1),    dof_handler (triangulation) {}

这是step3类的构造函数，它没有执行具体操作，只是调用了成员初始器对fe和dof_handler进行了初始化。fe是一个finite element对象，它接收1，1是多项式的次数，表明使用的是双线性插值的形函数。
之前的triangulation也被传递给了dof_handler。注意此时triangulation还没有具体建立网格，但dof_handler不介意，它只有在具体分配自由度时才关心网格。
下一步必须做的就是对计算域进行剖分，然后对每个顶点分配自由度。

1 2 3 4 5 6 7 8

void Step3::make_grid () {    GridGenerator::hyper_cube (triangulation, -1, 1);    triangulation.refine_global (5); std::cout << "Number of active cells: "        << triangulation.n_active_cells()        << std::endl; }

这个函数做的是第一步，即对计算域进行剖分，建立网格。计算域是[−1,1]x[−1,1]。
因为初次建立时只有一个单元，所以细化5次，形成1024个单元，这里用n_active_cells()验证一下个数。注意这里用的不是n_cells()函数，因为其还包涵父单元的概念。
下一步就是分配自由度：

1 2 3 4 5 6

void Step3::setup_system () {    dof_handler.distribute_dofs (fe); std::cout << "Number of degrees of freedom: "        << dof_handler.n_dofs()        << std::endl;

这里使用distribute_dofs()函数，接收的是fe，因为fe是线性插值，所以自由度是每个顶点上有一个。用n_dofs()输出一下，显示是1089。这是因为我们有32x32个网格，那么对应是33x33个节点。
然后创建稀疏矩阵

1 2 3 4

DynamicSparsityPattern dsp(dof_handler.n_dofs()); DoFTools::make_sparsity_pattern (dof_handler, dsp); sparsity_pattern.copy_from(dsp); system_matrix.reinit (sparsity_pattern);

注意，SparsityPattern和SparseMatrix不同，前者只存储元素的位置，后者则存储具体的数值。
然后建立右端项向量和解向量：

1 2 3

solution.reinit (dof_handler.n_dofs()); system_rhs.reinit (dof_handler.n_dofs()); }

下一步就是计算线性系统中的矩阵中的元素以及右端项的元素，这是每个有限元程序的核心部分！
组装矩阵和向量通常的方法是对所有单元进行循环，然后在每个单元上进行积分运算，得到该单元对整体的贡献。需要注意的是此时需要知道在每个真实单元上形函数在积分点位置上的值，但是！！形函数和积分点都是仅仅定义在参考单元上的，因此这些东西基本没用。那么关键问题来了，就是怎样将数据从参考单元上映射到真实单元上。这个活是由Mapping类来完成的，更加智能的是通常不需要人为指定它怎么做，它自动按标准双线性映射来操作。
现在我们需要处理三类东西：有限元finite element、积分quadrature和映射mapping。这些概念太多了，这里有一个类FEValues来将三者有机地整合起来。给它传进去这三个东西，它就能告诉你在真实单元上的形函数的值和梯度。
那么现在就开干吧：

1 2 3

void Step3::assemble_system () {    QGauss<2> quadrature_formula(2);

先确定在单元上的一套积分规则，这里使用的是Gauss数值积分方法，每个方向上选两个积分点。这套积分规则满足现在的问题。然后就可以生成FEValues的对象了：

1 2	FEValues<2> fe_values (fe, quadrature_formula,        update_values | update_gradients | update_JxW_values);

第一个参数告诉它参考单元是谁，第二个参数告诉它积分点及其权重(实际是一个Qudrature对象)，还有默认使用了双线性映射。最后告诉它需要在每个单元上计算什么，包括在积分点上的形函数值(为了计算右端项ϕi(xKq)f(xKq))、在积分点上的形函数梯度(为了计算矩阵元素∇ϕi(xKq)⋅∇ϕj(xKq))、Jacobian行列式。注意：这里都是积分点上的值。因为需要对单元进行循环，所有这些值需要更新，所以加上前缀update。这样显著地跟程序说明需要计算什么，就能加速运算，因为有的软件所有东西不管有用没用都一块计算，比如二阶导数、法向量等。注意最后用了按位或这一运算符，这在c语言中很常见，这里的意思就是我要计算谁“和”谁。

1 2	constunsignedint dofs_per_cell = fe.dofs_per_cell; constunsignedint n_q_points = quadrature_formula.size();

然后定义两个“快捷名”来称呼常用的两个变量：每个单元上的自由度个数和积分点个数。
现在终于开始一个单元一个单元地组装整体刚度矩阵和向量了。一种方法是直接将结果写入整体矩阵中，但这样对于大型矩阵的运算通常是很慢的。所以，这里是先在当前单元上计算该单元的单元刚度矩阵，然后将它的贡献叠加到整体上。计算右端项向量时也是这样的。
首先先创建以上两种单元矩阵和单元向量：

1 2	FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell); Vector<double> cell_rhs (dofs_per_cell);

当计算每个单元的贡献时，是对该单元上的自由度的局部标识进行循环，即从0到dofs_per_cell-1。
然而，当把结果传递给整体时，需要知道这些自由度的全局标识，这时需要一个临时的量来存储：

1	std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

来，现在真的要对单元进行循环：

1 2 3 4 5

DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),     endc = dof_handler.end(); for (; cell!=endc; ++cell) {

现在，我们整个人站在单元上，我们想要知道在积分点上的形函数的值及其梯度以及参考单元和真实单元之间变换的Jacobian行列式。因为这些值与每个单元的几何信息有关，所以必须在每个单元上都需要让FEValues重算一下这些东西：

1	fe_values.reinit (cell);

初始化单元的贡献为0，以便后面的赋值：

1 2	cell_matrix = 0; cell_rhs = 0;

然后开始积分，对积分点进行循环：

1 2	for (unsignedint q_index=0; q_index<n_q_points; ++q_index) {

首先是对矩阵的组装：

1 2 3 4 5

for (unsignedint i=0; i<dofs_per_cell; ++i) for (unsignedint j=0; j<dofs_per_cell; ++j)        cell_matrix(i,j) += (fe_values.shape_grad (i, q_index) *                fe_values.shape_grad (j, q_index) *                fe_values.JxW (q_index));

对于拉普拉斯方程，每个单元上的矩阵是形函数i和j的梯度的积分，程序中用quadrature代替integral，所以变成两个梯度的乘积乘以Jacobian行列式(这是为了从参考单元到真实单元)和权重。两个梯度是矢量，它们的乘积是一个点乘，是一个标量。
然后对右端项向量的组装：

1 2 3 4 5

for (unsignedint i=0; i<dofs_per_cell; ++i)    cell_rhs(i) += (fe_values.shape_value (i, q_index) * 1 *            fe_values.JxW (q_index)); }

这里的积分是形函数的值乘以f乘以JxW。
现在有了单元的矩阵，下一步就是把它组装到整体上。我们先得问问这个单元，它的某个局部标号的自由度的全局标识是多少：

1	cell->get_dof_indices (local_dof_indices);

后面就可以通过local_dof_indices[i]来获取全局标识。
然后再作循环叠加：

1 2 3 4 5 6 7 8

for (unsignedint i=0; i<dofs_per_cell; ++i) for (unsignedint j=0; j<dofs_per_cell; ++j) system_matrix.add (local_dof_indices[i],        local_dof_indices[j],        cell_matrix(i,j)); for (unsignedint i=0; i<dofs_per_cell; ++i) system_rhs(local_dof_indices[i]) += cell_rhs(i); }

这样，整个线性系统基本全做完了。But！还有一个重要的东西漏了：边界条件。事实上，如果该拉普拉斯方程不加上一个Dirichlet边界条件，它的解就不是惟一的，因为你可以在解上加一个任意的量。显然得解决这个问题：

1 2 3 4 5

std::map<types::global_dof_index,double> boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0,        ZeroFunction<2>(),        boundary_values);

这里使用的interpolate_boundary_values函数就是将函数值插入到特定位置的边界上，它需要的参数有：DoFHandler对象来获得边界上的自由度的全局标识;边界的一部分;边界上的值本身;输出对象。
”边界的一部分“意思是有时你只想在边界的一部分上赋予边界值，比如流体力学的入流和出流边界、变形体的固定部分等。这时可以对边界进行一部分一部分地标号，不同的标号施加不同的条件。默认条件下所有的边界标号都是0，这里也使用的是0，表明是对整个边界作用。描述边界值的函数这里用的是ZeroFunction，返回的是0，刚好适用现在的零边界条件。最后的输出对象boundary_values是一个列表，里面是成对的边界自由度全局标识及其对应的边界值。
知道上述信息后，再实际施加边界条件：

1 2 3 4 5 6 7 8 9 10

MatrixTools::apply_boundary_values (boundary_values,        system_matrix,        solution,        system_rhs); } 完全建立好线性系统了，终于该求解了。这个问题的变量个数是1089，实际是很小的量，通常的问题一般是10万个变量这个量级。传统的求解线性代数的算法，如Gauss消去或LU分解，对于大型系统不适用，这里用的是共轭梯度算法，即CG算法。 ```cpp void Step3::solve () { SolverControl solver_control(1000, 1e-12);

首先告诉CG算法何时停止计算，这里是创建一个SolverControl对象来控制：最多迭代1000步，精度为1e-12。通常这个精度值是迭代停止的判据。

1	SolverCG<> solver (solver_control);

然后创建一个CG的求解器。

1 2 3

solver.solve (system_matrix, solution, system_rhs,        PreconditionIdentity()); }

上面就是开始求解了。第四个参数是一个预条件子。求解完毕后，solution中就存储了解向量的离散值。
然后就是输出结果了：

1 2 3

void Step3::output_results () const {    DataOut<2> data_out;

首先让它知道从哪里提取数据：

1 2	data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution");

知道了目标数据后，离后面的可输出数据还隔了一个“中间数据”，因为deal.II是前后端分离的，这个中间数据像是个缓冲层，得到它的语句为：

1	data_out.build_patches ();

然后再输出成最终的可被可视化软件读取的数据：

1 2 3

std::ofstream output("solution.gpl"); data_out.write_gnuplot (output); }

上面一系列的成员函数通过run函数来按次序执行：

1 2 3 4 5 6 7 8

void Step3::run () {    make_grid ();    setup_system ();    assemble_system ();    solve ();    output_results (); }

该程序的main函数为：

1 2 3 4 5 6 7

intmain() {    deallog.depth_console (2);    Step3 laplace_problem;    laplace_problem.run (); return0; }

第一个语句是打印log信息到屏幕上，后面就是创建step3对象，然后执行。

计算结果

输出到屏幕上的信息有：

1 2 3 4

Number of active cells: 1024 Number of degrees of freedom: 1089 DEAL:cg::Starting value 0.121094 DEAL:cg::Convergence step 48 value 5.33692e-13

即，说明了单元个数1024、自由度个数1089，CG算法的起始残差是0.12，经过47步迭代后满足精度要求。
图形结果为：

转载本文请联系原作者获取授权，同时请注明本文来自亓欣波科学网博客。
链接地址：https://blog.sciencenet.cn/blog-441611-998557.html

上一篇：利用Mathematica进行有限元编程(一)：杆单元分析
下一篇：deal.II链接PETSc过程记录