dealii-step-8

求解弹性方程-ElasticProblem类

该类与step-6大体一致，唯一的变化是对fe变量使用了一个不同的类:我们现在使用的是一个更通用的类FESystem，而不是像FE_Q这样的具体的有限元类。事实上，FESystem本身并不是一个真正的有限元，因为它没有实现自己的形状函数。相反，它是一个可用于将多个其他元素堆叠在一起，形成一个向量值的有限元类。在我们的例子中，我们将组合FE_Q(1)对象的向量值元素，如下面的这个类的构造函数所示

  template <int dim>
  class ElasticProblem
  {
  public:
    ElasticProblem();
    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;
    DoFHandler<dim>    dof_handler;

    FESystem<dim> fe;

    AffineConstraints<double> constraints;

    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;

    Vector<double> solution;
    Vector<double> system_rhs;
  };

右端项

在转到主类的实现之前，我们声明并定义了描述右侧的函数。为了防止返回向量之前没有被设置为正确大小的情况，我们测试了这种情况，否则在函数开始时抛出一个异常：

  template <int dim>
void right_hand_side(const std::vector<Point<dim>> &points,
                     std::vector<Tensor<1, dim>> &  values)
{
  Assert(values.size() == points.size(),
         ExcDimensionMismatch(values.size(), points.size()));
  Assert(dim >= 2, ExcNotImplemented());

函数的其余部分实现计算载荷。我们将x方向上和y方向上使用一个单位力在，受力位于圆形点(0.5,0)和(-0.5,0);在3d中，这些中心的z分量也为零。
因此，我们首先定义两个表示这些区域中心的对象。注意，在构造Point对象时，所有值都设置为0，如果点在受力点0.2内，该点方向向量值变为1：

    Point<dim> point_1, point_2;
   point_1(0) = 0.5;
   point_2(0) = -0.5;

   for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
     {
       if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) ||
           ((points[point_n] - point_2).norm_square() < 0.2 * 0.2))
         values[point_n][0] = 1.0;
       else
         values[point_n][0] = 0.0;
         
       if (points[point_n].norm_square() < 0.2 * 0.2)
         values[point_n][1] = 1.0;
       else
         values[point_n][1] = 0.0;
     }
 }

ElasticProblem 类——ElasticProblem::ElasticProblem

下面是主类的构造函数。如前所述，我们希望构造一个向量值的有限元，它由若干个标量有限元(即，我们希望构建向量值元素，使其每个向量分量都包含一个标量元素的形状函数)。当然，我们想要叠加在一起的标量有限元元素的数量等于解函数的分量的数量，这是模糊的，因为我们考虑了在每个空间方向上的位移。FESystem类可以处理这个问题:我们将希望组成系统的有限元素传递给它，以及它应该重复多少次：

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

实际上，FESystem类有更多的构造函数，它们可以执行更复杂的操作，而不仅仅是将相同类型的几个标量有限元叠加成一个。

ElasticProblem::setup_system

建立方程组与步骤6的例子中使用的函数相同。DoFHandler类和这里使用的所有其他类都充分意识到，我们想要使用的有限元是向量值的，并且要注意有限元本身的向量值。

  template <int dim>
 void ElasticProblem<dim>::setup_system()
 {
   dof_handler.distribute_dofs(fe);
   solution.reinit(dof_handler.n_dofs());
   system_rhs.reinit(dof_handler.n_dofs());

   constraints.clear();
   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
   VectorTools::interpolate_boundary_values(dof_handler,
                                            0,
                                            Functions::ZeroFunction<dim>(dim),
                                            constraints);
   constraints.close();

   DynamicSparsityPattern dsp(dof_handler.n_dofs(), 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);
 }

ElasticProblem::assemble_system

这个程序最大的变化是矩阵和右手边的创建，因为它们是依赖于问题的。我们将逐步完成这个过程，因为它比前面的示例更复杂一些。但是，这个函数的第一部分与前面相同:设置一个合适的求积公式，初始化我们使用的(向量值)有限元和求积对象的FEValues对象，并声明一些辅助数组。此外，我们声明两个相同的缩写:n_q_points和dofs_per_cell。我们现在要求的是每个单元的自由度，而不是底层的标量Q1单元的自由度。

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

我们需要一个地方来存储单元格上所有求积点上的系数值。在目前的情况下，我们有两个系数lambda，mu

    std::vector<double> lambda_values(n_q_points);
    std::vector<double> mu_values(n_q_points);

我们也可以忽略上面的两个数组，因为我们将对和使用常系数，它可以这样声明。它们都表示总是返回常数值1.0的函数。

Functions::ConstantFunction<dim> lambda(1.), mu(1.);

像上面的两个常量函数一样，我们将对每个单元格只调用一次right_hand_side函数，以使事情变得更简单。

  std::vector<Tensor<1, dim>> rhs_values(n_q_points);

现在我们可以开始对所有单元循环,求出求积分点处的系数值:

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

        fe_values.reinit(cell);
        lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
        mu.value_list(fe_values.get_quadrature_points(), mu_values);
        right_hand_side(fe_values.get_quadrature_points(), rhs_values);

