dealii-step-6

使用网格自适应技术

步骤：
1.求解当前网格上的PDE;
2.使用一些指示误差的标准来估计每个单元的误差;
3.把那些误差大的细胞标记为“精化”，把那些误差特别小的细胞标记为“粗化”，其余的就不用管了;
4.细化和粗化标记的细胞，以获得一个新的网格;
5.在新网格上重复上述步骤，直到总体误差足够小。
自适应网格中四边形或者六面体网格会出现悬挂点的情况，；使用网格自适应主要是因为有些局部使用较稀疏的网格就能达到所需要的精度，没有必要增加网格数量，增加计算量。

dealii处理悬挂点

步骤：
1.创建一个AffineConstraints对象，该对象将存储有限元空间中的所有约束。在有悬挂点的情况下，这些约束可以解决方案空间的连续。
2.使用函数DoFTools::make_hanging_node_constraints()来确保有限元空间中的元素具有连续性。
3.使用AffineConstraints::distribute_local_to_global()将矩阵和右端项复制到全局对象中，这个函数的作用是确保悬挂点上的自由度是被约束的。通过将它们的行和列设为0并在对角线上添加一些东西以确保矩阵保持可逆，它们实际上从线性系统中被消除了。这个过程得到的矩阵对于我们解的拉普拉斯方程仍然是对称的正定的，所以我们可以继续使用共轭梯度法。
4.然后像往常一样求解线性系统，但是在这一步的最后，您需要确保悬挂点上的“自由度”得到正确的(受约束的)值，以便您随后以其他方式可视化或求值的解决方案实际上是连续的。这是通过调用AffineConstraints::distribute()来完成的。

如何得到局部细化的网格

    利用误差估计器对拉普拉斯算子数值解的能量误差进行估计，得到局部精细网格，实现它的类称为KellyErrorEstimator。该类计算一个向量时，其条目的数量与活动单元格的数量相同，并且每个条目都包含对该单元格上的错误的估计。这个估计值然后被用来细化网格的单元格:那些有较大误差的单元格将被标记为细化，那些有特别小的估计值的单元格将被标记为粗化。我们不必手动完成这些工作，一旦获得了误差估计向量，函数GridRefinement将为我们完成所有这些工作。

边界条件

AffineConstraints类也可以处理这样的约束，这使得我们可以使用与悬挂点约束相同对象来处理Dirichlet边界条件。这样，我们就不需要在组装之后再定义边界条件。我们所需要做的就是调用VectorTools:: interpolate_boundary_values()的变体，它在AffineConstraints对象中返回它的信息。

程序讨论

使用FE_Q定义双二次元素

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

定义悬挂点的连续性，保证全局的连续性，需要定义悬挂点的约束：

#include<deal.II/lac/affine_constraints.h>

为了在本地细化网格，我们需要来自库的一个函数，该函数根据我们计算的错误指示器来决定标记哪些单元格进行细化或粗化:

#include<deal.II/grid/grid_refinement.h>

我们需要一种简单的方法来实际计算基于误差估计的细化指标。一般来说，自适应能力是针对特定问题的，但以下文件中的错误指示符通常会生成非常适合于大量问题的网格:

#include<deal.II/numerics/error_estimator.h>

Step6类模板

template <int dim>
class Step6
{
public:
  Step6();

  void run();

private:
  void setup_system();
  void assemble_system();
  void solve();
  void refine_grid();
  void output_results(const unsigned int cycle) const;

  Triangulation<dim> triangulation;

  FE_Q<dim>       fe;
  DoFHandler<dim> dof_handler;

我们需要一个包含约束列表的对象来保存挂起节点和边界条件:

AffineConstraints<double> constraints;
SparseMatrix<double> system_matrix;
SparsityPattern    sparsity_pattern;
Vector<double> solution;
Vector<double> system_rhs;
};

系数变量

template <int dim>
double coefficient(const Point<dim> &p)
{
  if (p.square() < 0.5 * 0.5)
    return 20;
  else
    return 1;
}

Step6类实现

该类的构造函数与以前的构造函数基本相同，但这一次我们要使用二次元。为了做到这一点，我们只需要用所需的多项式次数(这里是2)替换构造函数参数(在前面的所有例子中都是1):

template <int dim>
Step6<dim>::Step6()
         :fe(2)
        ,dof_handler(triangulation)
{}

Step6::setup_system

首先是要处理悬挂点问题，还有就是二次单元具有9个自由度：

template <int dim>
void Step6<dim>::setup_system()
{
  dof_handler.distribute_dofs(fe);

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

现在，我们可以使用悬挂点约束填充affineconstraints对象。因为我们将在一个循环中调用这个函数，我们首先清除当前系统的约束，然后计算新的约束:

  constraints.clear();
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);

