基于改进波前法的曲面网格生成算法

覃先云+张见明

摘要： 为提高传统波前法（Advancing Front Method，AFM）的网格生成效率，利用多维搜索二叉树数据结构实现临近前沿和节点的快速查找，使整个网格生成的时间复杂度接近线性.针对周期曲面网格的生成，提出2种点修正算子，避免传统算法添加虚边界导致局部网格单元质量较差和虚边界计算复杂的问题.网格生成实例表明：多维搜索二叉树提高网格生成速度，引进点修正算子的波前法改善周期曲面网格质量.

关键词： 曲面网格生成； 波前法； 多维搜索二叉树； 点修正算子； 周期曲面

中图分类号： O241.82文献标志码： A

Abstract： To get better efficiency of mesh generation using Advancing Front Method（AFM）， the data structure of multidimensional binary search trees is used to implement fast query of the associative fronts and nodes， and the total meshing time complexity is nearly close to linearity. Two kinds of pointmodification operations are proposed to construct mesh on periodic surfaces， with which the problems generated by constructing additional virtual boundaries are avoided， such as poor local mesh quality and complicated calculation on virtual boundaries. The examples of mesh generation indicate that， the multidimensional binary search trees can improve mesh generation efficiency， and the AFM with pointmodification operations can improve the mesh quality of periodic surfaces.

Key words： surface mesh generation； advancing front method； multidimensional binary search tree； pointmodification operation； periodic surface

0引言

曲面网格生成是工程和科学计算及计算机图形学等领域的基础支撑技术.一方面，曲面网格可直接用于板和壳等结构的有限元分析以及作为边界元分析的实体表面网格，其网格质量直接影响数值计算精度；另一方面，曲面网格是有限元分析前处理三维实体网格生成的基础，很大程度地决定三维网格质量.另外，曲面网格在计算机图形学和地理信息系统等领域也有重要应用.[1]由此看来，研究曲面网格生成技术具有重要意义.

生成曲面网格的方法主要有Delaunay三角剖分算法[23]、波前法（Advancing Front Method， AFM）[410]和四叉树算法[11]等3种.AFM生成的网格单元质量最好，而且容易控制网格尺寸实现自适应网格生成，从而成为曲面网格生成技术的研究热点.[1，45]AFM的基本流程为：首先将待剖分曲面的边界离散为首尾相连的若干有向线段，得到初始前沿；然后在前沿集合中选择一个线段作为当前前沿，根据当前前沿插入一个新节点或连接已有节点，形成一个新单元；更新前沿，使前沿向待剖分区域内推进.循环选择当前前沿、插入节点（选择节点）、生成新单元和更新前沿操作直到前沿为空，网格生成结束.[1，67]

AFM的基本思想是依靠不断更新的前沿使网格逐渐向曲面待剖分区域内部生成.AFM不像Delaunay算法那样具有成熟的理论基础，但在程序实现上相对灵活.当然，实现高效的AFM算法需要更多的经验和技巧.为保证网格单元的质量和单元之间的正确连接关系，在新单元生成时需要进行大量的相交判断、包含判断和距离判断[1，10]，这些判断大约消耗AFM整个程序运行时间的80%以上.通过精心设计前沿管理数据结构，尽可能减少这些判断操作所花费的时间，可以大大提高计算效率.由于新单元的生成只受局部区域的影响，因此只需对单元和位于局部区域的临近前沿进行各种判断操作.该局部区域在曲面参数空间为一椭圆[1]，为便于实现快速查找，常用矩形代替椭圆[7，10].快速查找位于矩形内的前沿、节点和与矩形相交的前沿是加快整个判断操作的关键，也是实现高效AFM的一大难点.本文利用多维二叉树（简称kD树，k为维数）实现前沿和节点的快速查找，该数据结构由著名的计算机大师BENTLEY[12]提出.

实现高效AFM的另一大难点是如何提高周期性曲面网格生成质量、效率和稳定性.为利用常规的二维AFM对参数域进行剖分，通常对曲面人为地添加虚边界，使虚边界和原有曲面边界封闭整个曲面的参数空间.虚边界的引入将AFM扩展到各类封闭的周期曲面网格生成中.然而，对于复杂裁剪的周期曲面，特别是周期起始附近的曲面特征复杂（如小孔和复杂裁剪边界等），如何正确添加虚边界又成为难点.一般做法是：多次调用两线段求交算法，判断周期起始线与曲面原有边界是否相交，如果相交，就分割相应的边界，并改变原有边界的拓扑连接关系.新的边界拓扑连接关系的建立很复杂，而且容易出错，降低程序的稳定性.另外，如果在虚边界周围有孔和裁剪边界等特征，会导致虚边界附近的单元质量较差（即产生“虚边界问题”[5]），甚至不能用于有限元或边界元分析.为克服“虚边界问题”，本文在传统的AFM中引进点的修正算子：在进行各种判断操作时，首先调用点修正算子使相互判断的几何对象（点和线段等）都位于半个参数周期内，然后利用传统的相交、包含和距离判断算法.[1]

