基于改進波前法的曲面網格生成算法

覃先云+張見明

摘要： 為提高傳統波前法（Advancing Front Method，AFM）的網格生成效率，利用多維搜索二叉樹數據結構實現臨近前沿和節點的快速查找，使整個網格生成的時間復雜度接近線性.針對周期曲面網格的生成，提出2種點修正算子，避免傳統算法添加虛邊界導致局部網格單元質量較差和虛邊界計算復雜的問題.網格生成實例表明：多維搜索二叉樹提高網格生成速度，引進點修正算子的波前法改善周期曲面網格質量.

關鍵詞： 曲面網格生成； 波前法； 多維搜索二叉樹； 點修正算子； 周期曲面

中圖分類號： O241.82文獻標志碼： A

Abstract： To get better efficiency of mesh generation using Advancing Front Method（AFM）， the data structure of multidimensional binary search trees is used to implement fast query of the associative fronts and nodes， and the total meshing time complexity is nearly close to linearity. Two kinds of pointmodification operations are proposed to construct mesh on periodic surfaces， with which the problems generated by constructing additional virtual boundaries are avoided， such as poor local mesh quality and complicated calculation on virtual boundaries. The examples of mesh generation indicate that， the multidimensional binary search trees can improve mesh generation efficiency， and the AFM with pointmodification operations can improve the mesh quality of periodic surfaces.

Key words： surface mesh generation； advancing front method； multidimensional binary search tree； pointmodification operation； periodic surface

0引言

曲面網格生成是工程和科學計算及計算機圖形學等領域的基礎支撐技術.一方面，曲面網格可直接用于板和殼等結構的有限元分析以及作為邊界元分析的實體表面網格，其網格質量直接影響數值計算精度；另一方面，曲面網格是有限元分析前處理三維實體網格生成的基礎，很大程度地決定三維網格質量.另外，曲面網格在計算機圖形學和地理信息系統等領域也有重要應用.[1]由此看來，研究曲面網格生成技術具有重要意義.

生成曲面網格的方法主要有Delaunay三角剖分算法[23]、波前法（Advancing Front Method， AFM）[410]和四叉樹算法[11]等3種.AFM生成的網格單元質量最好，而且容易控制網格尺寸實現自適應網格生成，從而成為曲面網格生成技術的研究熱點.[1，45]AFM的基本流程為：首先將待剖分曲面的邊界離散為首尾相連的若干有向線段，得到初始前沿；然后在前沿集合中選擇一個線段作為當前前沿，根據當前前沿插入一個新節點或連接已有節點，形成一個新單元；更新前沿，使前沿向待剖分區域內推進.循環選擇當前前沿、插入節點（選擇節點）、生成新單元和更新前沿操作直到前沿為空，網格生成結束.[1，67]

AFM的基本思想是依靠不斷更新的前沿使網格逐漸向曲面待剖分區域內部生成.AFM不像Delaunay算法那樣具有成熟的理論基礎，但在程序實現上相對靈活.當然，實現高效的AFM算法需要更多的經驗和技巧.為保證網格單元的質量和單元之間的正確連接關系，在新單元生成時需要進行大量的相交判斷、包含判斷和距離判斷[1，10]，這些判斷大約消耗AFM整個程序運行時間的80%以上.通過精心設計前沿管理數據結構，盡可能減少這些判斷操作所花費的時間，可以大大提高計算效率.由于新單元的生成只受局部區域的影響，因此只需對單元和位于局部區域的臨近前沿進行各種判斷操作.該局部區域在曲面參數空間為一橢圓[1]，為便于實現快速查找，常用矩形代替橢圓[7，10].快速查找位于矩形內的前沿、節點和與矩形相交的前沿是加快整個判斷操作的關鍵，也是實現高效AFM的一大難點.本文利用多維二叉樹（簡稱kD樹，k為維數）實現前沿和節點的快速查找，該數據結構由著名的計算機大師BENTLEY[12]提出.