我们准备使用指示器0(整个边界)插入边界值，并将结果约束存储在我们的constraints对象中。注意，我们不像在前面的步骤中那样，在组装之后应用边界条件:而是在AffineConstraints对象中对函数空间施加所有约束。我们可以按任意顺序向AffineConstraints对象添加约束:如果两个约束冲突，则约束矩阵可以通过Assert宏中止或抛出异常:

  VectorTools::interpolate_boundary_values(dof_handler,
                                           0,
                                           Functions::ZeroFunction<dim>(),
                                           constraints);

在添加了所有约束之后，需要对它们进行排序和重新安排，以便更有效地执行某些操作。这个后处理是使用close()函数完成的，在此之后，就不能再添加任何约束了:

constraints.close();

现在我们首先构建我们的压缩稀疏模式，就像我们在前面的例子中所做的那样。然而，我们不会立即将其复制到最终的稀疏模式。注意，我们调用make_sparsity_pattern的一个变体，它将AffineConstraints对象作为第三个参数。我们通过将参数keep_constrained_dofs设为false(换句话说，我们将永远不会写入对应于受约束自由度的矩阵项)来让例程知道我们将永远不会写入由约束条件给出的位置。如果我们要在装配后压缩约束条件，我们将不得不传递true，因为我们将首先写入这些位置，然后在压缩期间再次将它们设置为零。

DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler,
                                  dsp,
                                  constraints,
                                  /*keep_constrained_dofs = */ false);

现在矩阵的所有非零元素都是已知的(即那些通过定期组装矩阵和那些通过消除约束引入的元素)。我们可以将中间对象复制到稀疏模式:

  sparsity_pattern.copy_from(dsp);
  system_matrix.reinit(sparsity_pattern);
}

Step6::assemble_system

我们要重新组装矩阵。与第5步相比，有两个代码变化:首先，我们必须使用高阶求积公式来解释有限元形状函数的高多项式次数。这很容易改变:QGauss类的构造函数接受每个空间方向上的正交点的数量。以前，对于双线性元素我们有两个点。现在我们应该用三个点来表示双二次元。其次，为了将每个单元上的局部矩阵和向量复制到全局系统中，我们不再使用手工编写的循环。相反，我们使用AffineConstraints::distribute_local_to_global()，它在内部执行这个循环，同时对与自由度上的约束度对应的行和列执行高斯消除。构成本地贡献的其余代码保持不变。