1AFM的改进

1.1快速查找临近前沿和节点

实现快速查找常用的数据结构是四叉树和ADT.[10]实践证明，运用这2种数据结构可以显著提高AFM效率.利用四叉树可以高效查找给定节点位于哪个树叶子内，其时间复杂度为O（log n），能够高效搜索给定节点的临近节点.但是，四叉树不便于实现任意给定矩形区域内节点的查询.ADT是基于二叉树的数据结构，可以提高查找效率[复杂度为O（log n）].可是，ADT规定节点插入准则，限制网格生成中任意节点的插入.本文采用基于二叉树的另一种数据结构kD树实现临近前沿和节点的快速查找.在kD树中，节点的插入和删除以及临近节点搜索和矩形区域查询操作的时间复杂度都为O（log n）[12]，而且对节点的插入没有条件限制，适合任意网格节点插入树中.由于网格生成在曲面的二维参数空间（u， v）中进行，因此只需要二维kD树（取k为2）.

在kD树中，规定方向辨别指示坐标u和v，在树的偶数层比较节点u值，而在奇数层比较节点v值.标志层数的值为方向辨别数.如果当前节点辨别数为偶数，那么所有u值小于或等于该节点u值的节点都位于该节点的左子树中，而所有u值大于该节点u值的节点都位于该节点的右子树中；如果当前节点辨别数为奇数，那么所有v值小于或等于该节点v值的节点都位于该节点的左子树中，而所有v值大于该节点v值的节点都位于该节点的右子树中.一棵完整的kD树构造包括节点插入、删除和优化，这些操作的具体实现见文献[12]和[13].图1显示节点A，B，C和D依次插入树的生成过程，对应节点的坐标分别为（0.3，0.4），（0.2，0.7），（0.6，0.6）和（0.5，0.8）.由随机生成的500个节点构造的一棵优化kD树见图2.

2改进的AFM流程

为方便控制网格生成的方式，把前沿分为3类：活动前沿、非活动前沿和拒绝前沿，分别存放在前沿管理数据结构对应的链表中.在网格生成过程中，只能从活动前沿中选择一个前沿作为当前前沿（基边）.当要求逐层生成单元时，把新生成的前沿存放在非活动前沿链表中，若活动前沿为空，则把非活动前沿链表中的值导入活动前沿链表中.拒绝前沿链表存放不能生成合适单元的基边.此外，为保证整个AFM的成功实施，对新单元的校核分为几何校核和拓扑校核等2类.几何校核是检查单元的形状质量是否满足给定值，而拓扑校核是检查单元间连接关系是否正确.综合前文分析，本文实现AFM的具体步骤如下.

（1）离散曲面边界，具体离散细节见文献[16].

（2）生成初始活动前沿.

（3）在活动前沿链表中依次选择基边，并基于该边生成三角形单元，逐渐剖分整个网格区域.

①准备新单元生成.首先在活动前沿中选择构成单元的基边，然后确定单元的尺寸，最后在参数空间确定搜索临近节点和前沿的矩形区域.[1， 7， 10]

②利用前文描述的kD树快速搜索位于矩形区域的节点和前沿.

③生成新单元：

（a）由单元尺寸确定理想节点的位置，使单元尽量在三维空间为等边三角形.需要注意的是，在确定理想节点时不仅要考虑单元尺寸，还要考虑附近前沿的分布特征[10].当前沿的特征阻止理想节点的生成时，转到（d）.

（b）判断搜索区域内是否存在节点可以生成有效的单元.搜索的节点为组成新单元的候选节点，并且依次校核这些节点的有效性.有效性校核[1]包括：保证由候选节点和组成基边的任一节点形成的新边没有跟前沿相交；保证由候选节点和基边组成形状合适的单元；保证候选节点和新边距离前沿为有效远的距离.其中，第一种称为拓扑校核，后2种统称为几何校核.如果有若干个单元满足校核，将选择质量最好的单元为新单元，然后转到④；否则转到（c）.