實現高效AFM的另一大難點是如何提高周期性曲面網格生成質量、效率和穩定性.為利用常規的二維AFM對參數域進行剖分，通常對曲面人為地添加虛邊界，使虛邊界和原有曲面邊界封閉整個曲面的參數空間.虛邊界的引入將AFM擴展到各類封閉的周期曲面網格生成中.然而，對于復雜裁剪的周期曲面，特別是周期起始附近的曲面特征復雜（如小孔和復雜裁剪邊界等），如何正確添加虛邊界又成為難點.一般做法是：多次調用兩線段求交算法，判斷周期起始線與曲面原有邊界是否相交，如果相交，就分割相應的邊界，并改變原有邊界的拓撲連接關系.新的邊界拓撲連接關系的建立很復雜，而且容易出錯，降低程序的穩定性.另外，如果在虛邊界周圍有孔和裁剪邊界等特征，會導致虛邊界附近的單元質量較差（即產生“虛邊界問題”[5]），甚至不能用于有限元或邊界元分析.為克服“虛邊界問題”，本文在傳統的AFM中引進點的修正算子：在進行各種判斷操作時，首先調用點修正算子使相互判斷的幾何對象（點和線段等）都位于半個參數周期內，然后利用傳統的相交、包含和距離判斷算法.[1]

1AFM的改進

1.1快速查找臨近前沿和節點

實現快速查找常用的數據結構是四叉樹和ADT.[10]實踐證明，運用這2種數據結構可以顯著提高AFM效率.利用四叉樹可以高效查找給定節點位于哪個樹葉子內，其時間復雜度為O（log n），能夠高效搜索給定節點的臨近節點.但是，四叉樹不便于實現任意給定矩形區域內節點的查詢.ADT是基于二叉樹的數據結構，可以提高查找效率[復雜度為O（log n）].可是，ADT規定節點插入準則，限制網格生成中任意節點的插入.本文采用基于二叉樹的另一種數據結構kD樹實現臨近前沿和節點的快速查找.在kD樹中，節點的插入和刪除以及臨近節點搜索和矩形區域查詢操作的時間復雜度都為O（log n）[12]，而且對節點的插入沒有條件限制，適合任意網格節點插入樹中.由于網格生成在曲面的二維參數空間（u， v）中進行，因此只需要二維kD樹（取k為2）.

在kD樹中，規定方向辨別指示坐標u和v，在樹的偶數層比較節點u值，而在奇數層比較節點v值.標志層數的值為方向辨別數.如果當前節點辨別數為偶數，那么所有u值小于或等于該節點u值的節點都位于該節點的左子樹中，而所有u值大于該節點u值的節點都位于該節點的右子樹中；如果當前節點辨別數為奇數，那么所有v值小于或等于該節點v值的節點都位于該節點的左子樹中，而所有v值大于該節點v值的節點都位于該節點的右子樹中.一棵完整的kD樹構造包括節點插入、刪除和優化，這些操作的具體實現見文獻[12]和[13].圖1顯示節點A，B，C和D依次插入樹的生成過程，對應節點的坐標分別為（0.3，0.4），（0.2，0.7），（0.6，0.6）和（0.5，0.8）.由隨機生成的500個節點構造的一棵優化kD樹見圖2.

2改進的AFM流程

為方便控制網格生成的方式，把前沿分為3類：活動前沿、非活動前沿和拒絕前沿，分別存放在前沿管理數據結構對應的鏈表中.在網格生成過程中，只能從活動前沿中選擇一個前沿作為當前前沿（基邊）.當要求逐層生成單元時，把新生成的前沿存放在非活動前沿鏈表中，若活動前沿為空，則把非活動前沿鏈表中的值導入活動前沿鏈表中.拒絕前沿鏈表存放不能生成合適單元的基邊.此外，為保證整個AFM的成功實施，對新單元的校核分為幾何校核和拓撲校核等2類.幾何校核是檢查單元的形狀質量是否滿足給定值，而拓撲校核是檢查單元間連接關系是否正確.綜合前文分析，本文實現AFM的具體步驟如下.

