dealii-step-4

最新推荐文章于 2021-10-08 14:53:36 发布

kkdut

最新推荐文章于 2021-10-08 14:53:36 发布

阅读量221

点赞数

分类专栏：有限元dealii 文章标签： c++

本文链接：https://blog.csdn.net/weixin_42312251/article/details/105890152

版权

有限元dealii 专栏收录该内容

7 篇文章 11 订阅

订阅专栏

进行任意维的编程（二维或者三维空间），即平面四边形或者六面体单元。
需要使用模板参数（template）：

template <int dim>
class Step4
{
public:
  Step4();
  void run();

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

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

  SparsityPattern      sparsity_pattern;
  SparseMatrix<double> system_matrix;

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

建立右端项和约束的默认构造函数：

template <int dim>
class RightHandSide : public Function<dim>
{
public:
  virtual double value(const Point<dim> & p,
                       const unsigned int component = 0) const override;
};

template <int dim>
class BoundaryValues : public Function<dim>
{
public:
  virtual double value(const Point<dim> & p,
                       const unsigned int component = 0) const override;
};

输入不同的维数对应的公式

在2D时 $4(x^4+y^4)$ , 在 3D时 $4(x^4+y^4+z^4)$

template <int dim>
double RightHandSide<dim>::value(const Point<dim> &p,
                                 const unsigned int /*component*/) const
{
  double return_value = 0.0;
  for (unsigned int i = 0; i < dim; ++i)
    return_value += 4.0 * std::pow(p(i), 4.0);

  return return_value;
}

边界关于维数的不同方程：在2D时 $x^2+y^2)$ , 在 3D时 $x^2+y^2+z^2)$

template <int dim>
double BoundaryValues<dim>::value(const Point<dim> &p,
                                  const unsigned int /*component*/) const
{
  return p.square();
}

接下来是使用上述函数的类模板的实现。与前面一样，我们将所有内容编写为模板，其中包含一个形式参数dim，我们假设在定义模板函数时该参数是未知的。只有在之后，编译器才会找到Step4<2>的声明(实际上是在主函数中)，然后用dim替换为2编译整个类，这个过程被称为“模板的实例化”。这样做时，它还会将RightHandSide的实例替换为RightHandSide<2>，并将类模板实例化为一个类。

实现Step4类

编译器还将在main()中发现一个声明Step4<3>。这将导致它再次返回到Step4模板，将所有出现的dim替换为3，并再次编译该类。注意，两个实例化Step4<2>和Step4<3>是完全独立的类;它们唯一的共同特征是它们都是由相同的通用模板实例化的，但是它们不能相互转换，例如，它们不共享代码(两个实例化都是完全独立编译的)。

Step4类的构造函数

它指定了所需的有限元素的多项式次数，并将DoFHandler与三角剖分关联起来：

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

创建网格

只有2维与3维相似，数据库可以进行抽象，本例中解决[−1,1]×[−1,1] 的 2D问题, 或者 [−1,1]×[−1,1]×[−1,1] 的3D问题，都可以使用GridGenerator::hyper_cube()类，

template <int dim>
void Step4<dim>::make_grid()
{
  GridGenerator::hyper_cube(triangulation, -1, 1);
  triangulation.refine_global(4);

  std::cout << "   Number of active cells: " << triangulation.n_active_cells()
            << std::endl
            << "   Total number of cells: " << triangulation.n_cells()
            << std::endl;
}

赋给矩阵自由度数

虽然形式相似，但是对于不同的维度所对应的自由度数却是相差巨大的：

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

组总纲

我们现在想要使用一个非常数的右侧函数和非零的边界值。这两个任务都是很容易实现的，只需在矩阵和右侧的汇编中添加几行新代码。我们希望有一个非常量的右端项，因此我们使用上面声明的类的一个对象来生成必要的数据。因为这个右端项只在当前函数中局部使用，所以我们在这里将它声明为一个局部变量：

 RightHandSide<dim> right_hand_side;

因为要计算非常量右端项，我们需要单元的积分点：

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

遍历所有单元并组装局部组件:

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

现在我们要把局部矩阵和右端项组装起来。这与上一个示例是完全一样的,但现在我们使用顺序循环同时处理局部矩阵和右端向量的循环，使计算快一点。组装右手端是因为其是变量：

      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

            const auto x_q = fe_values.quadrature_point(q_index);#auto由之前的定义自动分配x_q的变量类型
            cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
                            right_hand_side.value(x_q) *        // 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));

          system_rhs(local_dof_indices[i]) += cell_rhs(i);
        }

在这个例子中有非齐次的边界值，我们需要使用BoundaryValues替换函数::ZeroFunction，它描述了我们想要的边界值

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

求解器

这个步骤与维度无关，因此，跟以往的例子是一样的：

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

  std::cout << "   " << solver_control.last_step()
            << " CG iterations needed to obtain convergence." << std::endl;
}

输出结果

我们希望以VTK格式而不是gnuplot来编写输出。VTK格式是目前使用最广泛的一种格式，并得到了诸如Visit和Paraview等许多可视化程序的支持。要以这种格式写入数据，只需使用data_out.write_vtk替换data_out.write_gnuplot ：

template <int dim>
void Step4<dim>::output_results() const
{
  DataOut<dim> data_out;

  data_out.attach_dof_handler(dof_handler);
  data_out.add_data_vector(solution, "solution");

  data_out.build_patches();

  std::ofstream output(dim == 2 ? "solution-2d.vtk" : "solution-3d.vtk");
  data_out.write_vtk(output);
}

运行函数

template <int dim>
void Step4<dim>::run()
{
  std::cout << "Solving problem in " << dim << " space dimensions."
            << std::endl;
  make_grid();
  setup_system();
  assemble_system();
  solve();
  output_results();
}

主函数

增加大括号是为了在运行3d问题时，2d问题已经运行完了，可以释放了2d中的内存，从而为3d问题提供更多内存：

int main()
{
  deallog.depth_console(0);
  {
    Step4<2> laplace_problem_2d;
    laplace_problem_2d.run();
  }

  {
    Step4<3> laplace_problem_3d;
    laplace_problem_3d.run();
  }

  return 0;
}

kkdut

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
dealii-step-4

进行任意维的编程（二维或者三维空间），即平面四边形或者六面体单元。需要使用模板参数（template）：template <int dim>class Step4{public: Step4(); void run();private: void make_grid(); void setup_system(); void assemble_syste...
复制链接

扫一扫