基于区域分解技术的并行四面体网格生成算法

徐 权，崔 涛，刘青凯，曹小林

（1.北京应用物理与计算数学研究所，北京100088；2.中国科学院计算数学与科学工程计算研究所，北京100190）

0 引 言

网格生成是数值模拟的前处理过程，几十年来，串行网格生成的研究取得了巨大的进展，已经形成了很多高效的串行网格生成算法和软件。随着高性能计算的快速发展，实际应用日趋复杂，计算规模快速增长。对于非结构网格的应用，一些大规模数值模拟［1］在千万亿次计算机的数千上万核上模拟的网格数达到上亿规模甚至更多。由于单机内存和处理能力的限制，串行网格生成软件难以满足生成如此庞大的复杂网格的需求。因此，研究非结构网格的并行生成技术在实际应用中有着非常重要的意义。

并行网格生成方法主要分为任务并行和数据并行两种并行模式。任务并行一般是细粒度层次的，这种并行方式需要开发出串行算法中各个环节之间存在的并行性，然后利用并行开发工具重新设计算法，达到并行网格生成的目的［2］。但是，子任务之间往往存在相互依赖关系或者存在数据竞争，算法会产生大量的通信和同步开销，从而降低算法的并行效率。数据并行是先将区域分解为多个子区域，每个子区域被映射到不同的处理器上，然后调用串行的网格生成算法并行地生成子网格［3］。由于对串行算法的复用度高，实现难度低，数据并行模式是当前并行网格生成算法的主流。数据并行算法一般采用连续或者离散区域分解方式实现［4］，前者是直接分解点、线、面等元素表示的区域为多个子区域，后者是先在区域内生成稀疏的背景网格，然后利用网格划分工具［5，6］得到子网格。然而，在数据分解过程中，数据并行算法往往需要引入人工边界，当引入的人工边界和初始边界距离很近时，在其附近生成的网格质量一般会比较差，从而影响并行算法的稳定性。数据并行的网格生成算法主要面临在保证子区域交界面网格一致的情况下，如何提高网格的质量等问题［7］。

本文基于Ivanov E G［8，9］提出的并行网格生成算法，设计了一种针对三维复杂几何区域的并行四面体网格生成算法。该算法采用分而治之的策略，将复杂的三维几何区域分解成若干个子区域，在各个子区域上采用约束Delaunay算法分别生成四面体网格，同时采用迭代的技术解决子区域界面上网格的一致性问题。

1 算法总体框架

一般地，数据并行模式的网格生成算法主要由以下几个步骤：

（1）区域分解：按照一定的策略将三维几何区域剖分为没有重叠的子区域。

（2）生成子区域面网格：在子区域的所有面上并行的生成面网格。

（3）子区域体网格并行生成：对子区域并行地生成四面体网格。

对于复杂区域，该算法的主要困难在于：

（1）区域分解时，子区域的交界面很难自动生成。

（2）子区域网格生成时，子区域交界面上的网格很难保证一致性和协调性。对此常用的方法是对子区域面网格采用约束的Delaunay算法生成网格，但是这样无法保证交界面附近的体网格质量。

本文提出了一种改进的并行四面体网格生成算法。图1给出了算法的总体框架。它主要包含3个模块，首先进行全局面网格生成，然后是区域分解，即将全局面网格剖分为若干个没有重叠的子区域，最后采用迭代技巧对子区域并行地生成四面体网格。这样能够在保持子区域间交界面上网格一致性的同时，保证网格的质量。下面将详细介绍本文设计的区域分解和子网格并行生成算法。

图1 算法的总体框架

2 区域分解

区域分解是并行网格生成的前处理过程，对整个并行网格生成过程有着重要的影响，直接影响到并行网格生成算法的整体效率和最终的网格质量。下面将介绍本文的区域分解算法。

2.1 分割面的形成

区域分解的关键问题是确定分割面后求各个子区域的几何描述，特别是子区域间的交界面。本文首先采用Sutherland－Hodgman［10］提出的多边形裁剪算法计算分割面与待剖分区域的所有界面相交的边的集合，并设计如下算法基于边的集合求分割平面和待分割区域的交界面。

算法1：Find＿Interface

输入：三维几何区域三角面网格

输出：子区域交界面

（1）利用Sutherland－Hodgman多边形裁剪算法对待分割区域Ω中的每个三角面进行裁剪，标记所有在分割平面上的边为true。

