问题描述
本案例是综合应用有限元进行问题求解的第一个案例，求解的主要方程是泊松方程（Poission's Equation）。
强形式：
按照有限元的基本理论思想，将上述方程进行有限维近似以获取有限近似解。因此在解空间内，我们获得的方程弱形式如下：
通过分部积分，我们可以得到：
所以，我们可以看到当满足的本质边界条件的同时，上述弱形式同时满足了相关的自然边界条件。
事实上，我们需要在解空间寻找满足如下关系的解：
在对真实解的近似求解需要通过采用分段线性化的手段完成，有限元中将上述用于线性化的函数称之为形函数，对于传统的伽辽金方法弱形式中的权重函数（weighting function or test function）与试函数（trial solution）采用相同的近似函数来实现。
在有限元中，近似函数实际上是一系列的多项式组合而成，形式如下：
上式中，是多项式的系数，称之为自由度，
即为形函数。
这样一来，上述的解空间即被离散，形成了如下的离散化方程：
最终得到的含有雅可比式（Jacobi）的矩阵方程组如下：
其中：
计算左手矩阵（matrix）和右手向量（vector）
对于计算域离散后得到的矩阵及向量形式如下：
对象对于单个单元而言，在单元高斯点（Gauss quadrature）上的积分离散结果：
代码解析
头文件
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
创建实体并枚举自由度
#include <deal.II/grid/grid_generator.h>
对象划分网格（在实体的基础上进行细化）
#include <deal.II/fe/fe_q.h>
有限单元类型的选择，决定了高斯点的数目,网格细化并不决定高斯点的数目，单元类型才会决定。
#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>
输出工具及一般的C++输入输出流。
using namespace dealii;
包含deal.ii 的命名空间。
程序本体
本着C++面向对象编程的基本思想和理念，最大程度的封装函数及成员。
类
class Step3
{
public:
Step3();
voidrun();
对外的公共接口包含run及构造函数。
private:
voidmake_grid();
voidsetup_system();
voidassemble_system();
voidsolve();
voidoutput_results()const;
其余函数均为私有成员。
Triangulation<2>triangulation;
FE_Q<2>fe;
DoFHandler<2>dof_handler;
SparsityPatternsparsity_pattern;
SparseMatrix<double>system_matrix;
Vector<double>solution;
Vector<double>system_rhs;
};
注意！deal.ii一个显著的特点是类模板实现，因此在使用函数及变量的过程中，有很多现成的工具只需要将模板实例化即可。因此，上述列出的类成员中，均是实例化后的对象。
函数解析
构造函数
Step3::Step3()
: fe(1)
, dof_handler(triangulation)
{}
将FE单元类型与实体关联，其他函数默认构造函数即可。
网格生成
voidStep3::make_grid()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(5);
std::cout <<"Number of active cells: "<< triangulation.n_active_cells()
<< std::endl;
}
划分网格类型为hyper_cube，细化5次，并输出总单元数目。注意triangulation.n_active_cells()指的是细化后的最终单元数目， Triangulation::n_cells()则包含了所有的单元（即划分前后的子单元及父单元数目的和）。
建立体系
voidStep3::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout <<"Number of degrees of freedom: "<< dof_handler.n_dofs()
<< std::endl;
由于在构造函数中我们已经将实体和单元类型关联，但是由于实体并没有进行网格剖分，因此，体系的自由度是未知的。
这里将已经划分完网格的实体用于自由度分配，该函数将按照网格节点枚举自由度，建立体系为进一步的求解作好准备工作。
DynamicSparsityPatterndsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
正如我们在前面的例子中看到的，我们通过首先创建一个临时结构，标记那些可能是非零的条目，然后将数据复制到 SparsityPattern 对象，此时系统矩阵可以使用该对象来设置稀疏模式。
值得注意的是， SparsityPattern 对象不保存矩阵的值，它只存储条目所在的位置。 条目本身存储在 SparseMatrix 类型的对象中，我们的变量 system_matrix 就是其中之一。
system_matrix.reinit(sparsity_pattern);
稀疏模式和矩阵之间的区别是为了允许多个矩阵使用相同的稀疏模式。 这在这里似乎无关紧要，但是当考虑矩阵可以具有的大小以及构建稀疏模式可能需要一些时间时，这在大规模问题中变得很重要。
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
此函数中要做的最后一件事是将右侧向量和解向量的大小设置为正确的值。
组装体系
组装体系是以最终形成线性方程组左右手矩阵及向量为目的的，是有限元求解的核心部分。组装矩阵和向量的一般方法是遍历所有单元格，并在每个单元格上计算该单元格对全局矩阵和右侧的正交的贡献。 现在要意识到的是，我们需要实际单元格上正交点位置处的形状函数值。 但是，有限元形状函数和正交点都仅在参考单元上定义。 因此，它们对我们几乎没有帮助，实际上我们几乎不会直接从这些对象查询有关有限元形状函数或正交点的信息。
相反，需要的是一种将这些数据从参考单元格映射到实际单元格的方法。 可以这样做的类是从 Mapping 类派生的，尽管通常不必直接处理它们：库中的许多函数可以将映射对象作为参数，但是当省略它时，它们只是求助于标准 双线性 Q1 映射。
所以我们现在拥有的是三个类的集合来处理：有限元、正交和映射对象。 这太多了，所以有一种类可以协调这三者之间的信息交换：FEValues 类。 如果给定这些对象中每三个（或两个，以及隐式线性映射）的一个实例，它将能够为您提供有关真实单元格上正交点处形状函数的值和梯度的信息。
voidStep3::assemble_system()
{
我们需要一个求积公式来评估每个单元格的积分。 让我们采用每个方向有两个正交点的高斯公式，即因为我们在 2D 中，所以总共有四个点。 这个求积公式精确地整合了最多三个度数的多项式（在 1D 中）。
QGauss<2> quadrature_formula(fe.degree + 1);
然后初始化我们上面简要讨论过的对象。需要告诉使用哪个有限元，以及正交点及其权重（由正交对象共同描述）。如前所述，我们使用隐含的 Q1 映射，而不是自己明确指定一个。最后，我们必须告诉它希望它在每个单元格上计算什么：需要正交点处的形状函数值（对于右侧 (φi,f)）、它们的梯度（对于矩阵条目（ ∇φi,∇φj))，以及正交点的权重和从参考单元到实际单元的雅可比变换的行列式。
我们实际需要什么样的信息已在FEValues 的构造函数的第三个参数给出。由于这些值必须重新计算或更新，每次我们进入一个新单元格时，所有这些标志都以前缀 update_ 开头，然后指示我们想要更新的实际内容。如果我们想要计算形状函数的值，给出的标志是 update_values；对于梯度，给出的标志是 update_gradients。雅可比行列式和正交权重总是一起使用，所以只计算乘积（雅可比乘以权重，或简称 JxW）；由于我们需要它们，我们还必须列出 update_JxW_values：
FEValues<2>fe_values(fe,
quadrature_formula,
update_values|update_gradients|update_JxW_values);
这种方法的优点是我们可以在每个单元格上指定我们实际需要的信息类型。 很容易理解，与在每个单元上计算所有内容（包括二阶导数、单元的法向量等）的方法相比，这种方法可以显着加快有限元计算，而不管是否需要它们。
该函数的第三个参数用了C语言中的位操作符，相关概念及用法可参考对应C语言内容。
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
现在，我们要逐个单元地组装全局矩阵和向量。 我们可以将结果直接写入全局矩阵，但这不是很有效，因为访问稀疏矩阵的元素很慢。 相反，我们首先计算每个单元格在具有当前单元格自由度的小矩阵中的贡献，并且只有在该单元格的计算完成后才将它们转移到全局矩阵。 我们对右手边的向量做同样的事情。 所以让我们首先分配这些对象（这些是局部对象，所有自由度都与所有其他对象耦合，我们应该使用完整矩阵对象而不是稀疏对象进行局部操作；一切都会稍后转移到全局稀疏矩阵 上）。
FullMatrix<double>cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double>cell_rhs(dofs_per_cell);
在组装每个单元格的贡献时，我们使用自由度的本地编号（即从 0 到 dofs_per_cell-1 的数字）来完成此操作。 然而，当我们将结果转移到全局矩阵中时，我们必须知道自由度的全局数。 当我们查询它们时，我们需要这些数字的临时（临时）数组。
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
从激活的网格开始，在单元内对需要的数据进行计算。
for(constauto&cell : dof_handler.active_cell_iterators())
{
我们希望计算形状函数的值和梯度，以及在正交点处参考单元格和真实单元格之间映射的雅可比矩阵的行列式。 由于所有这些值都取决于单元格的几何形状，我们必须让 FEValues 对象在每个单元格上重新计算它们：
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
对单元进行积分，通过循环所有正交点来完成，我们将通过 q_index 对其进行编号。
for(constunsignedintq_index : fe_values.quadrature_point_indices())
{
对于拉普拉斯问题，矩阵的元素计算如下所示，形函数的梯度可以通过FE_Values来查询。
for(constunsignedinti : fe_values.dof_indices())
for(constunsignedintj : fe_values.dof_indices())
cell_matrix(i, j) +=
(fe_values.shape_grad(i, q_index) *// grad phi_i(x_q)
fe_values.shape_grad(j, q_index) *// grad phi_j(x_q)
fe_values.JxW(q_index));// dx
同样的，对于残余向量做同样的计算。
for(constunsignedinti : fe_values.dof_indices())
cell_rhs(i) += (fe_values.shape_value(i, q_index) *// phi_i(x_q)
1. *// f(x_q)
fe_values.JxW(q_index));// dx
}
完成每个单元的计算贡献后，我们必须将它转移到全局矩阵和右侧。 为此，我们首先必须找出该单元格上的自由度具有哪些全局数字。 让我们简单地向单元格询问该信息：
cell->get_dof_indices(local_dof_indices);
然后再次循环遍历所有形状函数 i 和 j 并将局部元素转移到全局矩阵。 全局数字可以使用 local_dof_indices[i] 获得：
for(constunsignedinti : fe_values.dof_indices())
for(constunsignedintj : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
同样，对残余向量也需要：
for(constunsignedinti : fe_values.dof_indices())
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
对于一开始的问题，截至目前，求解的准备工作几乎全部完成。但是，作为一个唯一可解的问题，我们必须将边界条件同时加上去，这是得到问题唯一解的闭合解空间必不可少的一步。（事实上，没有 Dirichlet 边界值的拉普拉斯方程甚至不是唯一可解的，因为您可以向离散解添加任意常数）。
为此，我们首先获得边界上的自由度列表以及形状函数在那里应具有的值。 为简单起见，我们只对边界值函数进行插值，而不是将其投影到边界上。 库中有一个函数就是这样做的：VectorTools::interpolate_boundary_values()。 它的参数是（省略存在默认值且我们不关心的参数）： DoFHandler 对象，用于获取边界上自由度的全局数； 应插入边界值的边界分量； 边值函数本身； 和输出对象。
边界的组成部分含义如下：在许多情况下，可能只想在边界的一部分上强加某些边界值。 例如，在流体动力学中具有流入和流出边界，或者在体的变形计算中具有约束和自由的部分。 然后，我们需要用指标表示边界的这些不同部分，并告诉 interpolate_boundary_values 函数仅计算边界特定部分（例如约束部分或流入边界）的边界值。 默认情况下，所有边界的边界指示符为 0，除非另有说明。 如果边界的部分具有不同的边界条件，则必须使用不同的边界指示符为这些部分编号。 然后，下面的函数调用将仅确定边界指示符实际上为零指定为第二个参数的那些边界部分的边界值。
描述边界值的函数是一个 Function 类型的对象或派生类的对象。派生类之一是 Functions::ZeroFunction，它描述（不出所料）处处为零的函数。我们就地创建这样一个对象并将其传递给 VectorTools::interpolate_boundary_values() 函数。
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
求解问题
voidStep3::solve()
{
首先，我们需要有一个对象知道如何告诉 CG 算法何时停止。 这是通过使用 SolverControl 对象完成的，作为停止标准，我们说：在最多 1000 次迭代后停止（这远远超过 1089 个变量所需的数量；请参阅结果部分以了解实际使用了多少），以及 如果残差的范数低于 10−12，则停止。
SolverControl solver_control(1000, 1e-12);
然后我们需要求解器本身。 SolverCG 类的模板参数是向量的类型，保留空尖括号表示我们采用默认参数（即 Vector<double>）。 但是，我们明确提到了模板参数：
SolverCG<Vector<double>> solver(solver_control);
现在求解方程组。 CG 求解器将预处理器作为其第四个参数。
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
输出结果
典型有限元程序的最后一部分是输出结果，并可能进行一些后处理（例如计算边界处的最大应力值，或穿过流出的平均通量等）。 本文之讨论将最终的计算结果文件输出。
voidStep3::output_results() const
{
DataOut<2> data_out;
现在我们必须告诉它从哪里获取它将写入的值。 我们告诉它使用哪个 DoFHandler 对象，以及解向量（以及解变量将出现在输出文件中的名称）。 如果我们想在输出中查看多个向量（例如右侧、每个单元格的错误等），我们也会添加它们：
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution,"solution");
在 DataOut 对象知道它要处理哪些数据之后，我们必须告诉它把它们处理成后端可以处理的东西。 原因是我们将前端（知道如何处理 DoFHandler 对象和数据向量）与后端（知道许多不同的输出格式）分离，并使用中间数据格式将数据从前端传输到后端 . 数据通过以下函数转换为这种中间格式：
data_out.build_patches();
现在我们已经为实际输出准备好了一切。 只需打开一个文件并将数据写入其中，使用 VTK 格式（我们在这里使用的 DataOut 类中有许多其他函数可以将数据写入 postscript、AVS、GMV、Gnuplot 或其他一些文件格式）：
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
voidStep3::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