（1）離散曲面邊界，具體離散細節見文獻[16].

（2）生成初始活動前沿.

（3）在活動前沿鏈表中依次選擇基邊，并基于該邊生成三角形單元，逐漸剖分整個網格區域.

①準備新單元生成.首先在活動前沿中選擇構成單元的基邊，然后確定單元的尺寸，最后在參數空間確定搜索臨近節點和前沿的矩形區域.[1， 7， 10]

②利用前文描述的kD樹快速搜索位于矩形區域的節點和前沿.

③生成新單元：

（a）由單元尺寸確定理想節點的位置，使單元盡量在三維空間為等邊三角形.需要注意的是，在確定理想節點時不僅要考慮單元尺寸，還要考慮附近前沿的分布特征[10].當前沿的特征阻止理想節點的生成時，轉到（d）.

（b）判斷搜索區域內是否存在節點可以生成有效的單元.搜索的節點為組成新單元的候選節點，并且依次校核這些節點的有效性.有效性校核[1]包括：保證由候選節點和組成基邊的任一節點形成的新邊沒有跟前沿相交；保證由候選節點和基邊組成形狀合適的單元；保證候選節點和新邊距離前沿為有效遠的距離.其中，第一種稱為拓撲校核，后2種統稱為幾何校核.如果有若干個單元滿足校核，將選擇質量最好的單元為新單元，然后轉到④；否則轉到（c）.

（c）根據（b）中的校核方法判斷理想節點與基邊是否組成合適的新單元.如果通過校核，那么生成新的單元后轉到④，否則轉到（d）.

（d）由理想點和基邊確定若干個候選節點，同樣用b的方法校核這些節點組成單元的有效性.如果沒有單元通過校核，那么把該基邊存放在拒絕前沿鏈表中；如果有若干個單元都滿足校核，那么選擇質量最好的單元為新單元，然后轉到④.

④更新前沿管理數據.首先刪掉不是前沿中的節點和邊，然后添加新的節點和邊到前沿管理數據結構中.當按層生成單元時，把新節點和邊存放在非活動前沿鏈表中，否則存放在活動前沿鏈表中；當活動前沿鏈表為空時，把非活動前沿鏈表中的值添加到活動前沿鏈表中.如果活動前沿鏈表為非空轉到①，否則轉到⑤.

⑤如果拒絕前沿鏈表為空，那么網格生成結束；否則，將拒絕前沿中的值添加到活動前沿鏈表中，再轉到①.先重復1次上面過程，然后考慮拒絕前沿鏈表是否為空，如不為空再到①，直到拒絕前沿鏈表為空為止.只對新單元進行拓撲校核，保證新邊沒有跟前沿相交就生成新單元.

3算例和分析

采用面對象的C++語言實現本文算法，編制C++類將kD樹和前沿管理的實現進行封裝.完成算例所用的PC機配置為雙核CPU，主頻為2.33 Hz，內存為3.25 GB.

3.1算例1：kD樹加速AFM

以邊長為2的正方形簡單區域的網格生成為例，比較采用和不采用kD樹這2種方案網格生成的速度.2種方案對不同數量的單元所消耗CPU時間見圖5（a）.圖5 （b）所示的網格單元為5 748個.從圖5（a）可以明顯看出，在不采用kD樹時CPU時間消耗隨單元數量增加呈平方趨勢增長，但在采用kD樹時CPU消耗時間隨單元數量增加呈近似線性增長.采用kD樹在單元數量相對較少的情況下優勢不明顯，但對大規模網格生成（10萬個以上）可以明顯提高效率.例如，在生成近60萬個單元時，采用kD樹只需要20 s左右，而不采用kD樹需要多達130 s.