（2）找出所有标记为true的边ei，记为集合E；

（3）提取E中所有有相连关系的边组成图G；找出图G中所有的简单回路，形成多边形集合P；

（4）对多边形集合P，找出其最小多边形子集p；（最小多边形子集指子集中的多边形覆盖整个图且相互不存在包含关系。）

（5）合并最小多边形子集合p中所有多边形得到图G，即为交界面的描述Γ。

2.2 整体区域分解

上文讲述了单步区域分解的过程，下面从整体上给出区域分解的方法。本文将采用分而治之的策略对区域进行分解，然后把子区域分发给每个处理器。首先按照单一方向确定分割平面，将区域分解为若干个子区域。根据求解问题的需要，如果区域还需要继续分解，则寻找第二分解方向，然后按照这个方向再次确定分割平面，对每个子区域继续进行分解。最后是按第三方向进行区域分解。

3 子区域网格的并行生成

对于子区域网格的并行生成，现有的算法基本都是先生成交界面网格，然后对每个子区域并行地生成约束Delaunay网格［8，9］。但是，这些算法很难同时保证人工界面上网格的一致性和交界面附近的网格质量。本文将采用迭代技巧对子区域并行地生成四面体网格，这样将可以同时实现上面两个目标。为方便描述，下面给出本文使用的记号。

计算区域记为Ω，生成的网格记为Th，子区域记为Ωi，其网格记为，红色子区域记为，其网格记为，黑色子区域记为，其网格记为，子区域和的交记为Γij，其网格记为。

子网格并行生成的算法描述如下：（流程图见图1）

算法2：Mesh＿Generation＿in＿Parallel

输入：子区域面网格

输出：子区域四面体网格

（1）对剖分后的子区域着色 （红／黑），使得同种颜色的区域的交为空或点；

（6）go to步 （2）。

因为串行的非结构网格生成方法已经比较成熟，且对复杂的边界和内部约束适应能力强，能够生成高质量的四面体网格，本算法在对每一个子区域生成Delaunay四面体网格时，采用现有的串行四面体网格生成软件－TetGen［11］。因此，该算法在保证子区域间网格的一致性的同时，也可以生成高质量的四面体网格。图2给出了一个机械器件四分区后的四面体网格结果。其中，左图是生成的四面体网格，右图是内部网格结果。

图2 机械器件四分区后四面体网格

4 实验算例及结果分析

本节将通过实验算例从不同的角度评价并行四面体网格生成算法的性能，并通过实验数据指出算法的优缺点。下面将分别从算法的收敛性，并行计算加速比，网格质量三方面评估算法的性能。本节使用的算例为长45，半径为5的圆柱体。使用的测试平台为曙光5000并行计算机系统，该系统共有384个计算节点，每个计算节点包含8颗主频为3.0GMHz Intel Xeon，内存为16GB的处理器内核。

4.1 收敛性说明

算例中将三维几何区域分别剖分为2，4，8，16个子区域，每个子区域分配一个处理器，然后并行地生成四面体网格。表1给出了并行网格生成时的迭代次数与体积控制参数 （最大单元体积），处理器数之间的关系。

通过表1可以看到，迭代次数会随着生成四面体网格时的控制参数和子区域数发生变化。一般情况下，对同样的体积控制参数，随着处理器数的增加，迭代次数会有所上升。理论上，体积控制参数足够小时，迭代过程会在有限步内停止。

表1 迭代次数与体积控制参数、CPU数间的关系

4.2 并行计算效率

并行计算的效率一般用并行效率和加速比衡量。他们的定义为

其中EN是并行效率，SN为并行计算加速比，N是并行计算使用的处理器个数，T1为单处理器计算所用的时间，TN为N个处理器并行计算所用的时间。

表2 给出了不同处理器数下分别生成105，106，107量级的四面体网格所用的时间。

表2 并行生成105，106，107量级网格的时间

从表2可以看出同串行的时间相比，随着处理器数的增加，对同样规模的四面体网格，并行生成网格的时间大大减少了。同时对于千万量级或更大规模的四面体网格，串行的网格生成软件不能生成，而并行的可以在短时间内生成。

图3 给出了不同规模下，并行地生成四面体网格所用时间的加速比。

图3 并行生成105，106，107量级网格的加速比

