基于區(qū)域分解技術(shù)的并行四面體網(wǎng)格生成算法

徐 權(quán)，崔 濤，劉青凱，曹小林

（1.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京100088；2.中國科學(xué)院計算數(shù)學(xué)與科學(xué)工程計算研究所，北京100190）

0 引 言

網(wǎng)格生成是數(shù)值模擬的前處理過程，幾十年來，串行網(wǎng)格生成的研究取得了巨大的進(jìn)展，已經(jīng)形成了很多高效的串行網(wǎng)格生成算法和軟件。隨著高性能計算的快速發(fā)展，實際應(yīng)用日趨復(fù)雜，計算規(guī)?？焖僭鲩L。對于非結(jié)構(gòu)網(wǎng)格的應(yīng)用，一些大規(guī)模數(shù)值模擬［1］在千萬億次計算機的數(shù)千上萬核上模擬的網(wǎng)格數(shù)達(dá)到上億規(guī)模甚至更多。由于單機內(nèi)存和處理能力的限制，串行網(wǎng)格生成軟件難以滿足生成如此龐大的復(fù)雜網(wǎng)格的需求。因此，研究非結(jié)構(gòu)網(wǎng)格的并行生成技術(shù)在實際應(yīng)用中有著非常重要的意義。

并行網(wǎng)格生成方法主要分為任務(wù)并行和數(shù)據(jù)并行兩種并行模式。任務(wù)并行一般是細(xì)粒度層次的，這種并行方式需要開發(fā)出串行算法中各個環(huán)節(jié)之間存在的并行性，然后利用并行開發(fā)工具重新設(shè)計算法，達(dá)到并行網(wǎng)格生成的目的［2］。但是，子任務(wù)之間往往存在相互依賴關(guān)系或者存在數(shù)據(jù)競爭，算法會產(chǎn)生大量的通信和同步開銷，從而降低算法的并行效率。數(shù)據(jù)并行是先將區(qū)域分解為多個子區(qū)域，每個子區(qū)域被映射到不同的處理器上，然后調(diào)用串行的網(wǎng)格生成算法并行地生成子網(wǎng)格［3］。由于對串行算法的復(fù)用度高，實現(xiàn)難度低，數(shù)據(jù)并行模式是當(dāng)前并行網(wǎng)格生成算法的主流。數(shù)據(jù)并行算法一般采用連續(xù)或者離散區(qū)域分解方式實現(xiàn)［4］，前者是直接分解點、線、面等元素表示的區(qū)域為多個子區(qū)域，后者是先在區(qū)域內(nèi)生成稀疏的背景網(wǎng)格，然后利用網(wǎng)格劃分工具［5，6］得到子網(wǎng)格。然而，在數(shù)據(jù)分解過程中，數(shù)據(jù)并行算法往往需要引入人工邊界，當(dāng)引入的人工邊界和初始邊界距離很近時，在其附近生成的網(wǎng)格質(zhì)量一般會比較差，從而影響并行算法的穩(wěn)定性。數(shù)據(jù)并行的網(wǎng)格生成算法主要面臨在保證子區(qū)域交界面網(wǎng)格一致的情況下，如何提高網(wǎng)格的質(zhì)量等問題［7］。

本文基于Ivanov E G［8，9］提出的并行網(wǎng)格生成算法，設(shè)計了一種針對三維復(fù)雜幾何區(qū)域的并行四面體網(wǎng)格生成算法。該算法采用分而治之的策略，將復(fù)雜的三維幾何區(qū)域分解成若干個子區(qū)域，在各個子區(qū)域上采用約束Delaunay算法分別生成四面體網(wǎng)格，同時采用迭代的技術(shù)解決子區(qū)域界面上網(wǎng)格的一致性問題。

1 算法總體框架

一般地，數(shù)據(jù)并行模式的網(wǎng)格生成算法主要由以下幾個步驟：

