平面波前法（AFT）三角剖分

在文章《AFT(波前法)三角网格剖分》中，已经对AFT方法做了较为详细的介绍。本文将在前文基础上，更进一步陈述细节，修复之代码中的一些不稳定因素。如Delaunay一样，波前法在2D和3D剖分中，思想不变。在AFT剖分曲面，仅需修改点的搜索方向，相交判定，投影等，但是相比2D来说，复杂很多。

1. 三角形数据结构

AFT方法不需要访问三角形的邻居信息，故不需要拓扑三角数据结果，只存储其三条边即可；当然也可以只存储三个顶点。在TRI中提供计算三角形质量的方法。如下图：

2.几何

几何仍然采用OCC的TopoDS_Face，此时要严格保证外边界逆时针，内孔顺时针的规则。如此才能保证正确的搜索方向。其它的和Delaunay中完全一致。

3.边界离散

将边界离散为短边，与Delaunay方法中的离散类似，只不过此时是将相邻两点连接，形成“前沿”。

4.网格划分

关于AFT的具体算法思路，文章《AFT(波前法)三角网格剖分》已经说过，此处不再赘述，仅描述关键细节。

4.1 I点和C点

AFT通过不断搜索区域的点，实现网格划分。点主要有两种，区域内部的点（I点，不存在），前沿上的点（C点，已存在）。I点为等边三角形的顶点，以当前“前沿”的中点，沿线段的法线方向偏移一定距离，形成I点。由于I点构成的三角形质量最优，所以，当I点搜索到时，不需要再搜索C点，加快代码效率。

对于C点的搜索，文章《AFT(波前法)三角网格剖分》中，直接搜索了所有前沿点。此处我们优化，只搜索距当前前沿距离小于3*len（len是网格大小）的点。并且，此点与前沿不共线，共线的三个点无法组成有效三角形。因为符合要求的C点，可能存在多个，故需判断那个最优。公式如下：

故搜索C点的代码结构如下：

if ((!isIPoint) ) {//当取I点时，三角形最优；|| (deltaOpti < 0.5) for (auto subIter = gene_front.begin(); subIter != gene_front.end(); subIter++) { //减少需要搜寻的点的个数，才是节省时间的王道；只搜索特定范围内的点；
if (iter != subIter) { // //找到一点； gp_Pnt p1 = (*subIter)->StartPoint(); //Handle(Geom_TrimmedCurve) base = *iter; //Extrema_ExtPC extrema(p1, base); Base::Point bP1((*iter)->StartPoint().X(),(*iter)->StartPoint().Y()); Base::Point bP2((*iter)->EndPoint().X(),(*iter)->EndPoint().Y());
Base::Point bP(p1.X(), p1.Y()); Base::Line L(bP1, bP2); Base::Orentation ori=L.Orien(bP);
          double distance = L.Distance(bP); if ((distance < 3.0 * len) && (ori!=Base::ON)) {//在三倍的范围内搜索新点。点与当前前沿共线，直接过滤掉，不能形成三角形。 //dis->visMesh(triangulars, gene_front, *iter, p1);
if ((p1.Distance((*iter)->EndPoint()) > 0.001) && (p1.Distance((*iter)->StartPoint()) > 0.001)) { //防止与当前前沿，共节点，点不能与当前前沿共线 //新生成三角形的两条边， Handle(Geom_TrimmedCurve) tempC1 = GC_MakeSegment((*iter)->EndPoint(), p1); Handle(Geom_TrimmedCurve) tempC2 = GC_MakeSegment(p1, (*iter)->StartPoint());
//判断这两条边是否和generation front 有交点。 TopoDS_Edge edge1 = BRepBuilderAPI_MakeEdge(tempC1); TopoDS_Edge edge2 = BRepBuilderAPI_MakeEdge((*subIter)); TopoDS_Edge edge3 = BRepBuilderAPI_MakeEdge(tempC2);
//判断搜索的方向是否正确； Handle(Geom_TrimmedCurve) tempCB = GC_MakeSegment((*iter)->StartPoint(), p1);
gp_Vec v1 = normalEdge(*iter);              gp_Vec v2((*iter)->StartPoint(), p1);
if ((!isIntercetion(*iter,tempC1, tempC2, iter,len)) && (isSameOfNormal(v1, v2))) {// //不相交 //以此点创建新的三角形； Handle(Geom_TrimmedCurve) triC1P1 = GC_MakeSegment((*iter)->StartPoint(), (*iter)->EndPoint()); tempTri.e1 = triC1P1; tempTri.e2 = tempC1; tempTri.e3 = tempC2;
if ((tempTri.detla() > deltaOpti) || (midPnt.Distance(p1) < 0.1)) { //select bigger one;如果两个值相近，优先取C点 (tempTri.detla(rC_1.getVal(), rC_2.getVal()) > deltaOpti) || (midPnt.Distance(p1) < 0.1) optiTri = tempTri; deltaOpti = optiTri.detla(); flag = 1.0;
cb_C = GC_MakeSegment(p1, (*iter)->EndPoint()); ac_C = GC_MakeSegment((*iter)->StartPoint(), p1); edge = triC1P1; finalP = p1;
isFind = true; }
              }
}
          }
} }
    }  