（c）根据（b）中的校核方法判断理想节点与基边是否组成合适的新单元.如果通过校核，那么生成新的单元后转到④，否则转到（d）.

（d）由理想点和基边确定若干个候选节点，同样用b的方法校核这些节点组成单元的有效性.如果没有单元通过校核，那么把该基边存放在拒绝前沿链表中；如果有若干个单元都满足校核，那么选择质量最好的单元为新单元，然后转到④.

④更新前沿管理数据.首先删掉不是前沿中的节点和边，然后添加新的节点和边到前沿管理数据结构中.当按层生成单元时，把新节点和边存放在非活动前沿链表中，否则存放在活动前沿链表中；当活动前沿链表为空时，把非活动前沿链表中的值添加到活动前沿链表中.如果活动前沿链表为非空转到①，否则转到⑤.

⑤如果拒绝前沿链表为空，那么网格生成结束；否则，将拒绝前沿中的值添加到活动前沿链表中，再转到①.先重复1次上面过程，然后考虑拒绝前沿链表是否为空，如不为空再到①，直到拒绝前沿链表为空为止.只对新单元进行拓扑校核，保证新边没有跟前沿相交就生成新单元.

3算例和分析

采用面对象的C++语言实现本文算法，编制C++类将kD树和前沿管理的实现进行封装.完成算例所用的PC机配置为双核CPU，主频为2.33 Hz，内存为3.25 GB.

3.1算例1：kD树加速AFM

以边长为2的正方形简单区域的网格生成为例，比较采用和不采用kD树这2种方案网格生成的速度.2种方案对不同数量的单元所消耗CPU时间见图5（a）.图5 （b）所示的网格单元为5 748个.从图5（a）可以明显看出，在不采用kD树时CPU时间消耗随单元数量增加呈平方趋势增长，但在采用kD树时CPU消耗时间随单元数量增加呈近似线性增长.采用kD树在单元数量相对较少的情况下优势不明显，但对大规模网格生成（10万个以上）可以明显提高效率.例如，在生成近60万个单元时，采用kD树只需要20 s左右，而不采用kD树需要多达130 s.

4结论

为实现高效、稳定的AFM生成曲面网格，对AFM具体实现进行改进.一方面，在AFM前沿数据管理中引进kD树以加快临近前沿和节点的查找，提高AFM网格生成速度；另一方面，在周期曲面网格生成中提出点修正算子并融入到AFM中，以避免添加虚边界导致网格质量较差和虚边界计算复杂的问题，从而改善网格质量和提高程序稳定性.

为方便控制网格生成方式和提高程序的稳定性，把前沿分为活动前沿、非活动前沿和拒绝前沿等3类，在算法中根据需求实现三者间的数据自动交互.对新单元的校核分为拓扑校核和几何校核，以保证整个区域网格的成功生成.基于这些改进，列出AFM实现的详细步骤.网格生成实例说明：kD树能极大地提高AFM网格生成效率，特别是对大规模网格生成效率更为明显；加入点修正算子的AFM可极大地改善周期曲面局部网格单元质量.

参考文献：

[1]关振群， 单菊林， 顾元宪. 基于黎曼度量的复杂参数曲面有限元网格生成方法[J]. 计算机学报， 2006， 29（10）： 18231833.

[2]熊英， 胡于进， 赵建军. 基于映射法和Delaunay方法的曲面三角网格划分算法[J]. 计算机辅助设计与图形学学报， 2002， 14（1）： 5660.

[3]孟宪海， 蔡强， 李吉刚， 等. 面向四面体网格生成的曲面Delaunay三角化算法[J]. 工程图学学报， 2006， 27（1）： 7681.

[4]黄晓东， 丁问司， 杜群贵. 基于波前法的参数曲面有限元网格生成算法[J]. 计算机辅助设计与图形学学报， 2010， 22（1）： 5259.

[5]杜群贵， 刘胜， 黄晓东. 闭曲面有限元网格生成的边界预调整方法[J]. 华南理工大学学报： 自然科学版， 2007， 35（2）：2732.

[6]LOHNER R. Automatic unstructured grid generators[J]. Finite Elements in Analysis and Design， 1997， 25（12）： 111134.

[7]LEE C K. Automatic adaptive metric advancing front triangulation over curved surfaces[J]. Eng Computations， 2000， 17（1）： 4874.

[8]梁义， 陈建军， 陈立岗， 等. 几何自适应参数曲面网格生成[J]. 计算机辅助设计与图形学学报， 2010， 22（2）： 327335.

