曹偉娟，李 明，高曙明

（浙江大學CAD＆CG國家重點實驗室，浙江 杭州 310027）

1 問題的提出

自然界和人造物體中存在很多旋轉對稱結構。在對具有旋轉對稱結構的模型進行有限元分析時，通常只取模型的一個對稱單元來提高計算效率和分析精度［1］。但是，對同一個模型來說，對稱單元的選擇并不唯一，由不同形狀的對稱單元生成的有限元網格的質量及規(guī)模也不同。圖1所示為同一零件的兩個不同對稱單元，圖1c所示的對稱單元與圖1b所示的對稱單元相比，網格節(jié)點數(shù)減少了約36%。在實際的工程應用中，分析人員一般通過在計算機輔助設計（Computer Aided Design，CAD）軟件中手工構造分割面直接對原始模型進行分割的方式得到對稱單元。手工構造對稱單元主要有兩個困難：①難以確定最終生成的網格的質量以及規(guī)模是否最優(yōu)；②難以構造形狀復雜的分割面。目前，自動構造對稱單元方面的研究相對較少，且只能保證局部最優(yōu)，效率也較低。雖然對稱單元的形狀決定了網格可能達到的質量和規(guī)模，但是目前還無法顯式地表示二者之間的關系，很難通過直接比較形狀來確定對稱單元的好壞。另外，因為這種先分割再生成網格的方法會引入多余的固定邊界，所以難免引入多余的網格節(jié)點。

針對以上問題，本文提出一種構造旋轉對稱模型最優(yōu)對稱單元的方法，選取輸入模型的Delaunay三角網格的某個對稱單元作為初始局部網格，通過對初始局部網格執(zhí)行帶對稱約束的局部Delaunay細分算法，得到一個滿足質量要求且規(guī)模最優(yōu)的局部網格，該局部網格即輸入模型的最優(yōu)對稱單元。方法生成的局部網格可以直接作為有限元分析的輸入。與經典Delaunay細分算法相比，該方法僅對局部網格進行細分，因此效率較高。本文對所得對稱單元的最優(yōu)性進行了論證。

2 相關工作

本文相關工作包括對稱單元的自動構造和基于Delaunay細分的網格生成兩部分內容。

Suresh［2］在其關于對稱性在工程分析中的自動應用研究中提出一種自動構造對稱單元的方法，該方法利用某些啟發(fā)式規(guī)則先篩選出一組對稱單元作為候選，再從中選擇生成網格質量最好的一個。這些啟發(fā)式規(guī)則篩選對稱單元的條件主要有：①對稱邊界到凹點的距離最遠；②對稱邊界經過凹點等。這種啟發(fā)式規(guī)則只能得到局部最優(yōu)的對稱單元，而且候選對稱單元間的比較也較耗時。

由于Delaunay三角剖分具有一些很好的幾何特性且簡單易實現(xiàn)，成為目前最流行的通用的全自動網格生成方法之一［3］。Ruppert［4］提出的 Delaunay細分算法能夠生成滿足給定尺寸比的高質量三角網格，并從理論上保證網格規(guī)模不超過最優(yōu)網格的常數(shù)倍。該方法通過不斷在初始Delaunay三角網格中添加細分點來恢復邊界和提高三角形質量，直到所有三角形都滿足給定尺寸比。細分點主要添加在被侵入邊界邊（以該邊為直徑的圓包含其他節(jié)點）的中點以及狹長三角形的外接圓圓心。Ungor［5］給出另一種叫做偏置圓心的新的細分點，該細分點能夠生成節(jié)點數(shù)更少、網格規(guī)模更小的三角網格。Zeng［6］在曲面重網格化研究中提出一種保對稱的Delaunay三角化方法，通過添加一組對稱的細分點來保證網格的對稱性。目前大部分網格化方法都難以保證對稱模型網格的對稱性，主要原因有數(shù)值誤差導致節(jié)點不對稱、退化情況的Delaunay三角化導致三角形的不對稱等。雖然Zeng通過添加對稱的細分點解除了數(shù)值誤差的影響，但仍無法避免退化的Delaunay三角形引起的不對稱。

3 方法概述