然后将局部刚度矩阵和右侧向量的项集合起来。这几乎是一对一的，与本例介绍中描述的模式相同。我们可以使用fe.system_to_component_index(i)计算数字comp(i)，即形状函数i的惟一非零向量组件的索引，组装局部矩阵。

        for (unsigned int i = 0; i < dofs_per_cell; ++i)
          {
            const unsigned int component_i =
              fe.system_to_component_index(i).first;

            for (unsigned int j = 0; j < dofs_per_cell; ++j)
              {
                const unsigned int component_j =
                  fe.system_to_component_index(j).first;

                for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
                  {
                    cell_matrix(i, j) +=
                                          (                                                  //
                        (fe_values.shape_grad(i, q_point)[component_i] * //
                         fe_values.shape_grad(j, q_point)[component_j] * //
                         lambda_values[q_point])                         //
                        +                                                //
                        (fe_values.shape_grad(i, q_point)[component_j] * //
                         fe_values.shape_grad(j, q_point)[component_i] * //
                         mu_values[q_point])                             //
                        +   
                            ((component_i == component_j) ?        
                           (fe_values.shape_grad(i, q_point) * //
                            fe_values.shape_grad(j, q_point) * //
                            mu_values[q_point]) :              //
                           0)                                  //
                        ) *                                    //
                      fe_values.JxW(q_point);                  //
                  }
              }
          }

组装右端项

        for (unsigned int i = 0; i < dofs_per_cell; ++i)
       {
         const unsigned int component_i =
           fe.system_to_component_index(i).first;

         for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
           cell_rhs(i) += fe_values.shape_value(i, q_point) *
                          rhs_values[q_point][component_i] *
                          fe_values.JxW(q_point);
       }

从局部自由度到全局矩阵和右手边向量的转换不依赖于所考虑的方程，因此与前面的所有例子都是一样的。一旦我们完成了整个线性系统的装配，同样的方法也适用于消除悬挂节点:

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

边界值的插值需要一个小的修改:因为解函数是向量值的，所以需要是边界值。函数Functions::ZeroFunction构造函数接受一个参数，该参数告诉它应该表示一个有这么多分量的向量值的常零函数。默认情况下，这个参数等于1，在这种情况下，function::ZeroFunction对象表示一个标量函数。由于解向量有dim分量，我们需要将dim作为分量的数量传递给零函数。

ElasticProblem::solve

解算器并不关心方程组从何而来，只要它保持正定和对称(这是CG解算器的要求)，方程组确实是这样的。因此，我们不需要改变任何事情

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

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

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

   constraints.distribute(solution);
 }

ElasticProblem::refine_grid

对网格进行细化的函数与step6中的示例相同。求积公式再次适用于线性元素。注意，在默认情况下，误差估计器将从有限元解的所有分量中得到的估计值相加，即，它使用在所有方向上的位移具有相同的权重。如果我们想让网格只适应x位移，我们可以给函数传递一个额外的参数，告诉它这样做，而不考虑误差指示器的所有其他方向的位移。然而，对于目前的问题，似乎应该考虑所有的位移分量的权重相等。

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

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

   triangulation.execute_coarsening_and_refinement();
 }

ElasticProblem::output_results

唯一的区别是解函数是向量值。DataOut类自动处理这个问题，但是我们必须为解决方案向量的每个组件提供一个不同的名称。为此，DataOut::add_vector()函数需要一个字符串向量。由于组件的数量与我们正在处理的维度的数量相同，所以我们使用下面的switch语句。

  template <int dim>
void ElasticProblem<dim>::output_results(const unsigned int cycle) const
{
  DataOut<dim> data_out;
  data_out.attach_dof_handler(dof_handler);

  std::vector<std::string> solution_names;
  switch (dim)
    {
      case 1:
        solution_names.emplace_back("displacement");
        break;
      case 2:
        solution_names.emplace_back("x_displacement");
        solution_names.emplace_back("y_displacement");
        break;
      case 3:
        solution_names.emplace_back("x_displacement");
        solution_names.emplace_back("y_displacement");
        solution_names.emplace_back("z_displacement");
        break;
      default:
        Assert(false, ExcNotImplemented());
    }

在为解决方案向量的不同组件设置名称之后，我们可以将解决方案向量添加到计划输出的数据向量列表中。请注意，下面的函数将字符串向量作为第二个参数，而我们在前面所有示例中使用的那个函数在这里接受一个字符串。(事实上，我们之前使用的函数会将单个字符串转换成只有一个元素的向量，并将其转发给另一个函数)

 data_out.add_data_vector(solution, solution_names);
 data_out.build_patches();

 std::ofstream output("solution-" + std::to_string(cycle) + ".vtk");
 data_out.write_vtk(output);

ElasticProblem::run

我们使用square [-1,1]^d作为域，并在开始第一次迭代之前全局加密两次,加密两次的原因:我们使用了QGauss求积公式，在每个方向上都有两个点，用于右手边的积分;这意味着在每个单元上有四个求积点(在2D中)。如果我们只对初始网格进行一次全局优化，那么在域上每个方向上只有四个求积点。然而,右边函数被选中,左边单元为约束，在这种情况下,用求积法计算的右手边的向量将只包含0(即使它当然是非0的，如果我们精确地用积分计算右手边的向量)，方程组的解是0向量，即一个处处为零的有限元函数。

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

      if (cycle == 0)
        {
          GridGenerator::hyper_cube(triangulation, -1, 1);
          triangulation.refine_global(4);
        }
      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

int main()
{
  try
    {
      Step8::ElasticProblem<2> elastic_problem_2d;
      elastic_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;
}