[9]WU B， WANG S. Automatic triangulation over threedimensional parametric surface based on advancing front method[J]. Finite Elements Anal & Des， 2005， 41（910）： 892910.

[10]FRYKESTIG J. Advancing front mesh generation techniques with application to the finite element method[R]. Gothenburg： Chalmers University of Technology， 1994.

[11]LIANG X， EBEIDA M S， ZHANG Y. Guaranteedquality allquadrilateral mesh generation with feature preservation[J]. Comput Methods Appl Mech & Eng， 2010（199）： 20722082.

[12]BENTLEY J. Multidimensional binary search trees used for associative searching[J]. Commun ACM， 1975， 18（9）： 509517.

[13]BERG M D， OTFRIED C， van KREVELD M. 计算几何：算法与应用[M]. 3版. 邓俊辉，译. 北京： 清华大学出版社， 2009.

[14]GUAN Zhenqun， SHAN Hulin， ZHENG Yao， et al. An extended advancing front technique for closed surfaces mesh generation[J]. Int J Numer Methods Eng， 2008， 74（4）： 642667.

[15]郭新强. 边界面法四边形网格生成研究与应用[D]. 长沙： 湖南大学， 2011.

[16]CUILLIRE J C. A direct method for the automatic discretization of 3D parametric curves[J]. Comput Aided Des， 1997， 29（9）： 639647.

（编辑武晓英）

为实现高效、稳定的AFM生成曲面网格，对AFM具体实现进行改进.一方面，在AFM前沿数据管理中引进kD树以加快临近前沿和节点的查找，提高AFM网格生成速度；另一方面，在周期曲面网格生成中提出点修正算子并融入到AFM中，以避免添加虚边界导致网格质量较差和虚边界计算复杂的问题，从而改善网格质量和提高程序稳定性.

为方便控制网格生成方式和提高程序的稳定性，把前沿分为活动前沿、非活动前沿和拒绝前沿等3类，在算法中根据需求实现三者间的数据自动交互.对新单元的校核分为拓扑校核和几何校核，以保证整个区域网格的成功生成.基于这些改进，列出AFM实现的详细步骤.网格生成实例说明：kD树能极大地提高AFM网格生成效率，特别是对大规模网格生成效率更为明显；加入点修正算子的AFM可极大地改善周期曲面局部网格单元质量.

参考文献：

[1]关振群， 单菊林， 顾元宪. 基于黎曼度量的复杂参数曲面有限元网格生成方法[J]. 计算机学报， 2006， 29（10）： 18231833.

[2]熊英， 胡于进， 赵建军. 基于映射法和Delaunay方法的曲面三角网格划分算法[J]. 计算机辅助设计与图形学学报， 2002， 14（1）： 5660.

[3]孟宪海， 蔡强， 李吉刚， 等. 面向四面体网格生成的曲面Delaunay三角化算法[J]. 工程图学学报， 2006， 27（1）： 7681.

[4]黄晓东， 丁问司， 杜群贵. 基于波前法的参数曲面有限元网格生成算法[J]. 计算机辅助设计与图形学学报， 2010， 22（1）： 5259.

[5]杜群贵， 刘胜， 黄晓东. 闭曲面有限元网格生成的边界预调整方法[J]. 华南理工大学学报： 自然科学版， 2007， 35（2）：2732.

[6]LOHNER R. Automatic unstructured grid generators[J]. Finite Elements in Analysis and Design， 1997， 25（12）： 111134.

[7]LEE C K. Automatic adaptive metric advancing front triangulation over curved surfaces[J]. Eng Computations， 2000， 17（1）： 4874.

[8]梁义， 陈建军， 陈立岗， 等. 几何自适应参数曲面网格生成[J]. 计算机辅助设计与图形学学报， 2010， 22（2）： 327335.

[9]WU B， WANG S. Automatic triangulation over threedimensional parametric surface based on advancing front method[J]. Finite Elements Anal & Des， 2005， 41（910）： 892910.

[10]FRYKESTIG J. Advancing front mesh generation techniques with application to the finite element method[R]. Gothenburg： Chalmers University of Technology， 1994.

[11]LIANG X， EBEIDA M S， ZHANG Y. Guaranteedquality allquadrilateral mesh generation with feature preservation[J]. Comput Methods Appl Mech & Eng， 2010（199）： 20722082.

[12]BENTLEY J. Multidimensional binary search trees used for associative searching[J]. Commun ACM， 1975， 18（9）： 509517.