對有限元分析來說，生成網格質量越高、網格數(shù)量越少的對稱單元越好。理論上，如果模型是對稱的，則生成的網格也應是對稱的，最優(yōu)對稱網格的對稱單元即為模型的最優(yōu)對稱單元。但是通過構造整體模型的對稱網格來確定對稱單元顯然不合實際，因為整體劃分代價較高，而且現(xiàn)有的網格劃分方法無法保證網格的對稱性。直接分割模型區(qū)域的方法無法保證網格的最優(yōu)性，由于引入了新的邊界，勢必增加多余的網格節(jié)點，而且會使原有模型的局部特征尺寸變小，網格規(guī)模進一步增大。

鑒于Delaunay細分算法對網格質量和規(guī)模的理論保證，本文提出通過基于Delaunay細分的網格劃分方法構造最優(yōu)對稱單元，其基本思想是：通過帶對稱約束的局部Delaunay細分構造輸入模型最優(yōu)對稱網格的對稱單元，細分過程中通過對稱操作維護局部網格的Delaunay性質以及對稱邊界間的對應關系。

方法的輸入是由平面直線圖（Planar Straight－Line Graph，PSLG）表示的二維旋轉對稱區(qū)域，并假定對稱階數(shù)已知，記為n；方法的輸出為由一組三角網格表示的最優(yōu)對稱單元。最優(yōu)的含義指在給定網格質量前提下，網格規(guī)模不超過最小對稱單元網格的常數(shù)倍。評價網格質量的標準為三角形角度值的下界，角度下界越小，網格質量越差；角度下界越大，網格規(guī)模越大。給定角度下界，可以預估所得網格的規(guī)模。方法流程如圖2所示，基本過程與經典Delaunay細分算法類似，主要有三點不同：①細分的對象不是整個Delaunay三角網格，而是它的一個對稱單元，即1／n局部網格；②細分點位置的確定；③添加新點后，通過帶對稱約束的邊翻轉操作更新Delaunay網格。下一章將分別針對這三點進行詳細闡述。

4 帶對稱約束的局部Delaunay細分

4.1 初始局部網格的生成

初始網格的選取遵循以下原則：①初始網格必須是模型的一個對稱單元網格；②初始網格必須是Delaunay網格。因為經典的Delaunay細分算法以模型的約束Delaunay三角割分（Constrained Delaunay Triangulation，CDT）作為初始網格，所以本文選取該網格的對稱單元作為初始局部網格，以保證細分得到的對稱單元網格具有同樣的最優(yōu)性。模型CDT的不同對稱單元具有相同的單元數(shù)，節(jié)點數(shù)因對稱邊界的不同而稍有差異，但對最終結果的影響不大，所得對稱單元網格都是最優(yōu)的。本文選擇對稱邊數(shù)目最少的對稱單元作為初始局部網格，因為這樣的網格節(jié)點數(shù)最少，并且需要移動三角形的次數(shù)通常也較少。其他的滿足Delaunay屬性的對稱單元網格不可能具有更少的節(jié)點數(shù)，這些多余的節(jié)點位置不一定是最合適的，而Delaunay細分只能增加節(jié)點，不能刪除節(jié)點，因此難以保證結果的最優(yōu)性。

由于輸入模型是對稱的，其CDT也是對稱的，可以找到一條由Delaunay邊連接而成的路徑，這條路徑和由它繞中心點旋轉360／n°得到的對稱路徑，以及模型的邊界所界定的網格，就是初始局部網格。這條路徑的末端節(jié)點分別屬于模型的外環(huán)和中心環(huán)（或中心點）。為避免退化Delaunay三角形引起的歧義，計算相關節(jié)點在以對稱中心為原點的極坐標系下的坐標，根據邊長和節(jié)點的極坐標確定唯一的三角剖分。

生成初始局部網格的步驟如下：①生成模型的CDT；②計算各節(jié)點到其他所有節(jié)點的最短路徑；③找出中心環(huán)節(jié)點和外環(huán)節(jié)點之間所有路徑中最短的一條，以及與其對稱的另一條路徑，標記兩路徑上的邊為對稱邊界，并關聯(lián)相應的對稱邊；④移除對稱邊界外的所有網格。

為快速確定CDT中邊的對稱邊，將所有對稱節(jié)點組織為環(huán)狀結構（如圖3b），將環(huán)號和環(huán)內編號作為各點的索引，環(huán)號從內到外依次增大，環(huán)內節(jié)點按逆時針方向編號。若某節(jié)點所在環(huán)的長度為l，編號為i，則可以確定其對稱點的編號為i＋l／n。