（1）區(qū)域分解：按照一定的策略將三維幾何區(qū)域剖分為沒有重疊的子區(qū)域。

（2）生成子區(qū)域面網(wǎng)格：在子區(qū)域的所有面上并行的生成面網(wǎng)格。

（3）子區(qū)域體網(wǎng)格并行生成：對子區(qū)域并行地生成四面體網(wǎng)格。

對于復(fù)雜區(qū)域，該算法的主要困難在于：

（1）區(qū)域分解時，子區(qū)域的交界面很難自動生成。

（2）子區(qū)域網(wǎng)格生成時，子區(qū)域交界面上的網(wǎng)格很難保證一致性和協(xié)調(diào)性。對此常用的方法是對子區(qū)域面網(wǎng)格采用約束的Delaunay算法生成網(wǎng)格，但是這樣無法保證交界面附近的體網(wǎng)格質(zhì)量。

本文提出了一種改進(jìn)的并行四面體網(wǎng)格生成算法。圖1給出了算法的總體框架。它主要包含3個模塊，首先進(jìn)行全局面網(wǎng)格生成，然后是區(qū)域分解，即將全局面網(wǎng)格剖分為若干個沒有重疊的子區(qū)域，最后采用迭代技巧對子區(qū)域并行地生成四面體網(wǎng)格。這樣能夠在保持子區(qū)域間交界面上網(wǎng)格一致性的同時，保證網(wǎng)格的質(zhì)量。下面將詳細(xì)介紹本文設(shè)計的區(qū)域分解和子網(wǎng)格并行生成算法。

圖1 算法的總體框架

2 區(qū)域分解

區(qū)域分解是并行網(wǎng)格生成的前處理過程，對整個并行網(wǎng)格生成過程有著重要的影響，直接影響到并行網(wǎng)格生成算法的整體效率和最終的網(wǎng)格質(zhì)量。下面將介紹本文的區(qū)域分解算法。

2.1 分割面的形成

區(qū)域分解的關(guān)鍵問題是確定分割面后求各個子區(qū)域的幾何描述，特別是子區(qū)域間的交界面。本文首先采用Sutherland－Hodgman［10］提出的多邊形裁剪算法計算分割面與待剖分區(qū)域的所有界面相交的邊的集合，并設(shè)計如下算法基于邊的集合求分割平面和待分割區(qū)域的交界面。

算法1：Find＿Interface

輸入：三維幾何區(qū)域三角面網(wǎng)格

輸出：子區(qū)域交界面

（1）利用Sutherland－Hodgman多邊形裁剪算法對待分割區(qū)域Ω中的每個三角面進(jìn)行裁剪，標(biāo)記所有在分割平面上的邊為true。

（2）找出所有標(biāo)記為true的邊ei，記為集合E；

（3）提取E中所有有相連關(guān)系的邊組成圖G；找出圖G中所有的簡單回路，形成多邊形集合P；

（4）對多邊形集合P，找出其最小多邊形子集p；（最小多邊形子集指子集中的多邊形覆蓋整個圖且相互不存在包含關(guān)系。）

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

2.2 整體區(qū)域分解

上文講述了單步區(qū)域分解的過程，下面從整體上給出區(qū)域分解的方法。本文將采用分而治之的策略對區(qū)域進(jìn)行分解，然后把子區(qū)域分發(fā)給每個處理器。首先按照單一方向確定分割平面，將區(qū)域分解為若干個子區(qū)域。根據(jù)求解問題的需要，如果區(qū)域還需要繼續(xù)分解，則尋找第二分解方向，然后按照這個方向再次確定分割平面，對每個子區(qū)域繼續(xù)進(jìn)行分解。最后是按第三方向進(jìn)行區(qū)域分解。

3 子區(qū)域網(wǎng)格的并行生成

