dealii-step-3

算例3介绍了泊松方程在dealii中的求解
泊松方程为边界为0，右端项非零,可以转化为AU=F：
使用类Step3开始程序的运行：

class Step3
{
public:
  Step3();
  void run();

使用一些私有对象完成对应的功能：

private:
  void make_grid();
  void setup_system();
  void assemble_system();
  void solve();
  void output_results() const;

跟step-1一样，定义网格节点变量：

  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中变量的值,默认构造函数会实现其他变量的值：

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

跟step-1一样，产生网格信息：

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

n_active_cells是指所有的网格，包括已加密的和原先的父单元。
然后将动态稀疏矩阵赋予sparsity_pattern中：

  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());

使用积分公式计算每个单元的积分：

  QGauss<2> quadrature_formula(fe.degree + 1);

需要每个单元的不同的形函数、倒数信息、雅可比行列式的权重：

  FEValues<2> fe_values(fe,
                        quadrature_formula,
                        update_values | update_gradients | 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);

使用局部编号来计算自由度，当需要把局部单元自由度耦合到全局时，需要知道该单元的全局编号，使用types::global_dof_index来定义一个临时矩阵存储该信息：

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

使用循环来遍历全部单元，因为不修改cell，可以把变量cell设为const：

  for (const auto &cell : dof_handler.active_cell_iterators())
    {

因为每个单元的雅克比的导数是不一样的，所以在遍历时，需要重置变量cell,以及全局矩阵和右端项:

      fe_values.reinit(cell);
      cell_matrix = 0;
      cell_rhs    = 0;

开始遍历所有的积分点，由于每个积分点都与其他点有关，需要两层循环，使用i和j，使用shape_grid代表积分点处的形函数导数，JxW表示权重,对右端项做同样的操作：

      for (unsigned int q_index = 0; q_index < n_q_points; ++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) +=
                (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 (unsigned int i = 0; i < dofs_per_cell; ++i)
            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);

遍历所有的单元，将单纲添加到总纲上：

      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));
#右端项
      for (unsigned int i = 0; i < dofs_per_cell; ++i)
        system_rhs(local_dof_indices[i]) += cell_rhs(i);
    }

需要把约束或者说边界条件加到方程中，否则方程有无数个解，首先获得边界的自由度和形函数，使用VectorTools::interpolate_boundary_values()插值边界函数，边界有很多种形式，不同的问题边界条件不一，需要单独添加，使用Functions::ZeroFunction（处处为0的方程）定义VectorTools::interpolate_boundary_values()，使用std::map绑定总纲和边界约束：

  std::map<types::global_dof_index, double> boundary_values;
  VectorTools::interpolate_boundary_values(dof_handler,
                                           0,
                                           Functions::ZeroFunction<2>(),
                                           boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
                                     system_matrix,
                                     solution,
                                     system_rhs);
}

求解方程，收敛条件为运行1000次或者残差范数小于：10^12，使用SolverCG模板求解，参数使用默认，CG求解器有一个前置调节器作为它的第四个参数。我们还没有准备好深入研究这个问题，所以我们使用PreconditionIdentity作为预处理：

void Step3::solve()
{
  SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> solver(solver_control);
  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}

输出结果，把数据dof_handler和solution传递给data_out，把前端数据转化为中间数据格式（后端数据），：

void Step3::output_results() const
{
  DataOut<2> data_out;
  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(solution, "solution");
  data_out.build_patches();
#写入数据
  std::ofstream output("solution.vtk");
  data_out.write_vtk(output);
}

调用类Step3中的其他函数运行：

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

主函数：

int main()
{
  deallog.depth_console(2);
  Step3 laplace_problem;
  laplace_problem.run();

  return 0;
}

dealii-step-3

推荐阅读更多精彩内容