ellipticequation.cpp

基于前面上传的变系数的椭圆型方程的例子

#include "EllipticEquation.h"

/**
 * 计算二次项系数的函数的表达式。
 */
double _a_(const double * p)

{
  if (p[0]*p[0] > p[1]*p[1])
    return 4;
  else
    return 1;
}

/**
 * 计算矩阵 Matrix 的单元刚度矩阵，未做改变。
 */
void
Matrix::getElementMatrix(const Element<double,DIM>& ele0,
                         const Element<double,DIM>& ele1,
                         const ActiveElementPairIterator<DIM>::State state)
{
  double volume = ele0.templateElement().volume();
  const QuadratureInfo<DIM>& quad_info =
    ele0.findQuadratureInfo(algebricAccuracy());
  std::vector<double> jacobian =
    ele0.local_to_global_jacobian(quad_info.quadraturePoint());
  int n_quadrature_point = quad_info.n_quadraturePoint();
  std::vector<Point<DIM> > q_point =
    ele0.local_to_global(quad_info.quadraturePoint());
  std::vector<std::vector<std::vector<double> > > basis_grad =
    ele0.basis_function_gradient(q_point);
  int n_ele_dof = eleDof0().size();
  for (int l = 0;l < n_quadrature_point;l ++) {
    double Jxw = quad_info.weight(l)*jacobian[l]*volume;
    double a_val = _a_(q_point[l]);
    for (int j = 0;j < n_ele_dof;j ++) {
      for (int k = 0;k < n_ele_dof;k ++) {
        elementMatrix(j,k) += Jxw*a_val*innerProduct(basis_grad[j][l],
                                                     basis_grad[k][l]);
      }
    }
  }
}

/**
 * 初始化函数也基本没有改变，但是在读入网格数据文件的时候，我们现在需
 * 要使用 MovingMesh 类提供的 readDomain 函数。这个函数会不但读入网格
 * 剖分文件，并且会对区域描述文件 D.d 和 D.log 进行分析。注意：D.d 中
 * 的节点附近的剖分尺度在这个部分是不会用到的。
 */
void EllipticEquation::initialize()
{
  triangle_template_geometry.readData("triangle.tmp_geo");
  triangle_coord_transform.readData("triangle.crd_trs");
  triangle_template_dof.reinit(triangle_template_geometry);
  triangle_template_dof.readData("triangle.1.tmp_dof");
  triangle_basis_function.reinit(triangle_template_dof);
  triangle_basis_function.readData("triangle.1.bas_fun");

  template_element.resize(1);
  template_element[0].reinit(triangle_template_geometry,
                             triangle_coord_transform,
                             triangle_template_dof,
                             triangle_basis_function);

  /// 读入网格数据文件和区域描述文件
  this->readDomain("D");
}

/**
 * 构造有限元空间。我们只需要做一次这个操作，即使是网格移动了以后，有
 * 限元空间也不必更新。事实上，网格移动了以后，有限元空间中每个自由度
 * 的插值点会不对。但是由于我们下面的计算并不使用到插值操作，因此就不
 * 必更新自由度的插值点了。
 */
void EllipticEquation::buildSpace()
{
  fem_space,reinit(*this, template_element);
        
  int n_element = mesh.n_geometry(DIM);
  fem_space.element().resize(n_element);
  for (int i = 0;i < n_element;i ++) {
    fem_space.element(i).reinit(fem_space,i,0);
  }

  fem_space.buildElement();
  fem_space.buildDof();
  fem_space.buildDofBoundaryMark();
}

/**
 * 构造线性系统，并求解，得到需要的有限元函数。在这个函数中，我们仅仅
 * 对边界条件的部分进行了更新，因为我们为了能够应用移动网格模块，更新
 * 了区域描述文件中的边界材料标识。
 */
void EllipticEquation::solve()
{
  Matrix stiff_matrix(fem_space);
  stiff_matrix.algebricAccuracy() = 3;
  stiff_matrix.build();

  u_h.reinit(fem_space);

  Vector<double> f_h;
  Operator::L2Discretize(&_f_, fem_space, f_h, 3);

  /// 使用狄氏边值条件。现在我们需要对四个边分别添加边值条件。
  BoundaryConditionAdmin<double,DIM> boundary_admin(fem_space);
  BoundaryFunction<double,DIM> boundary1(BoundaryConditionInfo::DIRICHLET,
                                         1, &_u_);
  BoundaryFunction<double,DIM> boundary2(BoundaryConditionInfo::DIRICHLET,
                                         2, &_u_);
  BoundaryFunction<double,DIM> boundary3(BoundaryConditionInfo::DIRICHLET,
                                         3, &_u_);
  BoundaryFunction<double,DIM> boundary4(BoundaryConditionInfo::DIRICHLET,
                                         4, &_u_);
  boundary_admin.add(boundary1);
  boundary_admin.add(boundary2);
  boundary_admin.add(boundary3);
  boundary_admin.add(boundary4);
  boundary_admin.apply(stiff_matrix, u_h, f_h);

  AMGSolver solver(stiff_matrix);
  solver.solve(u_h, f_h);      
}

/**
 * run 是整个程序运行的主要驱动引擎。这里是最外层的循环。
 */
