引子
这个例子的亮点：介绍了最简单的Laplace方程的组装矩阵和右端项的一键生成函数。
这个例子来求解仅考虑Neumann边界条件的Laplace问题：
如果有解的话，解必须满足协调条件：
这里考虑计算域是以原点为圆心，半径为1的圆，$f=-2,g=1$是满足协调条件的。
虽然协调条件允许有解，但解仍然是不定的：在方程中仅仅解的导数固定，解可以是任意常数。所以需要施加另外一个条件来固定这个常数。可以有以下可能的方法：
- 固定离散后的某个节点值为0或其他固定值。这意味着一个额外的条件$u_h(x_0)=0$。尽管这是通常的做法，但不是一个好方法，因为我们知道Laplace方程的解是在H1空间，它不允许定义某个点的值，因为这样不是连续函数的子集。因此，尽管固定某点的值在离散时是允许的，但这样不是连续函数，结果经常是数值上的一个尖峰。
- 固定计算域上的平均值是0或其他固定值。这满足了连续性。
- 固定计算域边界上的平均值是0或者其他固定值。这也满足了连续性。
这里选择最后一种，因为还想展示另一项技术。
具体到在程序中怎样解方程，除了额外的平均值限制，其他技术都已经涉及过了，这里需要把Dirichlet边界条件删掉，至于怎样映射，在绝大多数情况下，完全可以把它当成黑盒，程序已自动处理。唯一一点就是平均值限制。幸运的是，库中有一个类知道怎样处理这种限制。如果假定边界节点沿边界均匀分布，那么平均值限制：
就可以写为：
这个求和应该遍历所有在边界上的自由度。设$i_0$是边界上指标最小的自由度，那么：
这是ConstraintMatrix类所想要的形式。注意我们之前已经多次用到过这个类，比如悬点限制：中间节点上的值等于相邻节点的平均值。通常来说，ConstraintMatrix来处理均匀限制：
其中C是矩阵，U是节点值的向量。
程序解析
1 |
之前用过的头文件。1
这个是新的。1
2
3
4
C++的标准头文件。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24namespace Step11
{
using namespace dealii;
template <int dim>
class LaplaceProblem
{
public:
LaplaceProblem (const unsigned int mapping_degree);
void run ();
private:
void setup_system ();
void assemble_and_solve ();
void solve ();
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
MappingQ<dim> mapping;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
ConstraintMatrix mean_value_constraints;
Vector<double> solution;
Vector<double> system_rhs;
TableHandler output_table;
};
这个声明跟Step5差不多，仅有少许不同。它的构造函数中接收一个映射次数的参数。1
2
3
4
5
6
7
8
9
10
11template <int dim>
LaplaceProblem<dim>::LaplaceProblem (const unsigned int mapping_degree) :
fe (1),
dof_handler (triangulation),
mapping (mapping_degree)
{
std::cout << "Using mapping with degree " << mapping_degree << ":"
<< std::endl
<< "============================"
<< std::endl;
}
构造函数，使用1阶有限单元，即fe的参数是1，映射次数是后续手动传入的。1
2
3
4
5
6template <int dim>
void LaplaceProblem<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
然后是建立系统：首先建立一个DoFHandler对象，初始化解和右端项。1
2
3
4std::vector<bool> boundary_dofs (dof_handler.n_dofs(), false);
DoFTools::extract_boundary_dofs (dof_handler,
ComponentMask(),
boundary_dofs);
然后就是建立边界上自由度上的平均值为0的限制对象。这里首先得提取出边界上的自由度。DoFTools有一个函数能返回一个布尔值的列表，其中true代表节点在边界上。第二个参数是一个蒙版，选择矢量函数的某个分量。这里处理的是一个标量有限元，所以只有一个分量。1
2
3
4
5const unsigned int first_boundary_dof
= std::distance (boundary_dofs.begin(),
std::find (boundary_dofs.begin(),
boundary_dofs.end(),
true));
然后就是具体的限制。正如引子中所说，我们限制边界上某一个节点的值等于边界上其他所有自由度的值总和。所以，先找到第一个自由度，这里是查找第一个true值，然后计算它离第一个元素的距离从而确定它的指标。1
2
3
4
5
6
7mean_value_constraints.clear ();
mean_value_constraints.add_line (first_boundary_dof);
for (unsigned int i=first_boundary_dof+1; i<dof_handler.n_dofs(); ++i)
if (boundary_dofs[i] == true)
mean_value_constraints.add_entry (first_boundary_dof,
i, -1);
mean_value_constraints.close ();
然后创建一个限制对象。首先清除所有的之前的内容，然后加上这么一行限制，在其他所有自由度的总和前加上-1的权重。1
2
3
4
5
6
7DynamicSparsityPattern dsp (dof_handler.n_dofs(),
dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dsp);
mean_value_constraints.condense (dsp);
sparsity_pattern.copy_from (dsp);
system_matrix.reinit (sparsity_pattern);
}
然后就是创建稀疏矩阵。这里也是先使用动态稀疏模式来大致估计长度，然后再创建稀疏矩阵。
接下来是组装和求解：1
2
3
4
5
6
7
8
9
10
11
12
13template <int dim>
void LaplaceProblem<dim>::assemble_and_solve ()
{
const unsigned int gauss_degree
= std::max (static_cast<unsigned int>(std::ceil(1.*(mapping.get_degree()+1)/2)),
2U);
MatrixTools::create_laplace_matrix (mapping, dof_handler,
QGauss<dim>(gauss_degree),
system_matrix);
VectorTools::create_right_hand_side (mapping, dof_handler,
QGauss<dim>(gauss_degree),
ConstantFunction<dim>(-2),
system_rhs);
首先必须组装矩阵和右端项。因为Laplace方程比较简单，库函数直接提供了工具，将单元循环什么的都集成了起来，正如上面所写。1
2
3
4
5Vector<double> tmp (system_rhs.size());
VectorTools::create_boundary_right_hand_side (mapping, dof_handler,
QGauss<dim-1>(gauss_degree),
ConstantFunction<dim>(1),
tmp);
边界力的计算也是直接调用。1
system_rhs += tmp;
然后将边界上的贡献叠加到整体右端项中。注意这个地方需要显式地叠加，因为这两个create函数都是先清除输出变量，所以不能直接叠加。1
2
3
4mean_value_constraints.condense (system_matrix);
mean_value_constraints.condense (system_rhs);
solve ();
mean_value_constraints.distribute (solution);
然后施加限制条件，再求解。1
2
3
4
5
6
7
8
9
10
11
12Vector<float> norm_per_cell (triangulation.n_active_cells());
VectorTools::integrate_difference (mapping, dof_handler,
solution,
ZeroFunction<dim>(),
norm_per_cell,
QGauss<dim>(gauss_degree+1),
VectorTools::H1_seminorm);
const double norm = norm_per_cell.l2_norm();
output_table.add_value ("cells", triangulation.n_active_cells());
output_table.add_value ("|u|_1", norm);
output_table.add_value ("error", std::fabs(norm-std::sqrt(3.14159265358/2)));
}
这里是计算了范数并输出。1
2
3
4
5
6
7
8
9
10template <int dim>
void LaplaceProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
}
然后就是求解，跟Step5相同。
run函数：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18template <int dim>
void LaplaceProblem<dim>::run ()
{
GridGenerator::hyper_ball (triangulation);
static const SphericalManifold<dim> boundary;
triangulation.set_all_manifold_ids_on_boundary(0);
triangulation.set_manifold (0, boundary);
for (unsigned int cycle=0; cycle<6; ++cycle, triangulation.refine_global(1))
{
setup_system ();
assemble_and_solve ();
};
output_table.set_precision("|u|_1", 6);
output_table.set_precision("error", 6);
output_table.write_text (std::cout);
std::cout << std::endl;
}
}
最后是main函数：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
33int main ()
{
try
{
std::cout.precision(5);
for (unsigned int mapping_degree=1; mapping_degree<=3; ++mapping_degree)
Step11::LaplaceProblem<2>(mapping_degree).run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
};
return 0;
}
计算结果
1 | Using mapping with degree 1: |
其中一个有意思的地方是：使用线性映射的误差要比使用更高阶映射的三倍要大。因此在本例中使用更高阶映射，不是因为它提高了收敛阶数，而是它直接计算出了常数。另一方面，使用三次映射并没有显著提高精度。
