dealii-step-5

在均匀精细网格上计算

最明显的变化是make_grid_and_dofs函数已被删除，因为创建网格现在是在运行函数中完成的，其余功能现在位于setup_system中，除此之外，一切如常。

template <int dim>
class Step5
{
public:
  Step5();
  void run();

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

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

  SparsityPattern      sparsity_pattern;
  SparseMatrix<double> system_matrix;

  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;
}

Step5类实现

template <int dim>
Step5<dim>::Step5()
  : fe(1)
  , dof_handler(triangulation)
{}

这是前面示例中的make_grid_and_dofs函数，减去网格的生成。其他一切都没有改变:

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

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

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

  system_matrix.reinit(sparsity_pattern);

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

在前面的例子中，这个函数的功能并没有太大的变化，但是仍然有一些我们将展示的优化。为此，重要的是要注意，如果使用了有效的求解器(如预处理的CG方法)，组装矩阵和右侧可能需要相当长的时间，您应该考虑在某些地方使用一两个优化。函数的第一部分与之前完全相同:

template <int dim>
void Step5<dim>::assemble_system()
{
  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);

接下来是所有单元上的典型循环，计算局部贡献，然后将它们转移到全局矩阵和向量中。与第4步相比，这部分的惟一变化是我们将使用上面定义的系数函数来计算每个求积点的系数值:

  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) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
                              1.0 *                               // f(x_q)
                              fe_values.JxW(q_index));            // dx
            }
        }


      cell->get_dof_indices(local_dof_indices);
      for (unsigned int i = 0; i < dofs_per_cell; ++i)
        {
          for (unsigned int j = 0; j < dofs_per_cell; ++j)
            system_matrix.add(local_dof_indices[i],
                              local_dof_indices[j],
                              cell_matrix(i, j));

          system_rhs(local_dof_indices[i]) += cell_rhs(i);
        }
    }

零边界值：

  VectorTools::interpolate_boundary_values(dof_handler,
                                           0,
                                           Functions::ZeroFunction<dim>(),
                                           boundary_values);
  MatrixTools::apply_boundary_values(boundary_values,
                                     system_matrix,
                                     solution,
                                     system_rhs);
}

解决方案流程与前面的示例中类似。然而，我们现在将使用一个预条件共轭梯度算法。做出这种改变并不难。事实上，我们唯一需要改变的是，我们需要一个可以作为预调节器的物体。我们将使用SSOR(对称逐次过度松弛)，其松弛因子为1.2。为此，SparseMatrix类有一个执行SSOR步骤的函数，我们需要将这个函数的地址与它所作用的矩阵(即要反转的矩阵)和松弛因子打包到一个对象中。PreconditionSSOR类为我们做了这些。PreconditionSSOR类接受一个模板参数，该参数表示它应该处理的矩阵类型。默认值是SparseMatrix，这正是我们在这里需要的，所以我们只使用默认值，不指定尖括号中的任何内容。)请注意，在目前的情况下，SSOR并没有比大多数其他预调理器表现得更好(尽管比完全没有预调理要好)。在下一个教程的结果部分，即步骤6中，将对不同的预调剂进行简要的比较。这样，函数的其余部分就很简单了:我们不再使用之前创建的PreconditionIdentity对象，而是使用之前声明的pre调节器，CG求解器将为我们完成剩下的工作:

template <int dim>
void Step5<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);

  std::cout << "   " << solver_control.last_step()
            << " CG iterations needed to obtain convergence." << std::endl;
}

step5::output_results并设置输出标记:

template <int dim>
void Step5<dim>::output_results(const unsigned int cycle) const
{
  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);
}

这个程序中倒数第二件事是run()函数的定义。与以前的程序不同，我们将在每次迭代后全局细化的网格序列上进行计算。因此，该函数由一个超过6个周期的循环组成。在每个循环中，我们首先打印循环号，然后决定如何处理网格。如果这不是第一个循环，我们只需对现有的网格进行一次全局优化。然而，在运行这些循环之前，我们必须生成一个网格:在前面的示例中，我们已经使用了GridGenerator类中的一些函数。在这里，我们希望从存储单元格的文件中读取网格，这些单元格可能来自其他人，也可能是网格生成器工具的产品。为了从文件中读取网格，我们生成了一个数据类型为GridIn的对象，并将三角关系与之关联(例如，当我们请求它读取文件时，我们告诉它填充我们的三角关系对象)。然后打开相应的文件，并使用文件中的数据初始化三角关系

template <int dim>
void Step5<dim>::run()
{
  GridIn<dim> grid_in;
  grid_in.attach_triangulation(triangulation);
  std::ifstream input_file("circle-grid.inp");

我们现在要读取文件。然而，输入文件仅用于二维三角剖分，而此函数是用于任意维度的模板。由于这只是一个演示程序，我们不会为不同的维度使用不同的输入文件，但是如果我们不是在2D中，我们会快速地杀死整个程序。当然，由于下面的主函数假设我们是在二维空间中工作，所以在这个版本的程序中，我们可以跳过这个检查，而不会有任何不良影响。事实证明，超过90%的编程错误是无效的函数参数，如无效的数组大小等，因此我们在整个过程中大量使用断言。抓住这样的错误。对于这一点，Assert宏是一个很好的选择，因为它确保作为第一个参数给出的条件是有效的，如果不是，则抛出一个异常(它的第二个参数)，该异常通常会终止程序，提供错误发生的位置和原因:

  Assert(dim == 2, ExcInternalError());

ExcInternalError是一个全局定义的异常，每当出现严重错误时，就会抛出它。通常，我们希望使用更具体的异常，特别是在这种情况下，如果dim不等于2，我们当然会尝试做一些其他的事情，例如使用库函数创建一个网格。
如果我们通过了这个断言，我们就知道了dim==2，现在我们就可以读取网格了。它是UCD(非结构化单元数据)格式的(尽管惯例是对UCD文件使用inp后缀):

  grid_in.read_ucd(input_file);

文件中的网格描述了一个圆。因此，我们必须使用一个流形对象，它告诉三角剖分在网格细化时在边界上的何处放置新点。与步骤1不同的是，由于GridIn不知道域有一个圆形边界(与GridGenerator::hyper_shell不同)，我们必须在创建三角剖分后显式地将流形附加到边界上，以便在细化网格时获得正确的结果:

  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)
    {
      std::cout << "Cycle " << cycle << ':' << std::endl;

      if (cycle != 0)
        triangulation.refine_global(1);
      std::cout << "   Number of active cells: "  //
                << triangulation.n_active_cells() //
                << std::endl                      //
                << "   Total number of cells: "   //
                << triangulation.n_cells()        //
                << std::endl;

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

现在我们有了一个网格，我们写一些输出:

int main()
{
  Step5<2> laplace_problem_2d;
  laplace_problem_2d.run();
  return 0;
}