void EllipticEquation::run()
{
  initialize(); /// 完成准备工作
  buildSpace(); /// 建立有限元空间

  /**
   * 调用 MovingMesh 的 moveMesh 函数进行网格移动，主循环事实上隐藏在
   * 这个函数中。
   */
  moveMesh();
}

/**
 * 这个函数计算控制网格移动的控制函数。我们使用解的梯度的阶跃来作为控
 * 制量。由于我们现在使用最简单的线性逼近的方式，我们可以手工在计算这
 * 个量，避免要进行面单元上的数值积分带来的麻烦。请仔细理解下面这段代
 * 码。
 */
void EllipticEquation::getMonitor()
{
  u_int n_ele = n_geometry(DIM); /// 网格中所有的单元的个数
  u_int n_side = n_geometry(1); /// 网格中的所有边的条数。
  std::vector<bool> sflag(n_side, false); /// 加在每个边上的一个标志
  std::vector<double> sjump(n_side); /// 用来边上的阶跃
  std::vector<double> earea(n_ele); /// 用来存储每个单元的面积

  /// 对单元做循环
  FEMSpace<double,DIM>::ElementIterator the_ele = fem_space.beginElement();
  FEMSpace<double,DIM>::ElementIterator end_ele = fem_space.endElement();
  for (u_int i = 0;the_ele != end_ele;++ the_ele, ++ i) {
    /// 取出单元的几何体
    const GeometryBM& geo = the_ele->geometry();

    /// 下面是单元的三个顶点的坐标
    const Point<DIM>& p0 = point(geometry(0, geo.vertex(0)).vertex(0));
    const Point<DIM>& p1 = point(geometry(0, geo.vertex(1)).vertex(0));
    const Point<DIM>& p2 = point(geometry(0, geo.vertex(2)).vertex(0));

    /// 手工计算此单元的面积
    earea[i] = 0.5*((p1[0] - p0[0])*(p2[1] - p0[1]) -
                    (p1[1] - p0[1])*(p2[0] - p0[0]));

    /**
     * 计算单元上解函数的梯度。由于使用分片线性逼近，解函数的梯度是一
     * 个常数，因此在下面的函数中，我们将计算的坐标点设为 p0 是没有关
     * 系的。如果不是这样的逼近方式，下面的代码是不对的。
     */
    std::vector<double> u_h_grad = u_h.gradient(p0, *the_ele);

    /**
     * 下面对三角形的三条边进行循环，分别计算每条边上的梯度法向分量。
     * 并加到相应的边上去。
     */
    for (u_int j = 0;j < 3;++ j) {
      u_int k = geo.boundary(j); /// 这是边的序号
      const GeometryBM& sid = geometry(1, j); /// 取出边的几何体

      /// 下面是边的两个端点
      const Point<DIM>& ep0 = point(geometry(0, sid.vertex(0)).vertex(0));
      const Point<DIM>& ep1 = point(geometry(0, sid.vertex(1)).vertex(0));

      /// 计算梯度的法向分量乘以边的长度
      double norm_grad = (u_h_grad[0]*(ep0[1] - ep1[1]) +
                          u_h_grad[1]*(ep1[0] - ep0[0]));
      /**
       * 如果 sflag[k] 是 false，说明对面的单元上还没有计算过， 因此我 
       * 们将 norm_grad 直接设为该边上的值； 否则说明对面的单元上已经计 
       * 算过了，我们将本单元的 norm_grad 减去，就得到了边上的梯度阶跃 
       * 了。
       */
      if (sflag[k] == false) {
        sjump[k] = norm_grad;
        sflag[k] = true; /// 将标志改为 true，通知对面单元
      } else {
        sjump[k] -= norm_grad;
        sflag[k] = false; /// 将标志恢复为 false，表面这个边不在边界上
      }
    }
  }

  /// 重新开始对单元做循环
  the_ele = fem_space.beginElement();
  for (u_int i = 0;the_ele != end_ele;++ the_ele, ++ i) {
    const GeometryBM& geo = the_ele->geometry(); /// 单元的几何体
    monitor(i) = 0.0;
    /// 将三个边上的阶跃都算在这个单元头上
    for (u_int j = 0;j < 3;++ j) {
      u_int k = geo.boundary(j);
      const GeometryBM& sid = geometry(1, j);
      if (sflag[k] == true) continue; /// 如果在区域边界上，则跳过去
      monitor(i) += sjump[k]*sjump[k];
    }
    monitor(i) /= earea[i]; /// 除以面积，得到分布函数
  }
  smoothMonitor(2); /// 对控制函数进行磨光操作

  /**
   * beta 是用来调节对比度的重要参数，当这个数越大的时候，就使得网格的
   * 尺寸对比越大。我们目前还不知道如何选择这个参数，只能靠试验和经验
   * 了。
   */
  double beta = 1.0e+02;
  for (u_int i = 0;i < n_ele;++ i) { /// 将控制函数取调和形式
    monitor(i) = 1.0/sqrt(1.0 + beta*monitor(i));
  }
}

/**
 * 由于是静态问题，我们在这里先直接通过在新的网格上直接求解来获得新的
 * 解。
 */
void EllipticEquation::updateSolution()
{
  solve();
}

/// 输出解函数做可视化用
void EllipticEquation::outputSolution()
{
  u_h.writeOpenDXData("u_h.dx");
}