對于子區(qū)域網(wǎng)格的并行生成，現(xiàn)有的算法基本都是先生成交界面網(wǎng)格，然后對每個子區(qū)域并行地生成約束Delaunay網(wǎng)格［8，9］。但是，這些算法很難同時保證人工界面上網(wǎng)格的一致性和交界面附近的網(wǎng)格質(zhì)量。本文將采用迭代技巧對子區(qū)域并行地生成四面體網(wǎng)格，這樣將可以同時實現(xiàn)上面兩個目標(biāo)。為方便描述，下面給出本文使用的記號。

計算區(qū)域記為Ω，生成的網(wǎng)格記為Th，子區(qū)域記為Ωi，其網(wǎng)格記為，紅色子區(qū)域記為，其網(wǎng)格記為，黑色子區(qū)域記為，其網(wǎng)格記為，子區(qū)域和的交記為Γij，其網(wǎng)格記為。

子網(wǎng)格并行生成的算法描述如下：（流程圖見圖1）

算法2：Mesh＿Generation＿in＿Parallel

輸入：子區(qū)域面網(wǎng)格

輸出：子區(qū)域四面體網(wǎng)格

（1）對剖分后的子區(qū)域著色 （紅／黑），使得同種顏色的區(qū)域的交為空或點；

（6）go to步 （2）。

因為串行的非結(jié)構(gòu)網(wǎng)格生成方法已經(jīng)比較成熟，且對復(fù)雜的邊界和內(nèi)部約束適應(yīng)能力強，能夠生成高質(zhì)量的四面體網(wǎng)格，本算法在對每一個子區(qū)域生成Delaunay四面體網(wǎng)格時，采用現(xiàn)有的串行四面體網(wǎng)格生成軟件－TetGen［11］。因此，該算法在保證子區(qū)域間網(wǎng)格的一致性的同時，也可以生成高質(zhì)量的四面體網(wǎng)格。圖2給出了一個機械器件四分區(qū)后的四面體網(wǎng)格結(jié)果。其中，左圖是生成的四面體網(wǎng)格，右圖是內(nèi)部網(wǎng)格結(jié)果。

圖2 機械器件四分區(qū)后四面體網(wǎng)格

4 實驗算例及結(jié)果分析

本節(jié)將通過實驗算例從不同的角度評價并行四面體網(wǎng)格生成算法的性能，并通過實驗數(shù)據(jù)指出算法的優(yōu)缺點。下面將分別從算法的收斂性，并行計算加速比，網(wǎng)格質(zhì)量三方面評估算法的性能。本節(jié)使用的算例為長45，半徑為5的圓柱體。使用的測試平臺為曙光5000并行計算機系統(tǒng)，該系統(tǒng)共有384個計算節(jié)點，每個計算節(jié)點包含8顆主頻為3.0GMHz Intel Xeon，內(nèi)存為16GB的處理器內(nèi)核。

4.1 收斂性說明

算例中將三維幾何區(qū)域分別剖分為2，4，8，16個子區(qū)域，每個子區(qū)域分配一個處理器，然后并行地生成四面體網(wǎng)格。表1給出了并行網(wǎng)格生成時的迭代次數(shù)與體積控制參數(shù) （最大單元體積），處理器數(shù)之間的關(guān)系。

通過表1可以看到，迭代次數(shù)會隨著生成四面體網(wǎng)格時的控制參數(shù)和子區(qū)域數(shù)發(fā)生變化。一般情況下，對同樣的體積控制參數(shù)，隨著處理器數(shù)的增加，迭代次數(shù)會有所上升。理論上，體積控制參數(shù)足夠小時，迭代過程會在有限步內(nèi)停止。

表1 迭代次數(shù)與體積控制參數(shù)、CPU數(shù)間的關(guān)系

4.2 并行計算效率

并行計算的效率一般用并行效率和加速比衡量。他們的定義為

其中EN是并行效率，SN為并行計算加速比，N是并行計算使用的處理器個數(shù)，T1為單處理器計算所用的時間，TN為N個處理器并行計算所用的時間。

