算例3介紹了泊松方程在dealii中的求解
泊松方程為邊界為0立磁，右端項(xiàng)非零,可以轉(zhuǎn)化為AU=F：
使用類Step3開始程序的運(yùn)行：

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

使用一些私有對象完成對應(yīng)的功能：

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

跟step-1一樣，定義網(wǎng)格節(jié)點(diǎn)變量：

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

由拉普拉斯方程離散化得到的系統(tǒng)矩陣的稀疏模式和值的變量:

  SparsityPattern      sparsity_pattern;
  SparseMatrix<double> system_matrix;

方程右邊的解和變量：

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

傳入類Step3中變量的值,默認(rèn)構(gòu)造函數(shù)會(huì)實(shí)現(xiàn)其他變量的值：

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

跟step-1一樣，產(chǎn)生網(wǎng)格信息：

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是指所有的網(wǎng)格沽瞭，包括已加密的和原先的父單元。
然后將動(dòng)態(tài)稀疏矩陣賦予sparsity_pattern中：

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

為了區(qū)別對待稀疏矩陣與矩陣，轉(zhuǎn)換稀疏矩陣為矩陣,這對于大規(guī)模計(jì)算非常重要：

  system_matrix.reinit(sparsity_pattern);

最后是設(shè)置方程右手邊的向量和解向量的大惺└摇：

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

使用積分公式計(jì)算每個(gè)單元的積分：

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

需要每個(gè)單元的不同的形函數(shù)、倒數(shù)信息、雅可比行列式的權(quán)重：

  FEValues<2> fe_values(fe,
                        quadrature_formula,
                        update_values | update_gradients | update_JxW_values);

賦值一些常用的數(shù)值為變量僵娃，方便以后的重用以及更改,每個(gè)單元自由度個(gè)數(shù)以及求積分的點(diǎn)數(shù)目：

  const unsigned int dofs_per_cell = fe.dofs_per_cell;
  const unsigned int n_q_points    = quadrature_formula.size();

使用矩陣計(jì)算每個(gè)單元的單綱概作，最后再放到總綱的稀疏矩陣中,以及右邊的向量：

  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
  Vector<double>     cell_rhs(dofs_per_cell);

使用局部編號(hào)來計(jì)算自由度，當(dāng)需要把局部單元自由度耦合到全局時(shí)默怨，需要知道該單元的全局編號(hào)讯榕，使用types::global_dof_index來定義一個(gè)臨時(shí)矩陣存儲(chǔ)該信息：

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

使用循環(huán)來遍歷全部單元，因?yàn)椴恍薷腸ell匙睹，可以把變量cell設(shè)為const：

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

因?yàn)槊總€(gè)單元的雅克比的導(dǎo)數(shù)是不一樣的愚屁，所以在遍歷時(shí)，需要重置變量cell,以及全局矩陣和右端項(xiàng):

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

開始遍歷所有的積分點(diǎn)痕檬，由于每個(gè)積分點(diǎn)都與其他點(diǎn)有關(guān)霎槐，需要兩層循環(huán)，使用i和j谆棺，使用shape_grid代表積分點(diǎn)處的形函數(shù)導(dǎo)數(shù)栽燕，JxW表示權(quán)重,對右端項(xiàng)做同樣的操作：

      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
        }

組總綱，需要單元自由度的總數(shù)：

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

需要把約束或者說邊界條件加到方程中碍岔，否則方程有無數(shù)個(gè)解，首先獲得邊界的自由度和形函數(shù)朵夏，使用VectorTools::interpolate_boundary_values()插值邊界函數(shù)蔼啦，邊界有很多種形式，不同的問題邊界條件不一仰猖，需要單獨(dú)添加捏肢，使用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);
}

求解方程饥侵，收斂條件為運(yùn)行1000次或者殘差范數(shù)小于：10^12鸵赫，使用SolverCG模板求解，參數(shù)使用默認(rèn)躏升，CG求解器有一個(gè)前置調(diào)節(jié)器作為它的第四個(gè)參數(shù)辩棒。我們還沒有準(zhǔn)備好深入研究這個(gè)問題，所以我們使用PreconditionIdentity作為預(yù)處理：

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

輸出結(jié)果膨疏，把數(shù)據(jù)dof_handler和solution傳遞給data_out一睁，把前端數(shù)據(jù)轉(zhuǎn)化為中間數(shù)據(jù)格式（后端數(shù)據(jù)），：

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();
#寫入數(shù)據(jù)
  std::ofstream output("solution.vtk");
  data_out.write_vtk(output);
}

調(diào)用類Step3中的其他函數(shù)運(yùn)行：

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

主函數(shù)：

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

  return 0;
}