4.2 細分點位置的確定

如果網格中存在狹長三角形，則需要通過添加細分點來改善網格質量。記狹長三角形的外接圓圓心為o，則細分點的位置有四種情況：①當o位于局部網格內部時，為o；②當o位于模型區(qū)域外部或邊界邊上時，為距其最近的邊界邊的中點；③當o位于局部網格外部，但其某個對稱點落在局部網格內部時，為該對稱點；④當o位于對稱邊界上時，為o及其位于另一對稱邊界上的對稱點。前兩種情況下細分點位置的確定與經典Delaunay細分算法相同。后兩種情況的判別和處理如下：在定位o所在的三角形時，如果o落在某條對稱邊界邊的外側，則找到與該邊關連的另一條對稱邊界邊，繞中心點將o向該對稱邊旋轉360／n°得到o′，以o′代替o作為新的細分點，從對稱邊所在的三角形重新開始定位o′所在的三角形（如圖4）；如果o落在某條對稱邊界邊上，則通過旋轉將其對稱邊所在的三角形移動到該邊上（或將該邊所在的三角形移動到其對稱邊上），并更新對稱邊界。

4.3 Delaunay網格的更新

加入細分點后分裂細分點所在的三角形，并通過邊翻轉更新Delaunay網格。如果待翻轉邊為對稱邊界邊，則通過旋轉將其對稱邊所在的三角形移動到待翻轉邊上（或將待翻轉邊所在的三角形移動到其對稱邊上），并對移動后的三角形執(zhí)行翻轉操作（如圖5）。

在細分過程中，局部網格的對稱邊界在不斷變化，最終的對稱邊界可能不是最短的。為減少最終網格節(jié)點數(shù)，在細分結束后，通過最短路徑算法重新確定對稱邊界（如圖6）。

4.4 曲線邊界的處理

對于邊界為曲線的輸入模型，有兩種處理方法：①先對曲線進行離散，以直線段取代曲線，再進行細分；②采用能夠處理曲線邊界的細分算法［7］。先離散再細分生成的網格與輸入模型邊界的擬合度取決于離散的質量，而且網格規(guī)模也會受到影響。本文目前采用第一種方法。

5 最優(yōu)性證明

假設網格的角度下界為α。令θ＝360°／n，輸入模型X 的對稱群為Cn＝｛R0，R1，…，Rn－1｝，其中Ri（i∈N，0≤i＜n）為繞中心點逆時針旋轉i×θ角度的對稱變換。記局部細分算法生成的局部網格為LM，按照以下規(guī)則，通過整體細分向輸入模型的CDT中添加細分點構造整體網格M：局部細分每添加一個細分點p，整體細分便添加n個p在Cn下的對稱點qi＝Ri（p）（i∈N，0≤i＜n）。稱由狹長三角形引入的細分點為內部細分點，由被侵入邊界邊引入的細分點為邊界細分點，狹長三角形的外接圓和以被侵入邊界邊為直徑的圓為細分點的空域圓。

以下應用［4］的框架證明M 是規(guī)模最優(yōu)的，且LM＝M／n，因此LM也是規(guī)模最優(yōu)的。給定PSLG X，可以定義X所在平面上任意一點的局部特征尺寸lfs為以該點為圓心且與X的兩個邊界實體（邊或頂點）相交的最小圓的半徑［4］，如圖7所示。該函數(shù)滿足Lipschitz條件，即對X所在平面上的任意兩點p和q，有l(wèi)fs（q）≤lfs（p）＋｜pq｜。

定理1 LM＝M／n。

證明 采用歸納法證明。

（1）記初始局部網格為LM0，初始整體網格為M0，根據初始局部網格的構造方法，有LM0＝M0／n。

（2）令第i次細分得到的局部網格為LMi，第i次整體細分得到的整體網格為Mi，如果LMi＝Mi／n，則LMi＋1＝Mi＋1／n。

令M′＝R0（LMi＋1）∪R1（LMi＋1）∪…∪Rn（LMi＋1）。首先，根據整體細分的加點規(guī)則，M′的節(jié)點與Mi＋1的節(jié)點相同。其次，M′中位于Ri（LMi＋1）內部的邊都是局部Delaunay邊，而根據局部細分中對邊界邊翻轉的處理可知，Ri（LMi＋1）之間的對稱邊界邊也是局部Delaunay邊，因此M′的所有邊都是Delaunay邊，從而M′＝Mi＋1，LMi＋1＝Mi＋1／n。