从图3可以看出，该算法能够取得稳定的加速比。但是，随着网格规模的增大，加速比有所减小，主要是因为当生成的网格规模增大时，处理器间的迭代次数增加，以致增加了并行网格生成的时间。另一方面，在千万量级单元中，当处理器的数目增加时，加速效果不如处理器数少的情形，是因为当处理器数增加时，处理器间的迭代次数也相应的有所增加，从而增加了并行网格生成的时间。但是从总体上来说，相对于串行的网格生成，当子区域数较大时，并行网格生成在时间上大大减少了。

4.3 网格质量

本节将从外接球半径－最小边比和二面角两方面给出算法生成的网格的质量数据。图4，图5分别给出了4个CPU下并行生成四面体网格和串行生成四面体网格的数据统计比较。其中，图4给出了在不同的外接球半径－最小边比下并行和串行的生成的四面体数所占比例的对比关系；图5给出了在不同角度下并行和串行生成的四面体网格中二面角个数所占比例的对比关系。其中，串行和并行生成的四面体个数分别为1633504和1634014。

通过图4，图5可以看到，在相同的网格规模下，并行算法生成的四面体网格质量与串行算法生成的四面体网格质量基本一致，验证了本算法可以并行地生成高质量的网格。

5 结束语

本文面向三维复杂几何模型，基于区域分解技术，提出了一个并行的四面体网格生成算法。该算法能够自动对三维几何区域进行区域分解，然后子区域间并行地生成四面体网格。实验结果表明，同串行的四面体网格生成算法相比，该算法大大减少了网格生成的时间；在并行计算平台上可以生成千万量级甚至更大规模的四面体网格；在子区域网格并行生成的同时，不仅解决了子区域间交界面上网格的一致性问题，还保证了网格的质量。

［1］Lohner R.A 2nd generation parallel advancing front grid generator ［C］／／Proc of the 21st International Meshing Roundtable.Berlin：Springer Berlin Heidelberg，2013：457－474.

［2］Chrisochoides N，Nave D.Parallel Delaunay mesh generation kernel［J］.International Journal for Numerical Methods in Engineering，2003，58 （3）：161－176.

［3］LIANG Yi，CHEN Jianjun，CHEN Ligang，et al.Parallel planar Delaunay mesh generation ［J］.Journal of Zhejiang University （Engineering Science），2008，42 （4）：558－564 （in Chinese）.［梁义，陈建军，陈立岗，等.并行平面Delaunay网格生成 ［J］.浙江大学学报 （工学版），2008；42 （4）：558－564.］

［4］Chrisochoides N.A survey of parallel mesh generation methods［DB／OL］.［2011－08－11］.http：／／www.andrew.cmu.edu／user／sowen／pmesh＿survey.pdf.

［5］METIS.Family of multilevel partitioning algorithms ［CP／OL］.［2012－10－15］.http：／／glaros.dtc.umn.edt／gkhome／views／metis.

［6］Parmetis.Parallel graph partitioning and fill reducing matrix ordering［CP／OL］.［2012－12－22］.http：／／glaros.dtc.umn.edu／gkhome／metis／parmetis／overview.

［7］CHEN Jianjun.Unstructured mesh generation and its parallelization ［D］.Hangzhou：Zhejiang University，2006 （in Chinese）.［陈建军.非结构化网格生成及其并行化的若干问题研究 ［D］.杭州：浙江大学，2006.］

［8］Ivanov E G，Andra H，Kudryavtsev A N.Domain decomposition approach for automatic parallel generation of tetrahedral grids ［J］.Computational Methods in Applied Mathematics，2006，6 （2）：178－193.

［9］Andrae H，Ivanov E G，Glucheshenko O，et al.Automatic parallel generation of tetrahedral grids by using a domain decomposition approach ［J］.Journal of Computational Mathematics and Mathematic Physics，2008，48 （8）：1448－1457.

［10］Hearn D，Baker M P.Computer graphics with OpenGL［M］.3rd ed.CAI Shijie，SONG Jiqiang，CAI Min，et al transl.Beijing：Publishing House of Electronics Industry，2010：273－277 （in Chinese）.［赫恩，巴克.计算机图形学［M］.3版.蔡士杰，宋继强，蔡敏，等译.北京：电子工业出版社，2010：273－277.］

［11］Si H，A quality tetrahedral mesh generator and three－dimensional Delaunay triangulator ［CP／OL］.［2011－09－08］.http：／／tetgen.berlios.de.