[13]BERG M D， OTFRIED C， van KREVELD M. 计算几何：算法与应用[M]. 3版. 邓俊辉，译. 北京： 清华大学出版社， 2009.

[14]GUAN Zhenqun， SHAN Hulin， ZHENG Yao， et al. An extended advancing front technique for closed surfaces mesh generation[J]. Int J Numer Methods Eng， 2008， 74（4）： 642667.

[15]郭新强. 边界面法四边形网格生成研究与应用[D]. 长沙： 湖南大学， 2011.

[16]CUILLIRE J C. A direct method for the automatic discretization of 3D parametric curves[J]. Comput Aided Des， 1997， 29（9）： 639647.

（编辑武晓英）

为实现高效、稳定的AFM生成曲面网格，对AFM具体实现进行改进.一方面，在AFM前沿数据管理中引进kD树以加快临近前沿和节点的查找，提高AFM网格生成速度；另一方面，在周期曲面网格生成中提出点修正算子并融入到AFM中，以避免添加虚边界导致网格质量较差和虚边界计算复杂的问题，从而改善网格质量和提高程序稳定性.

为方便控制网格生成方式和提高程序的稳定性，把前沿分为活动前沿、非活动前沿和拒绝前沿等3类，在算法中根据需求实现三者间的数据自动交互.对新单元的校核分为拓扑校核和几何校核，以保证整个区域网格的成功生成.基于这些改进，列出AFM实现的详细步骤.网格生成实例说明：kD树能极大地提高AFM网格生成效率，特别是对大规模网格生成效率更为明显；加入点修正算子的AFM可极大地改善周期曲面局部网格单元质量.

参考文献：

[1]关振群， 单菊林， 顾元宪. 基于黎曼度量的复杂参数曲面有限元网格生成方法[J]. 计算机学报， 2006， 29（10）： 18231833.

[2]熊英， 胡于进， 赵建军. 基于映射法和Delaunay方法的曲面三角网格划分算法[J]. 计算机辅助设计与图形学学报， 2002， 14（1）： 5660.

[3]孟宪海， 蔡强， 李吉刚， 等. 面向四面体网格生成的曲面Delaunay三角化算法[J]. 工程图学学报， 2006， 27（1）： 7681.

[4]黄晓东， 丁问司， 杜群贵. 基于波前法的参数曲面有限元网格生成算法[J]. 计算机辅助设计与图形学学报， 2010， 22（1）： 5259.

[5]杜群贵， 刘胜， 黄晓东. 闭曲面有限元网格生成的边界预调整方法[J]. 华南理工大学学报： 自然科学版， 2007， 35（2）：2732.

[6]LOHNER R. Automatic unstructured grid generators[J]. Finite Elements in Analysis and Design， 1997， 25（12）： 111134.

[7]LEE C K. Automatic adaptive metric advancing front triangulation over curved surfaces[J]. Eng Computations， 2000， 17（1）： 4874.

[8]梁义， 陈建军， 陈立岗， 等. 几何自适应参数曲面网格生成[J]. 计算机辅助设计与图形学学报， 2010， 22（2）： 327335.

[9]WU B， WANG S. Automatic triangulation over threedimensional parametric surface based on advancing front method[J]. Finite Elements Anal & Des， 2005， 41（910）： 892910.

[10]FRYKESTIG J. Advancing front mesh generation techniques with application to the finite element method[R]. Gothenburg： Chalmers University of Technology， 1994.

[11]LIANG X， EBEIDA M S， ZHANG Y. Guaranteedquality allquadrilateral mesh generation with feature preservation[J]. Comput Methods Appl Mech & Eng， 2010（199）： 20722082.

[12]BENTLEY J. Multidimensional binary search trees used for associative searching[J]. Commun ACM， 1975， 18（9）： 509517.

[13]BERG M D， OTFRIED C， van KREVELD M. 计算几何：算法与应用[M]. 3版. 邓俊辉，译. 北京： 清华大学出版社， 2009.

[14]GUAN Zhenqun， SHAN Hulin， ZHENG Yao， et al. An extended advancing front technique for closed surfaces mesh generation[J]. Int J Numer Methods Eng， 2008， 74（4）： 642667.

[15]郭新强. 边界面法四边形网格生成研究与应用[D]. 长沙： 湖南大学， 2011.

[16]CUILLIRE J C. A direct method for the automatic discretization of 3D parametric curves[J]. Comput Aided Des， 1997， 29（9）： 639647.

（编辑武晓英）