表2 給出了不同處理器數(shù)下分別生成105，106，107量級的四面體網(wǎng)格所用的時間。

表2 并行生成105，106，107量級網(wǎng)格的時間

從表2可以看出同串行的時間相比，隨著處理器數(shù)的增加，對同樣規(guī)模的四面體網(wǎng)格，并行生成網(wǎng)格的時間大大減少了。同時對于千萬量級或更大規(guī)模的四面體網(wǎng)格，串行的網(wǎng)格生成軟件不能生成，而并行的可以在短時間內(nèi)生成。

圖3 給出了不同規(guī)模下，并行地生成四面體網(wǎng)格所用時間的加速比。

圖3 并行生成105，106，107量級網(wǎng)格的加速比

從圖3可以看出，該算法能夠取得穩(wěn)定的加速比。但是，隨著網(wǎng)格規(guī)模的增大，加速比有所減小，主要是因為當(dāng)生成的網(wǎng)格規(guī)模增大時，處理器間的迭代次數(shù)增加，以致增加了并行網(wǎng)格生成的時間。另一方面，在千萬量級單元中，當(dāng)處理器的數(shù)目增加時，加速效果不如處理器數(shù)少的情形，是因為當(dāng)處理器數(shù)增加時，處理器間的迭代次數(shù)也相應(yīng)的有所增加，從而增加了并行網(wǎng)格生成的時間。但是從總體上來說，相對于串行的網(wǎng)格生成，當(dāng)子區(qū)域數(shù)較大時，并行網(wǎng)格生成在時間上大大減少了。

4.3 網(wǎng)格質(zhì)量

本節(jié)將從外接球半徑－最小邊比和二面角兩方面給出算法生成的網(wǎng)格的質(zhì)量數(shù)據(jù)。圖4，圖5分別給出了4個CPU下并行生成四面體網(wǎng)格和串行生成四面體網(wǎng)格的數(shù)據(jù)統(tǒng)計比較。其中，圖4給出了在不同的外接球半徑－最小邊比下并行和串行的生成的四面體數(shù)所占比例的對比關(guān)系；圖5給出了在不同角度下并行和串行生成的四面體網(wǎng)格中二面角個數(shù)所占比例的對比關(guān)系。其中，串行和并行生成的四面體個數(shù)分別為1633504和1634014。

通過圖4，圖5可以看到，在相同的網(wǎng)格規(guī)模下，并行算法生成的四面體網(wǎng)格質(zhì)量與串行算法生成的四面體網(wǎng)格質(zhì)量基本一致，驗證了本算法可以并行地生成高質(zhì)量的網(wǎng)格。

5 結(jié)束語

本文面向三維復(fù)雜幾何模型，基于區(qū)域分解技術(shù)，提出了一個并行的四面體網(wǎng)格生成算法。該算法能夠自動對三維幾何區(qū)域進(jìn)行區(qū)域分解，然后子區(qū)域間并行地生成四面體網(wǎng)格。實驗結(jié)果表明，同串行的四面體網(wǎng)格生成算法相比，該算法大大減少了網(wǎng)格生成的時間；在并行計算平臺上可以生成千萬量級甚至更大規(guī)模的四面體網(wǎng)格；在子區(qū)域網(wǎng)格并行生成的同時，不僅解決了子區(qū)域間交界面上網(wǎng)格的一致性問題，還保證了網(wǎng)格的質(zhì)量。

［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網(wǎng)格生成 ［J］.浙江大學(xué)學(xué)報 （工學(xué)版），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）.［陳建軍.非結(jié)構(gòu)化網(wǎng)格生成及其并行化的若干問題研究 ［D］.杭州：浙江大學(xué)，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）.［赫恩，巴克.計算機圖形學(xué)［M］.3版.蔡士杰，宋繼強，蔡敏，等譯.北京：電子工業(yè)出版社，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.