4.2 前沿集合更新

当生成一个三角形时，会导致前沿集合变化，必须更新，解释如下截图：

4.3 相交判断

相交判断是AFT方法的核心，新形成的三角形不与前沿相交，才是合法。相交判断，转化为三角形两条边（除当前前沿边）与前沿是否相交的问题。平面两条直线相交本身很简单，但是由于存在端点相交情况，采用下左图的公式，和下右图代码，总是会出现意外情况，故此处先采用OCC两条线相交的代码。

相交判断代码如下：

bool AFT::isIntercetion(Handle(Geom_TrimmedCurve) cTri,Handle(Geom_TrimmedCurve) cTri1, Handle(Geom_TrimmedCurve) cTri2, vector<Handle(Geom_TrimmedCurve)>::iterator iter,double len)const {  //暂且将所有线段转化为2D，所有坐标点的Z为0；  //此处不相交的是已经生成的三角形  clock_t t = clock();  bool flag = false;
  Handle(Geom2d_TrimmedCurve) c2ds[2] = { curve3dTo2d(cTri1) ,curve3dTo2d(cTri2) };  //此处错误，应该是与Generation front 取求相交；  for (int i = 0; i < 2; i++) {
    Base::Point pB1(c2ds[i]->StartPoint().X(), c2ds[i]->StartPoint().Y());    Base::Point pB2(c2ds[i]->EndPoint().X(), c2ds[i]->EndPoint().Y());
    for (auto subIter = gene_front.begin(); subIter != gene_front.end(); subIter++) {

      Base::Point pB3((*subIter)->StartPoint().X(), (*subIter)->StartPoint().Y());      Base::Point pB4((*subIter)->EndPoint().X(), (*subIter)->EndPoint().Y());
      Base::Line lB1(pB1, pB2);      Base::Line lB2(pB3, pB4);
      double dis = lB1.Distance(lB2);      //Space_Function s_v(gp_Pnt(c2ds[i]->StartPoint().X(), c2ds[i]->StartPoint().Y(),0.0));      if (!isGeneFrontinTri(cTri, cTri1, cTri2, *subIter)) {//!isIn(*iter, cTri1, cTri2,*subIter)        if (dis < 2 * len) {          if (!isSame(*iter, *subIter)) {//  iter!=subIter            //此处使用OCCT求两条直线相交的方式太过复杂；            Handle(Geom2d_TrimmedCurve) c2d = curve3dTo2d(*subIter);            //此处增加判断，generation front 的点是否在新形成的三角形之内。
            Geom2dAPI_InterCurveCurve intersection(c2ds[i], c2d);            Standard_Integer nums = intersection.NbPoints();            if (nums > 0) {//              //只有一个交点              gp_Pnt2d intersectionP = intersection.Point(1);
              //此处使用以经存在的边的端点，来判断会有错误。因为，某个已存在三角形的顶点会落在新形成的边上，此时应该判定为相交；              gp_Pnt2d startP2D(c2ds[i]->StartPoint().X(), c2ds[i]->StartPoint().Y());              gp_Pnt2d endP2D(c2ds[i]->EndPoint().X(), c2ds[i]->EndPoint().Y());
              gp_Pnt2d startP2DEdge((*subIter)->StartPoint().X(), (*subIter)->StartPoint().Y());              gp_Pnt2d endP2DEdge((*subIter)->EndPoint().X(), (*subIter)->EndPoint().Y());
              if ((intersectionP.Distance(startP2D) < 0.01) || (intersectionP.Distance(endP2D) < 0.01)) {//前沿与三角形边，端点相交                //要判断此交点是否在三角形边上，                              if ((intersectionP.Distance(startP2DEdge) < 0.01) || (intersectionP.Distance(endP2DEdge) < 0.01)) {//交点是前沿端点                  flag = false;                }                else {                  //新的点落在了三角形边上，判定为相交；                  return true;                }              }              else {                //当交点，不与两个端点共节点的时候，说明相交                return true;              }
            }            else { //没有交点时，不相交              flag = false;            }
          }
        }      }      else {        return true;      }
    }  }   return flag;  }

4.4. 过滤

当搜索I点时，需要过滤掉某些离前沿线段，前沿点太近的点。这些点会形成形状不好的三角形。过滤只在搜索I点时使用，C点是已经存在的点，无法过滤。下图是是过滤和不过滤的对比，

bool AFT::isIPointOnEdgeOrClose(const gp_Pnt& p, double tollerance)const { //点是否离某个前沿很近 //遍历Generation front,缩小范围 bool flag = false; Base::Point bP(p.X(), p.Y()); for (auto iter = gene_front.begin(); iter != gene_front.end(); iter++) {    //此处只检查距p点，3*tollerance范围内的边    Base::Point bP1((*iter)->StartPoint().X(), (*iter)->StartPoint().Y()); Base::Point bP2((*iter)->EndPoint().X(), (*iter)->EndPoint().Y());
Base::Line L(bP1, bP2);
double distance = L.Distance(bP);
if (distance<tollerance/2.0) {      return true;    }  } return flag;}

