基于逐点插入的Delaunay四面体剖分并行算法研究

霍吉东

Delaunay四面体剖分凭借生成网格的高质量性和良好逼近性，其并行网格生成技术备受业界关注。以逐点插入思想的Delaunay四面体网格剖分串行算法为基础，采用“网格生成串行算法+新并行策略”的方式，提出一种基于数据并行的Delaunay四面体剖分并行算法。同时在Linux+MPI平台上实现上述并行算法，取得了良好的计算效率。

【关键词】Delaunay三角剖分 网格生成 并行算法 并行策略

1 引言

随着大型并行计算机软硬件技术的快速发展，网格剖分并行技术已成为科学工程计算领域研究的热点之一。Delaunay三角剖分是三维空间数值模拟阶段最基本的逼近单元和3D复杂对象可视化处理中最佳离散形式，剖分得到Delaunay三角网格具有良好的数学特性与优化特性。

基于逐点插入思想的Delaunay三角剖分，构成的网格唯一性、网格质量都较好，并且满足Delaunay三角剖分的空圆准则，具有较高的执行效率。而基于逐点插入的Delaunay四面体剖分内部的并行，耦合性是制约其并行效率的主要瓶颈，例如BW并行算法中插入点的冲突问题导致处理器之间较高的通信耗时，这是决定BW并行算法高低的主要因素。Yagawa等提出的自由网格法（free mesh method.FMM），有效的规避了耦合性的限制，充分利用网格的局部特性，适合大规模并行计算、负载均衡，不过局部网格生成的质量是决定剖分优劣的关键因素。地球物理勘探中，野外地层块实体断层之间耦合性很小（如图1所示地震层块体显示），并且可以通过野外放炮、检波一系列手段获取各个层面的数据点坐标，针对于此本文结合逐点插入算法和自由网格方法，提出了一种基于数据并行的Delaunay四面体剖分并行算法，此算法有效缩短了数据点同时插入时通信耗时，提高了网格剖分效率。

2 逐点插入Delaunay四面体剖分串行算法设计

本文提出的并行算法基于逐点插入算法，在此首先给出基于逐点插入的四面体剖分串行算法的具体实现过程。

2.1 数据结构

定义三种数据结构"Point"结构、"Tetrahedral"结构、"Circumscribed sphere"结构，数据结构具体定义如下：

2.1.1 点Point

三维情况下的数据点用二维数组node[i]存储：记录第i个点的横坐标、纵坐标、竖坐标。

2.1.2 四面体Tetrahedral

四面体信息存储在四面体信息矩阵中：tera_info[i]，分别存储第i个四面体四个顶点编号和第i个四面体四个外邻四面体编号。

2.1.3 外接球Circumscribed sphere

三维Delaunay逐点插入算法在执行的过程中要不断搜索四面体外接球的信息，需要记录下外接球的圆心坐标与半径，用二维数组存储：tetra_circum[]，存放第i个四面体外接球的球心横坐标、纵坐标、竖坐标、四面体外接球半径。

2.2 算法流程图

三维逐点插入Delaunay四面体剖分串行算法在本文中用于子块体网格剖分，具体实现过程如图2所示。

3 逐点插入Delaunay四面体剖分并行算法设计

3.1 并行算法基本思想

首先对三维数据点限定在一个规则的长方体，然后将大块体分割为多个子块体，每一个子块体含有上下两层数据点（相邻两个子块共享一层数据），对每一子块体针对上下两层数据点采用三维逐点插入Delaunay剖分串行算法进行四面体网格生成。然后合并子块剖分之后得到的局部四面体网格，得到整体Delaunay四面体网格，如图3所示。

3.2 并行算法采用的并行策略

将大块体按层分解为多个子块体，同时每一层数据点存储在一个数据文件中，以下给出相关实现。

MPI_Comm_rank（MPI_COMM_WORLD，&rank）；

MPI_Comm_size（MPI_COMM_WORLD，&size）；

//为实现负载均衡，采用交叉数据分解方法

for（int i=rank；i

//layer是大块体被分成子块之后包含的总层数，proceSize启动的进程的总数量。

{ Delaunay_3d：：read_data（fileName[i]，node）；

//fileName[]中存储的是多个数据文件名称，一个数据文件储存一层三维点坐标信息。

Delaunay_3d：：read_data（fileName[i+1]，node）

//读取第i与第i+1两相邻层数据，将数据点信息存储于二维数组node[][]中

Delaunay_3d：：delaunay（node，node_num_new，i+1）；

//包含i与i+1层数据的子块体采用前面给出的三维Delaunay四面体剖分串行算法进行网格剖分。

}

MPI_Finalize（）；//結束并行进程。

3.3 并行算法效率分析

实验数据：如表所示，在野外采集七个层的4486个三维点坐标信息，并按层将其存放在格式为dig的7个文件中，同时启动6个进程执行。可以看出来并行算法，比串行算法执行时间上明显的进步，有较高执行效率。

4 结论

Delaunay四面体网格剖分并行算法，通信耗时是限制效率的主要原因。本文基于数据并行结合逐点插入算法和自由网格法局部网格合成整体网格策略提出一种并行算法，此算法有效缩了通信耗时，并通过实验验证了并行算法的可行性与高效性。本课题只是采用了基于MPI编程编程模式的并行策略，可考虑向多混合编程模型的方向发展，可以选择GPU作为切入点，采用MPI+OpenMP+CUDA的三级混合编程模型，充分发挥各个并行模式的优势。

参考文献

[1]Marc Vigo，Nuria Pla，Computing directional constrained Delaunay triangulations[J].Computers&Graphics，2000（24）：181-190.

[2]Brassel K.E and Reif.Procedure to generate thiessen polygons[J]. Geographical Analysis，1979（11）：289-303.

[3]Dwyer R.a faster divide and conquer algorithm for constructing Delaunay triangulations[J].Algorithmica， 1987，2（1/4）：137-151.

[4]Bowyer A.Computing Dirichlet tessellations[J].The Computer Journal， 1981，24（02）：162-166.

[5]Watson D F.Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes[J]. The Computer Journal，1981，24（02）：167-172.

[6]Chrisochoides N.Parallel mesh generation[M].Bruaset AM，Tveito A. Numerical solution of partial differential equations on parallel computers.Heidelberg： Springer，2006：237-264.

[7]Yagawa G，Yamada T.Free mesh method： a new meshless finite element method[J].Computational Mechanics， 1996，18（05）：383-386.

作者單位

1.山东省计算中心（国家超级计算济南中心） 山东省济南市 250101

2.山东省计算机网络重点实验室 山东省济南市 250101