libmesh中自适应求解的有限元例子

阅读libmesh官方自带的自适应的例子，主要是学习libmesh中API的使用。

阅读example/adaptivity/ex1中的例子。分析如下:

求解问题

其实不知道整个方程是如何离散化的，也就是说对应的形函数和高斯积分是如何确定的.

因为在力学上只是针对特定的PDE方程进行离散化，所以整个形函数的确定，高斯点的布置都流程都一模一样的，然后想确定PDE方程的离散化，是学习这类程序的关键，不仅仅只适用于固体的有限元程序。

代码阅读

把整个程序都会使用中文注释的方式，进行解释和说明。其实也有好多目前不太清楚的。相信了解的越多对框架的理解就越深。

// 使用自适应mesh对象进行问题求解 (求解一维问题偏微分方程问题 使用自适应单元类型)

//该问题是求解如下的一维问题：
// 问题的公式:
//  - \epsilon * u''(x) + u(x) =  1 ;  (0<= x <=1)
// 边界条件 u(0) = u(1) = 0

// 需要的头文件
// Libmesh includes
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/edge_edge3.h"
#include "libmesh/gnuplot_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/fe.h"
#include "libmesh/getpot.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/dof_map.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/error_vector.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/enum_solver_package.h"
// 自行添加的头文件 目的后面输出VTK格式的文件
#include "libmesh/vtk_io.h"

using namespace libMesh;

// 一维问题的整体单元的组装
void assemble_1D(EquationSystems & es,
                 const std::string & system_name);
// 主函数入口
int main(int argc, char ** argv)
{
  // 初始化libmesh 环境 其中包括MPI环境,petsc环境 等其他环境
  LibMeshInit init (argc, argv);
    // 这个例子是需要线性方程求解的，所以必须libMesh默认的求解包不是为空的就可以.（当在配置LIBEMSH使用了petsc时，默认使用的求解器包为petsc）
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");

  // 要求libmesh库支持自适应 (需要在配置libmesh的过程中指定，由于是自适应的程序例子 所以必须支持自适应)
#ifndef LIBMESH_ENABLE_AMR
  libmesh_example_requires(false, "--enable-amr");
#else
  
  // 其中有一个初始化 特别是并行计算时
  // 创建mesh对象 
  Mesh mesh(init.comm());

  GetPot command_line (argc, argv);// 解析命令行参数

  int n = 4;
  if (command_line.search(1, "-n"))// 查看命令行参数是否有-n 获得-n后面的值 不然默认n的值为4
    n = command_line.next(n);

  
  // 创建1d单元 在0-1范围之间创建N个单元 (单元类型为EDG3，二次型单元)
  MeshTools::Generation::build_line(mesh, n, 0., 1., EDGE3);
  // 自己添加的函数 查看mesh对象的基本属性
  // 使用一个进程运行程序 并且没有输入参数:  ./example-dbg  
  // 输入信息为:
  /*
  Mesh Information:
  elem_dimensions()={1}  # mesh对象的维度 1维
  spatial_dimension()=1  # mesh对象的空间维度 1维
  n_nodes()=9            # mesh 对象一共有 9个节点  4个单元每个单元3个节点
    n_local_nodes()=9    # mesh 对象中 当前进程 存放的节点 (由于只有一个进程)
  n_elem()=4             # mesh 对象中 总共的单元数 4 个  (参与计算一般都是使用激活的单元)
    n_local_elem()=4     # mesh 对象中 当前进程的单元数
    n_active_elem()=4    # mesh 对象中 激活的单元  (激活的单元不一定等于总共的单元)
  n_subdomains()=1       # mesh 对象中 总的子域个数 (这个干什么用 我也不知道)
  n_partitions()=1       # mesh 对象中 整个数据被切分成多少份 (分布式存放需要统计 (DistributeMesh))
  n_processors()=1       # mesh 对象中 整个程序运行的进程个数 
  n_threads()=1          # mesh 对象中 环境是否开启了多线程 线程的个数
  processor_id()=0       # mesh 对象中 输出当前信息 属于第几个进程id编号
  */
  mesh.print_info();

  //针对mesh对象定义一个方程系统 
  // 整个方程系统中，可以包括多个系统，而且这些系统可以随机被激活 ( 固体系统，温度场系统 流体场系统等)
  // 系统之间可以独立 也可以耦合
  EquationSystems equation_systems(mesh);
  
  // 在方程系统中，添加一个隐式线性求解的系统 (目的为了求解上述的方程) 并且给这个系统名称为'1D'
  LinearImplicitSystem & system = equation_systems.add_system
    <LinearImplicitSystem>("1D");

  // 上述创建了线性隐式系统 给系统添加一个新的变量 （只有有变量最后才能形成方程组）
  // ('U' 而且变量使用二阶的方程进行迭代  二阶什么之类的 目前没有搞懂)
  system.add_variable("u", SECOND);

  // 由于系统中，不同边界条件之类的情况，需要指定整体矩阵矩阵函数 同时边界条件的处理 (最终才能求解方程)
  system.attach_assemble_function(assemble_1D);

  // 定义自适应mesh对象 (参数为mesh对象)表示给该对象进行自适应
  MeshRefinement mesh_refinement(mesh);
  
  //这些参数决定了将被精炼和粗化的单元的比例。
  // 任何单元的最大误差的30%以内的任何单元都将被细化，
  // 并且任何单元上的最小误差的30%以内的任何元素都可能被粗化。 
  mesh_refinement.refine_fraction()  = 0.7;
  mesh_refinement.coarsen_fraction() = 0.3;
  // 保证所有单元都不会被细化5次
  mesh_refinement.max_h_level()      = 5;

  // 初始化方程系统 (当存在多个系统时，对里面每个系统都需要初始化)
  equation_systems.init();
  // 自己添加的函数 查看方程系统的基本信息
  // equation_systems.print_info();// 输出方程系统的基本信息

  // Refinement parameters 自适应网格参数  最大迭代步数为5次
  const unsigned int max_r_steps = 5; 

  // 开始自适应 每次根据误差决定哪些部分进行网格加密 一共求解5次
  for (unsigned int r_step=0; r_step<=max_r_steps; r_step++)
    {
      /// 求解名称为‘1D’的系统方程 中间应该会调用分配整体矩阵
      equation_systems.get_system("1D").solve();

      // 确保mesh对象在最后一次迭代的过程中，不会被细化. 因为当最后一次迭代的过程中，不需要进行网格加密了。
      if (r_step != max_r_steps)
        {
          // 误差检测向量
          ErrorVector error;// 该类包含各种统计方法的实现 实际上就是std::vector数据加上一些统计函数
          KellyErrorEstimator error_estimator;// 误差检测器

          // 这个误差向量的值 是什么值？？？
          error_estimator.estimate_error(system, error);

          // 输出误差向量的数据
          libMesh::out << "Error estimate\nl2 norm = "
                       << error.l2_norm()
                       << "\nmaximum = "
                       << error.maximum()
                       << std::endl;

          // 根据计算出来的误差向量 进行mesh对象的单元自适应
          mesh_refinement.flag_elements_by_error_fraction (error);

          // 对mesh对象进行粗化和细化
          mesh_refinement.refine_and_coarsen_elements();

          // 针对新的mesh 对象，方程系统需要重新进行初始化
          equation_systems.reinit();
        }
    }

  // 使用GnuPlot进行输出即可
  GnuPlotIO plot(mesh, "Adaptivity Example 1", GnuPlotIO::GRID_ON);

  plot.write_equation_systems("gnuplot_script", equation_systems);
  
  // 自己定义的VTK格式输出
  //VTKIO vtk(mesh);
  //vtk.set_compression(false);
  //vtk.write("out1put.pvtu");
#endif // #ifndef LIBMESH_ENABLE_AMR

  return 0;
}