template <int dim>
void Step6<dim>::assemble_system()
{
  const QGauss<dim> quadrature_formula(fe.degree + 1);

  FEValues<dim> fe_values(fe,
                          quadrature_formula,
                          update_values | update_gradients |
                            update_quadrature_points | update_JxW_values);

  const unsigned int dofs_per_cell = fe.dofs_per_cell;
  const unsigned int n_q_points    = quadrature_formula.size();

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

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

  for (const auto &cell : dof_handler.active_cell_iterators())
    {
      cell_matrix = 0;
      cell_rhs    = 0;

      fe_values.reinit(cell);

      for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
        {
          const double current_coefficient =
            coefficient<dim>(fe_values.quadrature_point(q_index));
          for (unsigned int i = 0; i < dofs_per_cell; ++i)
            {
              for (unsigned int j = 0; j < dofs_per_cell; ++j)
                cell_matrix(i, j) +=
                  (current_coefficient *              // a(x_q)
                   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

              cell_rhs(i) += (1.0 *                               // f(x)
                              fe_values.shape_value(i, q_index) * // phi_i(x_q)
                              fe_values.JxW(q_index));            // dx
            }
        }

      cell->get_dof_indices(local_dof_indices);
      constraints.distribute_local_to_global(
        cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
    }

Step6::solve

求解线性系统的函数再次使用了SSOR预调节器，并且再次保持不变，只是必须包含悬挂点约束。如前所述，通过对矩阵的行和列进行特殊处理，可以从线性系统中去除与悬挂点约束和边界值对应的AffineConstraints对象的自由度。这样，这些自由度的值就错了，没有约束悬挂点。然后我们要做的是使用约束来为它们分配它们应该有的值。这个过程称为分布约束，它从无约束节点的值计算受约束节点的值，并且只需要在这个函数末尾找到一个附加的函数调用：

template <int dim>
void Step6<dim>::solve()
{
  SolverControl            solver_control(1000, 1e-12);
  SolverCG<Vector<double>> solver(solver_control);

  PreconditionSSOR<SparseMatrix<double>> preconditioner;
  preconditioner.initialize(system_matrix, 1.2);

  solver.solve(system_matrix, solution, system_rhs, preconditioner);

  constraints.distribute(solution);#增加的一项

Step6::refine_grid

我们使用一个复杂的误差估计方案来细化网格，而不是全局细化。我们将使用KellyErrorEstimator类，它实现了拉普拉斯方程的误差估计;估计器的工作方式是将描述自由度的DoFHandler对象和每个自由度的值向量作为输入，并为三角化的每个活动单元计算单个指示值(即每个活动单元一个值)。

template <int dim>
void Step6<dim>::refine_grid()
{
  Vector<float> estimated_error_per_cell(triangulation.n_active_cells());

  KellyErrorEstimator<dim>::estimate(dof_handler,
                                     QGauss<dim - 1>(fe.degree + 1),
                                     {},
                                     solution,
                                     estimated_error_per_cell);

上面的函数为estimated_error_per_cell数组中的每个单元格返回一个错误指示符值。精化现在按照以下步骤进行:精化具有最高错误值的30%的单元格，粗化具有最低错误值的3%的单元。人们很容易验证,如果第二个数为零,这将大约导致加倍的细胞每一步在两维空间,因为每个单元的30%,四个新也会被取代,而剩下的70%的单元保持不变。在实践中，通常会产生更多的单元，因为不允许一个细胞单元两次，而相邻的单元不精炼;在这种情况下，相邻的单元也会被细化。在许多应用程序中，要粗化的单元格数量将设置为比3%更大的值。非零值非常有用，特别是在初始(粗化)网格已经相当精细的情况下。在这种情况下，可能有必要在某些区域对其进行细化，而在其他区域进行粗化是有用的。在我们的例子中，初始网格非常粗糙，因此粗化只在少数可能发生过精化的区域才有必要。因此，一个小的、非零的值在这里是合适的。下面的函数现在获取这些细化指示器，并使用上面描述的方法标记三角化的一些单元格，以便细化或粗化。它来自实现了几种不同算法的类，这些算法基于单元格方向的错误指示器来细化三角剖分。

  GridRefinement::refine_and_coarsen_fixed_number(triangulation,
                                                  estimated_error_per_cell,
                                                  0.3,
                                                  0.03);

在前一个函数退出后，一些单元格被标记为细化，另一些标记为粗化。然而，细化或粗化本身现在还没有执行，因为在某些情况下，对这些标志的进一步修改是有用的。在这里，我们不想做任何这样的事情，因此我们可以告诉三角定位来执行标记单元格的操作:

 triangulation.execute_coarsening_and_refinement();

Step6::output_results

在每个网格的计算结束时，在我们继续下一个网格细化循环之前，我们希望输出这个循环的结果。我们已经在步骤1中看到了如何实现网格本身。我们惟一需要更改的是文件名的生成，因为它应该包含作为参数提供给这个函数的当前细化周期的数量。为此，我们只需将细化周期的编号作为字符串追加到文件名:

template <int dim>
void Step6<dim>::output_results(const unsigned int cycle) const
{
  {
    GridOut               grid_out;
    std::ofstream         output("grid-" + std::to_string(cycle) + ".gnuplot");
    GridOutFlags::Gnuplot gnuplot_flags(false, 5);
    grid_out.set_flags(gnuplot_flags);
    MappingQGeneric<dim> mapping(3);
    grid_out.write_gnuplot(triangulation, output, &mapping);
  }

  {
    DataOut<dim> data_out;
    data_out.attach_dof_handler(dof_handler);
    data_out.add_data_vector(solution, "solution");
    data_out.build_patches();

    std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
    data_out.write_vtu(output);
  }
}

Step6::run

函数主循环中的第一个块处理网格生成。如果这是程序的第一个周期，而不是像前面的示例那样从磁盘上的文件读取网格，那么我们现在再次使用库函数创建它。域也是一个圆，这就是为什么我们还要提供一个合适的边界对象。我们将圆心放在原点，半径为1(这是函数的两个隐藏参数，具有默认值)。通过查看粗网格，您将注意到它的质量比我们在前一个示例中从文件中读取的网格要差:单元格的形状不太均匀。但是，使用库函数这个程序可以在任何空间维度下工作，这是以前没有的情况。如果我们发现这不是第一个周期，我们想要细化网格。与上一个示例程序中使用的全局优化不同，我们现在使用上面描述的自适应过程:

template <int dim>
void Step6<dim>::run()
{
  for (unsigned int cycle = 0; cycle < 8; ++cycle)
    {
      std::cout << "Cycle " << cycle << ':' << std::endl;

      if (cycle == 0)
        {
          GridGenerator::hyper_ball(triangulation);
          triangulation.refine_global(1);
        }
      else
        refine_grid();


      std::cout << "   Number of active cells:       "
                << triangulation.n_active_cells() << std::endl;

      setup_system();

      std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
                << std::endl;

      assemble_system();
      solve();
      output_results(cycle);
    }
}

main()函数

与前面的示例相比，主函数的功能没有改变，但是我们采取了额外的谨慎步骤。有时，会出现一些问题(例如写入输出文件时磁盘空间不足，试图分配向量或矩阵时内存不足，或者由于某种原因无法读写文件)，在这种情况下，库会抛出异常。由于这些是运行时问题，而不是可以一次性修复的编程错误，所以这种异常不会在优化模式下关闭，这与我们用来针对编程错误进行测试的Assert宏形成了对比。如果未捕获，这些异常将把调用树传播到主函数，如果也未捕获，则程序将中止。在许多情况下，比如没有足够的内存或磁盘空间，我们不能做任何事情，但至少可以打印一些文本来解释程序失败的原因。:

int main()
{
  try
    {
      Step6<2> laplace_problem_2d;
      laplace_problem_2d.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;
}