deal.ii step-3 Lecture 10笔记

一个Laplace solver
Ω=[0,1]²
f(x)=1

类 Step3

class Step3
{
  public:
    Step3 ();
    void run ();
  private:
    void make_grid ();
    void setup_system ();
    void assemble_system ();
    void solve ();
    void output_results () const;
    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()
  : fe(1)
  , dof_handler(triangulation)//此时triangulation还未生成mesh，但不妨，它在distribute dof时才需要，现在声明是可以的
{}

void Step3::make_grid() 在之前已经介绍过

setup_system中根据sparsity_pattern矩阵创建了矩阵，并且重新初始化等式右侧解的向量，尺寸与dof总数相等
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());

void Step3::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());
}

QGuass 选择的为2d，并且因为FE_Q的fe为1，这里高斯积分的点应该是2，即每个坐标轴上两个点，共4个点

FEValues将会计算整个矩阵有等式右侧值，update_values为shape函数值，gradients of shape函数。
FEValues充当的是一个中间人的角色，

  QGauss<2> quadrature_formula(fe.degree + 1);
  FEValues<2> fe_values(fe,
                        quadrature_formula,
                        update_values | update_gradients | update_JxW_values);
  
}

定义好局部矩阵和等式右侧向量，用于后续的mapping和计算

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

边界条件，interpolate_boundary_values后面的参数为：dof_handler提取全局坐标，会被设置边界条件的边界部分，边界值函数，输出对象
一般边界部分都被指定为0

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

solve里面，1000最大的迭代次数，1e-12为tolerance。
SolverControl solver_control(1000, 1e-12);

PreconditionIdentity()表示没有Precondition

Conjugate Gradient algorithm
为SolverCG<Vector<double>> solver(solver_control)
关于扩展
using ConstantFunction<2>(1)可以表示为1的dirichlet边界值

在使用GridGenerator::hyper_cube()之后，加入下面的代码，可以使上下层边界部分指定为1

for (auto &face : triangulation.active_face_iterators())
  if (std::fabs(face->center()(1) - (-1.0)) < 1e-12 ||
      std::fabs(face->center()(1) - (1.0)) < 1e-12)
    face->set_boundary_id(1);

VectorTools::interpolate_boundary_values(dof_handler,
                                         1,
                                         ConstantFunction<2>(1.),
                                         boundary_values);