// 定义一维问题的整体矩阵分配函数
void assemble_1D(EquationSystems & es,
                 const std::string & system_name)
{
  // Ignore unused parameter warnings when !LIBMESH_ENABLE_AMR.
  libmesh_ignore(es, system_name);

#ifdef LIBMESH_ENABLE_AMR

  // It is a good idea to check we are solving the correct system
  libmesh_assert_equal_to (system_name, "1D");

  // Get a reference to the mesh object
  const MeshBase & mesh = es.get_mesh();

  // The dimension we are using, i.e. dim==1 获得mesh对象的维度
  const unsigned int dim = mesh.mesh_dimension();

  // Get a reference to the system we are solving 获得方程对象的系统
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("1D");

  // Get a reference to the DofMap object for this system. The DofMap object
  // handles the index translation from node and element numbers to degree of
  // freedom numbers. DofMap's are discussed in more detail in future examples.
  // 获得这个系统中自由度管理的引用， DOfMap对象是用来处理从节点和单元编号到自由度编号的索引转换。DofMap将在以后的示例中进行更详细的讨论。
  const DofMap & dof_map = system.get_dof_map();

  // Get a constant reference to the Finite Element type for the first
  // (and only) variable in the system.
  // 获得这个系统中，第一个变量的有限元类型
  FEType fe_type = dof_map.variable_type(0);

  // Build a finite element object of the specified type. The build
  // function dynamically allocates memory so we use a std::unique_ptr in this case.
  // A std::unique_ptr is a pointer that cleans up after itself. See examples 3 and 4
  // for more details on std::unique_ptr.

  // 生成指定类型的有限元对象。build函数动态分配内存，因此在本例中使用std：：unique_ptr。
  // 一个std：：unique_ptr是一个指针，它在自身之后进行清理。有关std：：unique\u ptr的更多详细信息，请参见示例3和4。
  std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));

  // Tell the finite element object to use fifth order Gaussian quadrature
  // 告诉有限元对象使用五阶高斯求积分
  QGauss qrule(dim, FIFTH);
  fe->attach_quadrature_rule(&qrule);

  // 在这里，我们定义了一些对单元格特定数据的引用，这些数据将用于分配线性系统
  // Here we define some references to cell-specific data that will be used to
  // assemble the linear system.

  // The element Jacobian * quadrature weight at each integration point.
  //  单元雅克比*积分点上的权重(每个积分点)
  const std::vector<Real> & JxW = fe->get_JxW();

  // The element shape functions evaluated at the quadrature points.
  // 单元的形函数 在高斯点上的值
  const std::vector<std::vector<Real>> & phi = fe->get_phi();

  // The element shape function gradients evaluated at the quadrature points.
  // 单元形函数的梯度在高斯点上的值
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

  // Declare a dense matrix and dense vector to hold the element matrix
  // and right-hand-side contribution
  // 声明单元的刚度矩阵 和 右手向量
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  // This vector will hold the degree of freedom indices for the element.
  // These define where in the global system the element degrees of freedom
  // get mapped.
  // 此向量将保留单元的自由度索引。这些定义了单元自由度在全局系统中的映射位置。
  std::vector<dof_id_type> dof_indices;// 自由度到全局之间的映射关系

  // 对每个进程而言，开始对所有本地单元进行遍历为了计算单元的刚度矩阵和右手项
  // We now loop over all the active elements in the mesh in order to calculate
  // the matrix and right-hand-side contribution from each element. Use a
  // const_element_iterator to loop over the elements. We make
  // el_end const as it is used only for the stopping condition of the loop.
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // 获取当前单元的自由度索引。这些定义了这个单元在全局矩阵和右手项。
      // Get the degree of freedom indices for the current element.
      // These define where in the global matrix and right-hand-side this
      // element will contribute to.
      dof_map.dof_indices(elem, dof_indices);// 获得单元的映射关系 (单元节点第一个节点 对应 整体矩阵的位置 )
      
      // 计算当前单元特定的单元数据，这其中包括计算当前单元积分点上和形函数的值

      // Compute the element-specific data for the current element. This
      // involves computing the location of the quadrature points (q_point)
      // and the shape functions (phi, dphi) for the current element.
      fe->reinit(elem);

      // Store the number of local degrees of freedom contained in this element
      //存储此单元中包含的局部自由度的个数
      const unsigned int n_dofs =
        cast_int<unsigned int>(dof_indices.size());
      libmesh_assert_equal_to (n_dofs, phi.size());

      // We resize and zero out Ke and Fe (resize() also clears the matrix and
      // vector). In this example, all elements in the mesh are EDGE3's, so
      // Ke will always be 3x3, and Fe will always be 3x1. If the mesh contained
      // different element types, then the size of Ke and Fe would change.
      Ke.resize(n_dofs, n_dofs);
      Fe.resize(n_dofs);


      // Now loop over quadrature points to handle numerical integration
      // 对单元的高斯积分点进行循环遍历
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          // Now build the element matrix and right-hand-side using loops to
          // integrate the test functions (i) against the trial functions (j).
          for (unsigned int i=0; i != n_dofs; i++)
            {
              Fe(i) += JxW[qp]*phi[i][qp];

              for (unsigned int j=0; j != n_dofs; j++)
                {
                  Ke(i,j) += JxW[qp]*(1.e-3*dphi[i][qp]*dphi[j][qp] +
                                      phi[i][qp]*phi[j][qp]);
                }
            }
        }
      // 此时我们已经完成了单元的刚度矩阵和右手项的高斯积分点上求和。最后一步是应用边界条件，在本例中是u（0）=u（1）=0的简单Dirichlet条件。

      // At this point we have completed the matrix and RHS summation. The
      // final step is to apply boundary conditions, which in this case are
      // simple Dirichlet conditions with u(0) = u(1) = 0.

      // Define the penalty parameter used to enforce the BC's
      double penalty = 1.e10;// 大数法进行处理固定边界条件

      // Loop over the sides of this element. For a 1D element, the "sides"
      // are defined as the nodes on each edge of the element, i.e. 1D elements
      // have 2 sides.
      // 循环单元的边界，对于1维单元而言，边界定义为每条边上的节点 1D单元有两个边界
      for (auto s : elem->side_index_range())
        {
          // If this element has a nullptr neighbor, then it is on the edge of the
          // mesh and we need to enforce a boundary condition using the penalty
          // method.
          // 如果这个单元有一个nullptr邻居，那么说明它就在网格的边缘，我们需要使用惩罚方法来强制一个边界条件。
          if (elem->neighbor_ptr(s) == nullptr)
            {
              Ke(s,s) += penalty;
              Fe(s)   += 0*penalty;
            }
        }

      // This is a function call that is necessary when using adaptive
      // mesh refinement. See Adaptivity Example 2 for more details.
      // 这是使用自适应网格细化时必需的函数调用。有关详细信息，请参见自适应示例2。
      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);

      // Add Ke and Fe to the global matrix and right-hand-side.
     //将单元刚度矩阵 和 右手项 组装到整体矩阵中去
      system.matrix->add_matrix(Ke, dof_indices);
      system.rhs->add_vector(Fe, dof_indices);
    }
#endif // #ifdef LIBMESH_ENABLE_AMR
}

该例子目前还有不清楚的问题有:

1. 不知道一维整体矩阵函数分配和组装过程中，如何对该问题进行离散化的，因为这个不是传统的有限元问题，所有在离散后，不是简单的 进行计算而已。
2. 装配函数具体调用发生在哪个过程中。