4結論

為實現高效、穩定的AFM生成曲面網格，對AFM具體實現進行改進.一方面，在AFM前沿數據管理中引進kD樹以加快臨近前沿和節點的查找，提高AFM網格生成速度；另一方面，在周期曲面網格生成中提出點修正算子并融入到AFM中，以避免添加虛邊界導致網格質量較差和虛邊界計算復雜的問題，從而改善網格質量和提高程序穩定性.

為方便控制網格生成方式和提高程序的穩定性，把前沿分為活動前沿、非活動前沿和拒絕前沿等3類，在算法中根據需求實現三者間的數據自動交互.對新單元的校核分為拓撲校核和幾何校核，以保證整個區域網格的成功生成.基于這些改進，列出AFM實現的詳細步驟.網格生成實例說明：kD樹能極大地提高AFM網格生成效率，特別是對大規模網格生成效率更為明顯；加入點修正算子的AFM可極大地改善周期曲面局部網格單元質量.

參考文獻：

[1]關振群， 單菊林， 顧元憲. 基于黎曼度量的復雜參數曲面有限元網格生成方法[J]. 計算機學報， 2006， 29（10）： 18231833.

[2]熊英， 胡于進， 趙建軍. 基于映射法和Delaunay方法的曲面三角網格劃分算法[J]. 計算機輔助設計與圖形學學報， 2002， 14（1）： 5660.

[3]孟憲海， 蔡強， 李吉剛， 等. 面向四面體網格生成的曲面Delaunay三角化算法[J]. 工程圖學學報， 2006， 27（1）： 7681.

[4]黃曉東， 丁問司， 杜群貴. 基于波前法的參數曲面有限元網格生成算法[J]. 計算機輔助設計與圖形學學報， 2010， 22（1）： 5259.

[5]杜群貴， 劉勝， 黃曉東. 閉曲面有限元網格生成的邊界預調整方法[J]. 華南理工大學學報： 自然科學版， 2007， 35（2）：2732.

[6]LOHNER R. Automatic unstructured grid generators[J]. Finite Elements in Analysis and Design， 1997， 25（12）： 111134.

[7]LEE C K. Automatic adaptive metric advancing front triangulation over curved surfaces[J]. Eng Computations， 2000， 17（1）： 4874.

[8]梁義， 陳建軍， 陳立崗， 等. 幾何自適應參數曲面網格生成[J]. 計算機輔助設計與圖形學學報， 2010， 22（2）： 327335.

[9]WU B， WANG S. Automatic triangulation over threedimensional parametric surface based on advancing front method[J]. Finite Elements Anal & Des， 2005， 41（910）： 892910.

[10]FRYKESTIG J. Advancing front mesh generation techniques with application to the finite element method[R]. Gothenburg： Chalmers University of Technology， 1994.

[11]LIANG X， EBEIDA M S， ZHANG Y. Guaranteedquality allquadrilateral mesh generation with feature preservation[J]. Comput Methods Appl Mech & Eng， 2010（199）： 20722082.

[12]BENTLEY J. Multidimensional binary search trees used for associative searching[J]. Commun ACM， 1975， 18（9）： 509517.

[13]BERG M D， OTFRIED C， van KREVELD M. 計算幾何：算法與應用[M]. 3版. 鄧俊輝，譯. 北京： 清華大學出版社， 2009.

[14]GUAN Zhenqun， SHAN Hulin， ZHENG Yao， et al. An extended advancing front technique for closed surfaces mesh generation[J]. Int J Numer Methods Eng， 2008， 74（4）： 642667.

[15]郭新強. 邊界面法四邊形網格生成研究與應用[D]. 長沙： 湖南大學， 2011.

[16]CUILLIRE J C. A direct method for the automatic discretization of 3D parametric curves[J]. Comput Aided Des， 1997， 29（9）： 639647.