（3）當局部細分結束時，整體細分也結束，根據歸納法，定理1成立。

引理1 令p和q為某次細分時加入的n個對稱細分點中的兩個相鄰點，記p的空域圓半徑為r，如果｜pq｜＜r，則必有θ＜60°，且｜pq｜＞2rsin（θ／2）。

證明 圖8中的兩個實心圓分別為p和q的空域圓，根據整體細分的對稱性可知兩圓全等，半徑均為r，線段op和oq間的夾角為θ。根據Delaunay三角形的空外接圓性質，中心點o不可能落在兩實心圓內部，因此有β≥θ。當｜pq｜＜r時，兩圓相互包含對方圓心，｜pq｜＝2rsin（β／2）≥2rsin（θ／2）。解不等式r＞2rsin（θ／2），得θ＜60°。

引理2 存在常數(shù)CT和CS，使得以下結論成立：

（1）在尚未加入細分點時，每個輸入節(jié)點p到其最近鄰接點的距離至少為lfs（p）。

（2）當加入內部細分點p及其對稱點時，p到其最近節(jié)點的距離至少為lfs（p）／CT。

（3）當加入邊界細分點p及其對稱點時，p到其最近鄰接點的距離至少為lfs（p）／CS。

證明 根據局部特征尺寸的定義，結論（1）自然成立。根據引理1，當θ≥60°時，細分點p的最近節(jié)點永遠是引入p的狹長三角形的節(jié)點或被侵入邊界邊的端點，根據文獻［4］中的推導，引理2成立。當θ＜60°時，p的最近節(jié)點可能是其相鄰對稱點q。下面證明θ＜60°時引理也成立。

結論（2）：p為狹長三角形T的外接圓圓心，T的外接圓半徑為r。按照文獻［4］中引理2結論1的推導，當CS≥CT≥1時，有r≥lfs（p）／（1＋2CSsinα）。當｜pq｜≥r時，p的最近點為T的節(jié)點，距離為r；當｜pq｜＜r時，p的最近點為q，距離為｜pq｜。根據引理1，｜pq｜＞2rsin（θ／2）≥2sin（θ／2）lfs（p）／（1＋2CSsinα）。因此，當CT≥（1＋2CSsin α）／（2sin（θ／2））時，結論（2）成立。

結論（3）：p為被侵入邊界邊S的中點，S的長度為2r。按照文獻［4］中引理2對情況2的推導，當CS≥1＋20.5CT時，有r≥lfs（p）／（1＋20.5＋21.5CSsinα）。當｜pq｜≥r時，p的最近點為S的端點，距離為r；當｜pq｜＜r時，p的最近點為q，距離為｜pq｜。根據引理1，｜pq｜＞2rsin（θ／2）≥2sin（θ／2）lfs（p）／（1＋20.5＋21.5CSsinα）。因此，當CS≥（1＋20.5＋21.5CSsinα）／（2sin（θ／2））時，結論（3）成立。

當sinα＜20.5sin（θ／2）／2時，選擇CS＝ （1＋20.5）／（2sin（θ／2）－21.5sinα），CT＝ （1＋2CSsinα）／（2sin（θ／2）），可以同時滿足方框中的條件。

引理3 M中任意節(jié)點p到其最近鄰居q的距離不小于lfs（p）／（CS＋1）。

證明 如果p在q之后被加入或兩者同時被加入，則根據引理2有｜pq｜≥lfs（p）／CS。如果p在節(jié)點q之前加入網格，則有｜pq｜≥lfs（q）／CS，根據Lipschitiz條件，｜pq｜≥（lfs（p）－｜pq｜）／CS，因此，｜pq｜≥lfs（p）／（CS＋1）。

定理2 M是規(guī)模最優(yōu)的。

證明 根據文獻［4］的論證，如果三角形半徑邊長比下限為A的網格M中任意節(jié)點p到其最近鄰居的距離不小于lfs（p）／（CS＋1），則M 的規(guī)模不超過輸入模型Ω的任一三角形半徑邊長比下限為A 的網格規(guī)模的C倍，其中C＝O（（CS＋1）2A）。因此，M是規(guī)模最優(yōu)的。

定理3 LM是規(guī)模最優(yōu)的。

