霍吉東

Delaunay四面體剖分憑借生成網(wǎng)格的高質(zhì)量性和良好逼近性，其并行網(wǎng)格生成技術(shù)備受業(yè)界關(guān)注。以逐點插入思想的Delaunay四面體網(wǎng)格剖分串行算法為基礎(chǔ)，采用“網(wǎng)格生成串行算法+新并行策略”的方式，提出一種基于數(shù)據(jù)并行的Delaunay四面體剖分并行算法。同時在Linux+MPI平臺上實現(xiàn)上述并行算法，取得了良好的計算效率。

【關(guān)鍵詞】Delaunay三角剖分 網(wǎng)格生成 并行算法 并行策略

1 引言

隨著大型并行計算機軟硬件技術(shù)的快速發(fā)展，網(wǎng)格剖分并行技術(shù)已成為科學工程計算領(lǐng)域研究的熱點之一。Delaunay三角剖分是三維空間數(shù)值模擬階段最基本的逼近單元和3D復雜對象可視化處理中最佳離散形式，剖分得到Delaunay三角網(wǎng)格具有良好的數(shù)學特性與優(yōu)化特性。

基于逐點插入思想的Delaunay三角剖分，構(gòu)成的網(wǎng)格唯一性、網(wǎng)格質(zhì)量都較好，并且滿足Delaunay三角剖分的空圓準則，具有較高的執(zhí)行效率。而基于逐點插入的Delaunay四面體剖分內(nèi)部的并行，耦合性是制約其并行效率的主要瓶頸，例如BW并行算法中插入點的沖突問題導致處理器之間較高的通信耗時，這是決定BW并行算法高低的主要因素。Yagawa等提出的自由網(wǎng)格法（free mesh method.FMM），有效的規(guī)避了耦合性的限制，充分利用網(wǎng)格的局部特性，適合大規(guī)模并行計算、負載均衡，不過局部網(wǎng)格生成的質(zhì)量是決定剖分優(yōu)劣的關(guān)鍵因素。地球物理勘探中，野外地層塊實體斷層之間耦合性很小（如圖1所示地震層塊體顯示），并且可以通過野外放炮、檢波一系列手段獲取各個層面的數(shù)據(jù)點坐標，針對于此本文結(jié)合逐點插入算法和自由網(wǎng)格方法，提出了一種基于數(shù)據(jù)并行的Delaunay四面體剖分并行算法，此算法有效縮短了數(shù)據(jù)點同時插入時通信耗時，提高了網(wǎng)格剖分效率。

2 逐點插入Delaunay四面體剖分串行算法設(shè)計

本文提出的并行算法基于逐點插入算法，在此首先給出基于逐點插入的四面體剖分串行算法的具體實現(xiàn)過程。

2.1 數(shù)據(jù)結(jié)構(gòu)

定義三種數(shù)據(jù)結(jié)構(gòu)"Point"結(jié)構(gòu)、"Tetrahedral"結(jié)構(gòu)、"Circumscribed sphere"結(jié)構(gòu)，數(shù)據(jù)結(jié)構(gòu)具體定義如下：

2.1.1 點Point

三維情況下的數(shù)據(jù)點用二維數(shù)組node[i]存儲：記錄第i個點的橫坐標、縱坐標、豎坐標。

2.1.2 四面體Tetrahedral

四面體信息存儲在四面體信息矩陣中：tera_info[i]，分別存儲第i個四面體四個頂點編號和第i個四面體四個外鄰四面體編號。

2.1.3 外接球Circumscribed sphere

三維Delaunay逐點插入算法在執(zhí)行的過程中要不斷搜索四面體外接球的信息，需要記錄下外接球的圓心坐標與半徑，用二維數(shù)組存儲：tetra_circum[]，存放第i個四面體外接球的球心橫坐標、縱坐標、豎坐標、四面體外接球半徑。

2.2 算法流程圖

三維逐點插入Delaunay四面體剖分串行算法在本文中用于子塊體網(wǎng)格剖分，具體實現(xiàn)過程如圖2所示。

3 逐點插入Delaunay四面體剖分并行算法設(shè)計

3.1 并行算法基本思想

首先對三維數(shù)據(jù)點限定在一個規(guī)則的長方體，然后將大塊體分割為多個子塊體，每一個子塊體含有上下兩層數(shù)據(jù)點（相鄰兩個子塊共享一層數(shù)據(jù)），對每一子塊體針對上下兩層數(shù)據(jù)點采用三維逐點插入Delaunay剖分串行算法進行四面體網(wǎng)格生成。然后合并子塊剖分之后得到的局部四面體網(wǎng)格，得到整體Delaunay四面體網(wǎng)格，如圖3所示。

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

將大塊體按層分解為多個子塊體，同時每一層數(shù)據(jù)點存儲在一個數(shù)據(jù)文件中，以下給出相關(guān)實現(xiàn)。

MPI_Comm_rank（MPI_COMM_WORLD，&rank）；

MPI_Comm_size（MPI_COMM_WORLD，&size）；

//為實現(xiàn)負載均衡，采用交叉數(shù)據(jù)分解方法

for（int i=rank；i

//layer是大塊體被分成子塊之后包含的總層數(shù)，proceSize啟動的進程的總數(shù)量。

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

//fileName[]中存儲的是多個數(shù)據(jù)文件名稱，一個數(shù)據(jù)文件儲存一層三維點坐標信息。

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

//讀取第i與第i+1兩相鄰層數(shù)據(jù)，將數(shù)據(jù)點信息存儲于二維數(shù)組node[][]中

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

//包含i與i+1層數(shù)據(jù)的子塊體采用前面給出的三維Delaunay四面體剖分串行算法進行網(wǎng)格剖分。

}

MPI_Finalize（）；//結(jié)束并行進程。

3.3 并行算法效率分析

實驗數(shù)據(jù)：如表所示，在野外采集七個層的4486個三維點坐標信息，并按層將其存放在格式為dig的7個文件中，同時啟動6個進程執(zhí)行?？梢钥闯鰜聿⑿兴惴ǎ却兴惴▓?zhí)行時間上明顯的進步，有較高執(zhí)行效率。

4 結(jié)論

Delaunay四面體網(wǎng)格剖分并行算法，通信耗時是限制效率的主要原因。本文基于數(shù)據(jù)并行結(jié)合逐點插入算法和自由網(wǎng)格法局部網(wǎng)格合成整體網(wǎng)格策略提出一種并行算法，此算法有效縮了通信耗時，并通過實驗驗證了并行算法的可行性與高效性。本課題只是采用了基于MPI編程編程模式的并行策略，可考慮向多混合編程模型的方向發(fā)展，可以選擇GPU作為切入點，采用MPI+OpenMP+CUDA的三級混合編程模型，充分發(fā)揮各個并行模式的優(yōu)勢。

參考文獻

[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.山東省計算機網(wǎng)絡(luò)重點實驗室 山東省濟南市 250101