基于逐點插入的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