4.5 注意事项

当仅采用线段相交判断三角形合法，是不够的。会出现以下状况：这三种情况均非法，需要避免。

前两种情况判断，前沿是否有一顶点在三角形内部，如有，非法。

第三种情况通过判断，某两条前沿，是否与三角形边共线。通过前沿点是否在三角形边上，且此点不是端点。代码如下：

bool AFT::isGeneFrontinTri(Handle(Geom_TrimmedCurve) cTri1, Handle(Geom_TrimmedCurve) cTri2, Handle(Geom_TrimmedCurve) cTri3, Handle(Geom_TrimmedCurve) edge)const {  //判断前沿会不会在当前三角形内部,如果Edge的两个点，有一个在三角形内部，即是在三角形内部  gp_Pnt p0 = cTri1->StartPoint();  gp_Pnt p1 = cTri2->StartPoint();  gp_Pnt p2 = cTri3->StartPoint();
  gp_Pnt pS = edge->StartPoint();  gp_Pnt pEnd = edge->EndPoint();
  Base::Point bP0(p0.X(),p0.Y());  Base::Point bP1(p1.X(), p1.Y());  Base::Point bP2(p2.X(), p2.Y());
  Base::Point bPS(pS.X(), pS.Y());  Base::Point bPEnd(pEnd.X(), pEnd.Y());
  Base::Vector BV1(bP0, bP1);  Base::Vector BV2(bP1, bP2);  Base::Vector BV3(bP2, bP0);
  if ((bPS != bP0) && (bPS != bP1) && (bPS != bP2)) {//不是端点    Base::Vector BV1Temp(bP0, bPS);    Base::Vector BV2Temp(bP1, bPS);    Base::Vector BV3Temp(bP2, bPS);
    //同时成立,点在三角形内部    if ((BV1.CrossProduct(BV1Temp) > 0.000001) && (BV2.CrossProduct(BV2Temp) > 0.000001) && (BV3.CrossProduct(BV3Temp) > 0.000001))  return true;  }
  if ((bPEnd != bP0) && (bPEnd != bP1) && (bPEnd != bP2)) {//不是端点    Base::Vector BV1Temp(bP0, bPEnd);    Base::Vector BV2Temp(bP1, bPEnd);    Base::Vector BV3Temp(bP2, bPEnd);    if ((BV1.CrossProduct(BV1Temp) > 0.000001) && (BV2.CrossProduct(BV2Temp) > 0.000001) && (BV3.CrossProduct(BV3Temp) > 0.000001))  return true;  }
  //前沿在三角形边上，但是不共节点:直接判断点是否在边上，非端点  Base::Line BLine1(bP0, bP1);  Base::Line BLine2(bP1, bP2);  Base::Line BLine3(bP2, bP0);
  if (BLine1.isOnLine(bPS)) return true;  if (BLine2.isOnLine(bPS)) return true;  if (BLine3.isOnLine(bPS)) return true;
  if (BLine1.isOnLine(bPEnd)) return true;  if (BLine2.isOnLine(bPEnd)) return true;  if (BLine3.isOnLine(bPEnd)) return true;
  return false;}

bool isOnLine(const Point& p) { //不包含端点 if (Orien(p) != ON) return false; //排除端点 if (p.Distance(_p1) < 0.000001) return false; if (p.Distance(_p2) < 0.000001) return false;
if (std::abs(p.Distance(_p1) + p.Distance(_p2) - _p1.Distance(_p2))<0.000001) return true; return false;    }

5.其他计算

5.1 点到直线的距离

上述代码中，用到很多点到直线距离计算。工程上，点到直线距离定义如下：

代码如下：

    double Distance(const Point& p, int flag = 0)const {      Vector line_v(_p1, _p2);      Vector unit = line_v.unitV();//会调用operator=；
      Vector AP(_p1, p);      double alpha = AP.DotProduct(unit);      //判断距离的正负      double sign = 0;      Orentation ori = Orien(p);      if (flag) {        if (ori == LEFT) sign = 1.0;        if (ori == RIGHT) sign = -1.0;
      }      else {        sign = 1;      }
      if (alpha < 0) {        //垂点落在了BA之外        return _p1.Distance(p)*sign;      }      else if (alpha > line_v.Magnitude()) {        //垂点落在了AB之外；        return _p2.Distance(p)*sign;      }      else {        //通过点向直线做垂线，交点落在了直线内；        Vector CP = AP - unit * alpha;
        return CP.Magnitude()*sign;      }
    }

5.2 点在三角形内部

点在三角形内部，通过点在三条有向边的左边判断。

也可以使用类似的方法判断，三角形、四边形是否按逆时针编号。

6.实例

下图是不同尺寸的划分实例，不尽如人意的是速度还是太慢（与之前的Delaunay相比）。

路漫漫其修远兮，共勉。

参考：

Finite Element Mesh Generation