證明 由定理1和定理2可得LM是規(guī)模最優(yōu)的。

6 實現(xiàn)與結果

本文方法的實現(xiàn)基于開源程序Triangle［8］。對稱邊界上的邊與其他邊界邊一樣表示為邊界子段，利用子段的標記屬性來識別和區(qū)分中心環(huán)、外環(huán)和其他邊界邊。對稱邊界上的對應子段通過指針相互關聯(lián)。因此，方法的核心是對稱邊界的處理，適用于任何基于Delaunay細分的網格劃分算法，而不依賴于具體的細分算法，也可基于其他Delaunay細分算法來實現(xiàn)。

圖9所示為若干二維模型的整體網格和兩種對稱單元網格的比較，角度下界為30°。第一列為對整體模型應用經典網格劃分算法生成的網格；第二列為對直接分割得到的對稱單元應用經典網格劃分算法生成的網格，分割線為輸入模型的CDT邊；第三列為由本文方法基于細分構造的對稱單元網格。圖后的數(shù)字表示網格節(jié)點數(shù)。兩組對稱單元網格的規(guī)模相差不大，主要原因有：①選擇CDT邊作為分割線，使得基于細分的初始局部網格與基于直接分割的初始局部網格相同；②細分過程中對稱邊界的變化不大，因此對內部網格的影響也很小。由于二維模型形狀相對簡單，兩種對稱單元網格規(guī)模差別不明顯，但是對三維模型，其對稱單元的分割邊界可能會很復雜，這種差別將會更為顯著。

7 結束語

本文提出一種基于Delaunay細分構造旋轉對稱模型的最優(yōu)對稱單元的方法，該方法具有以下特點：①能夠從理論上保證對稱單元的最優(yōu)性；②無需對原始模型進行分割和比較，且僅對局部網格進行細分，因此效率更高；③生成的對稱單元網格可以直接作為有限元分析的輸入。

未來工作主要有以下幾點：①方法在三維情形的實現(xiàn)；②改進對曲線邊界模型的處理；③考慮采用其他單元類型的最優(yōu)對稱單元的構造方法。本文方法主要面向三角形和四面體網格。對四邊形或其他單元類型，目前并沒有保證網格規(guī)模的理論依據，因此，本文方法生成的對稱單元可能不是最優(yōu)的。未來可以從其他網格化方法入手，考慮最優(yōu)對稱單元的構造。

［1］ BOSSAVIT A.Boundary value problems with symmetry and their approximation by finite elements［J］.SIAM Journal on Applied Mathematics，1993，53（5）：1352－1380.

［2］ SURESH K，SIRPOTDAR A.Automated symmetry exploitation in engineering analysis［J］.Engineering with Computers，2006，21（4）：304－311.

［3］ GUAN Zhenqun，SONG Chao，GU Yuanxian，et al.Recent advances of research on finite element mesh generation methods［J］.Journal of Computer Aided Design ＆ Computer Graphics，2003，15（1）：1－14（in Chinese）.［關振群，宋 超，顧元憲，等.有限元網格生成方法研究的新進展［J］.計算機輔助設計與圖形學學報，2003，15（1）：1－14.］

［4］ RUPPERT J.A Delaunay refinement algorithm for quality 2－dimensional mesh generation ［J］.Journal of Algorithms，1995，18（3）：548－585.

［5］ UNGOR A.Off－centers：a new type of Steiner points for computing size－optimal quality－guaranteed Delaunay triangulations［J］.Computational Geometry：Theory and Applications，2009，42（2）：109－118.

［6］ ZENG W，SHI R，GU 段 .Global surface remeshing using symmetric delaunay triangulation in uniformization spaces［C］／／Proceedings of the 2011 8th International Symposium on Voronoi Diagrams in Science and Engineering.Washington，D.C.，USA：IEEE Computer Society，2011：160－169.

［7］ ANTHONY 段 ，F(xiàn)LATAU F A.Implementing ruppert's algorithm for generic curves in 2D ［EB／OL］.［2012－03－05］.http：／／www.imr.sandia.gov／papers／imr19／RNAnthony.pdf.

［8］ SHEWCHUK 段 .Triangle：engineering a 2Dquality mesh generator and delaunay triangulator［M］／／Applied Computational Geometry：Towards Geometric Engineering.Berlin，Germany：Springer－Verlag，1996：203－222.