（編輯武曉英）

為實現高效、穩定的AFM生成曲面網格，對AFM具體實現進行改進.一方面，在AFM前沿數據管理中引進kD樹以加快臨近前沿和節點的查找，提高AFM網格生成速度；另一方面，在周期曲面網格生成中提出點修正算子并融入到AFM中，以避免添加虛邊界導致網格質量較差和虛邊界計算復雜的問題，從而改善網格質量和提高程序穩定性.

為方便控制網格生成方式和提高程序的穩定性，把前沿分為活動前沿、非活動前沿和拒絕前沿等3類，在算法中根據需求實現三者間的數據自動交互.對新單元的校核分為拓撲校核和幾何校核，以保證整個區域網格的成功生成.基于這些改進，列出AFM實現的詳細步驟.網格生成實例說明：kD樹能極大地提高AFM網格生成效率，特別是對大規模網格生成效率更為明顯；加入點修正算子的AFM可極大地改善周期曲面局部網格單元質量.

參考文獻：

[1]關振群， 單菊林， 顧元憲. 基于黎曼度量的復雜參數曲面有限元網格生成方法[J]. 計算機學報， 2006， 29（10）： 18231833.

[2]熊英， 胡于進， 趙建軍. 基于映射法和Delaunay方法的曲面三角網格劃分算法[J]. 計算機輔助設計與圖形學學報， 2002， 14（1）： 5660.

[3]孟憲海， 蔡強， 李吉剛， 等. 面向四面體網格生成的曲面Delaunay三角化算法[J]. 工程圖學學報， 2006， 27（1）： 7681.

[4]黃曉東， 丁問司， 杜群貴. 基于波前法的參數曲面有限元網格生成算法[J]. 計算機輔助設計與圖形學學報， 2010， 22（1）： 5259.

[5]杜群貴， 劉勝， 黃曉東. 閉曲面有限元網格生成的邊界預調整方法[J]. 華南理工大學學報： 自然科學版， 2007， 35（2）：2732.

[6]LOHNER R. Automatic unstructured grid generators[J]. Finite Elements in Analysis and Design， 1997， 25（12）： 111134.

[7]LEE C K. Automatic adaptive metric advancing front triangulation over curved surfaces[J]. Eng Computations， 2000， 17（1）： 4874.

[8]梁義， 陳建軍， 陳立崗， 等. 幾何自適應參數曲面網格生成[J]. 計算機輔助設計與圖形學學報， 2010， 22（2）： 327335.

[9]WU B， WANG S. Automatic triangulation over threedimensional parametric surface based on advancing front method[J]. Finite Elements Anal & Des， 2005， 41（910）： 892910.

[10]FRYKESTIG J. Advancing front mesh generation techniques with application to the finite element method[R]. Gothenburg： Chalmers University of Technology， 1994.

[11]LIANG X， EBEIDA M S， ZHANG Y. Guaranteedquality allquadrilateral mesh generation with feature preservation[J]. Comput Methods Appl Mech & Eng， 2010（199）： 20722082.

[12]BENTLEY J. Multidimensional binary search trees used for associative searching[J]. Commun ACM， 1975， 18（9）： 509517.

[13]BERG M D， OTFRIED C， van KREVELD M. 計算幾何：算法與應用[M]. 3版. 鄧俊輝，譯. 北京： 清華大學出版社， 2009.

[14]GUAN Zhenqun， SHAN Hulin， ZHENG Yao， et al. An extended advancing front technique for closed surfaces mesh generation[J]. Int J Numer Methods Eng， 2008， 74（4）： 642667.

[15]郭新強. 邊界面法四邊形網格生成研究與應用[D]. 長沙： 湖南大學， 2011.

[16]CUILLIRE J C. A direct method for the automatic discretization of 3D parametric curves[J]. Comput Aided Des， 1997， 29（9）： 639647.

