基于區域分解技術的并行四面體網格生成算法

徐 權，崔 濤，劉青凱，曹小林

（1.北京應用物理與計算數學研究所，北京100088；2.中國科學院計算數學與科學工程計算研究所，北京100190）

0 引 言

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

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

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

1 算法總體框架

一般地，數據并行模式的網格生成算法主要由以下幾個步驟：

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

（2）生成子區域面網格：在子區域的所有面上并行的生成面網格。

（3）子區域體網格并行生成：對子區域并行地生成四面體網格。

對于復雜區域，該算法的主要困難在于：

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

（2）子區域網格生成時，子區域交界面上的網格很難保證一致性和協調性。對此常用的方法是對子區域面網格采用約束的Delaunay算法生成網格，但是這樣無法保證交界面附近的體網格質量。

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

圖1 算法的總體框架

2 區域分解

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

2.1 分割面的形成

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

算法1：Find＿Interface

輸入：三維幾何區域三角面網格

輸出：子區域交界面

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

（2）找出所有標記為true的邊ei，記為集合E；

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

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

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

2.2 整體區域分解

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

3 子區域網格的并行生成

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

計算區域記為Ω，生成的網格記為Th，子區域記為Ωi，其網格記為，紅色子區域記為，其網格記為，黑色子區域記為，其網格記為，子區域和的交記為Γij，其網格記為。

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

算法2：Mesh＿Generation＿in＿Parallel

輸入：子區域面網格

輸出：子區域四面體網格

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

（6）go to步 （2）。

因為串行的非結構網格生成方法已經比較成熟，且對復雜的邊界和內部約束適應能力強，能夠生成高質量的四面體網格，本算法在對每一個子區域生成Delaunay四面體網格時，采用現有的串行四面體網格生成軟件－TetGen［11］。因此，該算法在保證子區域間網格的一致性的同時，也可以生成高質量的四面體網格。圖2給出了一個機械器件四分區后的四面體網格結果。其中，左圖是生成的四面體網格，右圖是內部網格結果。

圖2 機械器件四分區后四面體網格

4 實驗算例及結果分析

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

4.1 收斂性說明

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

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

表1 迭代次數與體積控制參數、CPU數間的關系

4.2 并行計算效率

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

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

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

表2 并行生成105，106，107量級網格的時間

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

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

圖3 并行生成105，106，107量級網格的加速比

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

4.3 網格質量

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

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

5 結束語

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

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