基于逐點(diǎn)插入的Delaunay四面體剖分并行算法研究

霍吉東

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

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

1 引言

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

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

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

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

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 點(diǎn)Point

三維情況下的數(shù)據(jù)點(diǎn)用二維數(shù)組node[i]存儲(chǔ)：記錄第i個(gè)點(diǎn)的橫坐標(biāo)、縱坐標(biāo)、豎坐標(biāo)。

2.1.2 四面體Tetrahedral

四面體信息存儲(chǔ)在四面體信息矩陣中：tera_info[i]，分別存儲(chǔ)第i個(gè)四面體四個(gè)頂點(diǎn)編號(hào)和第i個(gè)四面體四個(gè)外鄰四面體編號(hào)。

2.1.3 外接球Circumscribed sphere

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

2.2 算法流程圖

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

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

3.1 并行算法基本思想

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

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

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

MPI_Comm_rank（MPI_COMM_WORLD，&rank）；

MPI_Comm_size（MPI_COMM_WORLD，&size）；

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

for（int i=rank；i

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

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

//fileName[]中存儲(chǔ)的是多個(gè)數(shù)據(jù)文件名稱(chēng)，一個(gè)數(shù)據(jù)文件儲(chǔ)存一層三維點(diǎn)坐標(biāo)信息。

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

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

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

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

}

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

3.3 并行算法效率分析

實(shí)驗(yàn)數(shù)據(jù)：如表所示，在野外采集七個(gè)層的4486個(gè)三維點(diǎn)坐標(biāo)信息，并按層將其存放在格式為dig的7個(gè)文件中，同時(shí)啟動(dòng)6個(gè)進(jìn)程執(zhí)行。可以看出來(lái)并行算法，比串行算法執(zhí)行時(shí)間上明顯的進(jìn)步，有較高執(zhí)行效率。

4 結(jié)論

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

參考文獻(xiàn)

[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.山東省計(jì)算中心（國(guó)家超級(jí)計(jì)算濟(jì)南中心） 山東省濟(jì)南市 250101

2.山東省計(jì)算機(jī)網(wǎng)絡(luò)重點(diǎn)實(shí)驗(yàn)室 山東省濟(jì)南市 250101