（編輯武曉英）

為實現高效、穩定的AFM生成曲面網格，對AFM具體實現進行改進.一方面，在AFM前沿數據管理中引進kD樹以加快臨近前沿和節點的查找，提高AFM網格生成速度；另一方面，在周期曲面網格生成中提出點修正算子并融入到AFM中，以避免添加虛邊界導致網格質量較差和虛邊界計算復雜的問題，從而改善網格質量和提高程序穩定性.

為方便控制網格生成方式和提高程序的穩定性，把前沿分為活動前沿、非活動前沿和拒絕前沿等3類，在算法中根據需求實現三者間的數據自動交互.對新單元的校核分為拓撲校核和幾何校核，以保證整個區域網格的成功生成.基于這些改進，列出AFM實現的詳細步驟.網格生成實例說明：kD樹能極大地提高AFM網格生成效率，特別是對大規模網格生成效率更為明顯；加入點修正算子的AFM可極大地改善周期曲面局部網格單元質量.

參考文獻：

[1]關振群， 單菊林， 顧元憲. 基于黎曼度量的復雜參數曲面有限元網格生成方法[J]. 計算機學報， 2006， 29（10）： 18231833.

[2]熊英， 胡于進， 趙建軍. 基于映射法和Delaunay方法的曲面三角網格劃分算法[J]. 計算機輔助設計與圖形學學報， 2002， 14（1）： 5660.

[3]孟憲海， 蔡強， 李吉剛， 等. 面向四面體網格生成的曲面Delaunay三角化算法[J]. 工程圖學學報， 2006， 27（1）： 7681.

[4]黃曉東， 丁問司， 杜群貴. 基于波前法的參數曲面有限元網格生成算法[J]. 計算機輔助設計與圖形學學報， 2010， 22（1）： 5259.

[5]杜群貴， 劉勝， 黃曉東. 閉曲面有限元網格生成的邊界預調整方法[J]. 華南理工大學學報： 自然科學版， 2007， 35（2）：2732.

[6]LOHNER R. Automatic unstructured grid generators[J]. Finite Elements in Analysis and Design， 1997， 25（12）： 111134.

[7]LEE C K. Automatic adaptive metric advancing front triangulation over curved surfaces[J]. Eng Computations， 2000， 17（1）： 4874.

[8]梁義， 陳建軍， 陳立崗， 等. 幾何自適應參數曲面網格生成[J]. 計算機輔助設計與圖形學學報， 2010， 22（2）： 327335.

[9]WU B， WANG S. Automatic triangulation over threedimensional parametric surface based on advancing front method[J]. Finite Elements Anal & Des， 2005， 41（910）： 892910.

[10]FRYKESTIG J. Advancing front mesh generation techniques with application to the finite element method[R]. Gothenburg： Chalmers University of Technology， 1994.

[11]LIANG X， EBEIDA M S， ZHANG Y. Guaranteedquality allquadrilateral mesh generation with feature preservation[J]. Comput Methods Appl Mech & Eng， 2010（199）： 20722082.

[12]BENTLEY J. Multidimensional binary search trees used for associative searching[J]. Commun ACM， 1975， 18（9）： 509517.

[13]BERG M D， OTFRIED C， van KREVELD M. 計算幾何：算法與應用[M]. 3版. 鄧俊輝，譯. 北京： 清華大學出版社， 2009.

[14]GUAN Zhenqun， SHAN Hulin， ZHENG Yao， et al. An extended advancing front technique for closed surfaces mesh generation[J]. Int J Numer Methods Eng， 2008， 74（4）： 642667.

[15]郭新強. 邊界面法四邊形網格生成研究與應用[D]. 長沙： 湖南大學， 2011.

[16]CUILLIRE J C. A direct method for the automatic discretization of 3D parametric curves[J]. Comput Aided Des， 1997， 29（9）： 639647.

（編輯武曉英）