【gmsh源码阅读】边的网格划分流程-CSDN博客

本文链接：https://blog.csdn.net/loveoobaby/article/details/131158922

1. 边划分的作用

在有限元分析中，梁、桁架是一维单元，由线划分而来，同时线的网格划分也是面网格划分的基础。在划分面之前，需要先将其边界换分。如下图，矩形板进行网格划分之前，需要对四个边界进行划分，生成若干节点，以这些节点为基础生成三角形或四边形单元。

$p=p(u)=[x(u), y(u), z(u)], u\in [0, 1]$

2. 微分几何：曲线论

在阅读gmsh边网格划分代码之前，需了解一些微分几何的知识。微分几何中，曲线表示成参数方程:

$p=p(u)=[x(u), y(u), z(u)] , u\in[0, 1]$

一阶导矢： $\dot{p}=\frac{dp}{du}=[\frac{dx(u))}{du},\frac{dy(u))}{du}, \frac{dz(u))}{du} ]$

曲线的弧长： $s(u)=\int_{u{0}}^{u}\left | \dot{p}(u))\right |du$

我们需要了解一段曲线在微分几何中，一般用参数方程来表示，参数有取值范围，0代表曲线的起点，1代表终点。曲线可以求导获得n阶导矢，对一阶导矢的模进行积分可以求得弧长。了解这几个概念基本可以读懂gmsh的代码。

3. gmsh曲线划分源码

gmsh划分曲线的代码在文件meshGEdge.cpp中实现的，其入口是void meshGEdge::operator()(GEdge *ge)函数。

void meshGEdge::operator()(GEdge *ge)
{

  ge->model()->setCurrentMeshEntity(ge); //设置当前mesh的对象

  if(ge->degenerate(1)) { // 如果是退化边，直接标记完成
    ge->meshStatistics.status = GEdge::DONE;
    return;
  }

  // compute bounds
  Range<double> bounds = ge->parBounds(0); // 范围边的uv坐标范围0-1
  double t_begin = bounds.low(); // 0
  double t_end = bounds.high();  // 1

  // if a BL is ending at one of the ends, then create specific points
  std::vector<MVertex *> _addBegin, _addEnd;

  // integration to get length, number of points, etc
  int N;
  std::vector<IntPoint> Points;
  double a;
  int filterMinimumN;
  meshGEdgeProcessing(ge, t_begin, t_end, N, Points, a, filterMinimumN); // 计算需要划分的点数量N

  // look if we are doing the STL triangulation
  std::vector<MVertex *> &mesh_vertices = ge->mesh_vertices;

  // 如果不存在起始点就退出
  GPoint beg_p, end_p;
  if(!ge->getBeginVertex() || !ge->getEndVertex()) {
    Msg::Warning("Skipping curve with no begin or end point");
    return;
  }
  // 获取曲线的起点、终点
  else if(ge->getBeginVertex() == ge->getEndVertex() &&
          ge->getBeginVertex()->edges().size() == 1) {
    end_p = beg_p = ge->point(t_begin);
    Msg::Debug("Meshing periodic closed curve (tag %i, %i points)", ge->tag(),
               N);
  }
  else {
    // 多数情况下曲线有两个端点
    MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
    MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
    beg_p = GPoint(v0->x(), v0->y(), v0->z());
    end_p = GPoint(v1->x(), v1->y(), v1->z());
  }

  // do not consider the first and the last vertex (those are not
  // classified on this mesh edge)
  if(N > 1) {
    const double b = a / static_cast<double>(N - 1);
    int count = 1, NUMP = 1;
    IntPoint P1, P2;
    mesh_vertices.resize(N - 2); // 内部点的数量是N-2

      // 通过插值方式插入中间节点
    while(NUMP < N - 1) {
      P1 = Points[count - 1];
      P2 = Points[count];
      const double d = (double)NUMP * b;
      if((std::abs(P2.p) >= std::abs(d)) && (std::abs(P1.p) < std::abs(d))) {
        double const dt = P2.t - P1.t;
        double const dlc = P2.lc - P1.lc;
        double const dp = P2.p - P1.p;
        double const t = P1.t + dt / dp * (d - P1.p);
        SVector3 der = ge->firstDer(t);
        const double d = norm(der);
        double lc = d / (P1.lc + dlc / dp * (d - P1.p));
        GPoint V = ge->point(t);
        // 新生成的点保存下来
        mesh_vertices[NUMP - 1] =
          new MEdgeVertex(V.x(), V.y(), V.z(), ge, t, 0, lc);
        NUMP++;
      }
      else {
        count++;
      }
    }
    mesh_vertices.resize(NUMP - 1);
  }


  if(CTX::instance()->mesh.algo2d != ALGO_2D_BAMG &&
     CTX::instance()->mesh.algo2d != ALGO_2D_QUAD_QUASI_STRUCT &&
     CTX::instance()->mesh.algo2d != ALGO_2D_PACK_PRLGRMS)
    if(_addBegin.empty() && _addEnd.empty())
      filterPoints(ge, filterMinimumN - 2);

  std::vector<MLine *> &lines = ge->lines;

  // 将新生成的边保存下来
  for(std::size_t i = 0; i <= mesh_vertices.size(); i++) {
    MVertex *v0 =
      (i == 0) ? ge->getBeginVertex()->mesh_vertices[0] : mesh_vertices[i - 1];

    MVertex *v1 = (i == mesh_vertices.size()) ?
                    ge->getEndVertex()->mesh_vertices[0] :
                    mesh_vertices[i];
    lines.push_back(new MLine(v0, v1));
  }

  if(ge->getBeginVertex() == ge->getEndVertex() &&
     ge->getBeginVertex()->edges().size() == 1) {
    MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
    v0->x() = beg_p.x();
    v0->y() = beg_p.y();
    v0->z() = beg_p.z();
  }

  Msg::Info("Meshing curve %d (%s): %li interior vertices", ge->tag(),
             ge->getTypeString().c_str(), ge->mesh_vertices.size());

  // 标记划分结束
  ge->meshStatistics.status = GEdge::DONE;
}

其流程是通过积分得到线的长度以及内部节点的数量，通过插值方式生成内部节点。曲线长度及内部点的数量都是通过Integration函数进行的数值积分，其原理是递归第使用直线代替曲线，近似逼近真实曲线，并且递归深度最小是7层，即使是直线也得使用129个点进行计算。获得需要插入的点数量之后，进行线性插值生成内部节点。