基于子域分解的全六面體網(wǎng)格生成方法

汪 攀 張見(jiàn)明 韓 磊 鞠傳明

池寶濤湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙，410082

0 引言

網(wǎng)格生成是數(shù)值模擬的主要性能瓶頸，其自動(dòng)生成算法的研究一直被廣泛關(guān)注。有限元方法的成功應(yīng)用在很大程度上取決于對(duì)分析對(duì)象進(jìn)行有限元網(wǎng)格劃分的合理性。相對(duì)于四面體網(wǎng)格，六面體網(wǎng)格在計(jì)算精度、劃分?jǐn)?shù)量、抗畸變程度以及重劃分次數(shù)等方面均具有優(yōu)勢(shì)[1]，因此六面體網(wǎng)格也稱為黃金網(wǎng)格。但是，由于六面體網(wǎng)格復(fù)雜的拓?fù)潢P(guān)系以及較差的幾何自適應(yīng)能力，故針對(duì)復(fù)雜三維實(shí)體的全六面體單元網(wǎng)格全自動(dòng)生成方法，目前仍處于探索階段。目前有代表性的全六面體網(wǎng)格自動(dòng)生成方法有：超限映射法[2]、掃描法[3-4]、基于柵格法[5-7]、前沿推進(jìn)法[8-9]和多子區(qū)域法[10]。其中超限映射法的優(yōu)點(diǎn)是算法簡(jiǎn)單、速度快、單元質(zhì)量好、密度可控制，并且可與形狀優(yōu)化算法集成，因此，映射法在眾多的商業(yè)有限元分析軟件中占有重要的地位。但是，映射法一般只能直接處理簡(jiǎn)單的單連通域問(wèn)題，并且對(duì)于帶有倒圓角、小孔等小特征的三維實(shí)體，其生成的六面體網(wǎng)格質(zhì)量不好；對(duì)于復(fù)雜的多連通域問(wèn)題，通常需要首先手工將待剖分域分解成幾何形狀規(guī)則的可映射子區(qū)域，然后在每個(gè)子區(qū)域上應(yīng)用超限映射法，這一缺陷不利于大規(guī)模全自動(dòng)網(wǎng)格劃分的實(shí)現(xiàn)。為了解決這一問(wèn)題，PRICE等[11]提出用中面法將三維待剖分域分解成簡(jiǎn)單子區(qū)域，但是現(xiàn)有中面算法一般需要進(jìn)行大量幾何與代數(shù)計(jì)算，自動(dòng)化程度和幾何適應(yīng)能力也有待于提高，且對(duì)于存在多凹面的幾何體來(lái)說(shuō)中軸面生成存在困難。

針對(duì)上述問(wèn)題，本文提出了一種基于子域分解的全六面體網(wǎng)格自動(dòng)生成方法，該方法首先提取實(shí)體模型的分解特征，然后通過(guò)分解特征形成分解面，最后利用分解面將復(fù)雜實(shí)體分解為多個(gè)可映射的簡(jiǎn)單子域，并在各子域上用超限映射法生成六面體網(wǎng)格。

1 算法原理

1.1 相關(guān)定義

1.1.1 與邊相關(guān)的定義

(1)二面角：與邊相連的兩個(gè)面的夾角,設(shè)為α。

(2)凸邊：當(dāng)π/2-ε<α<π/2+ε時(shí)，則該邊為凸邊，其中ε取π/6，下同。

(3)自然邊：當(dāng)π-ε<α<π+ε時(shí)，則該邊為自然邊。

(4)凹邊：當(dāng)3π/2-ε<α<3π/2+ε時(shí)，則該邊為凹邊。

(5)反轉(zhuǎn)邊：當(dāng)3π/2+ε<α?xí)r，則該邊為反轉(zhuǎn)邊。

(6)特征邊：凹邊和反轉(zhuǎn)邊統(tǒng)稱為特征邊。

(7)虛邊：用于形成封閉的分解環(huán)時(shí)構(gòu)成出來(lái)的邊。

1.1.2 與環(huán)相關(guān)的定義

(1)特征環(huán)：由相連的特征邊構(gòu)成的環(huán)。特征環(huán)有可能是開(kāi)環(huán)，也有可能是封閉的環(huán)。

(2)分解環(huán)：由特征環(huán)和虛邊組成的連續(xù)封閉的環(huán)。用于形成分割面，對(duì)三維實(shí)體進(jìn)行分割。

1.1.3 與面相關(guān)的定義

分解面：由分解環(huán)得到的虛擬面，用于分割三維實(shí)體。

1.2 算法流程

本文所提出算法的主要步驟如下。

(1)識(shí)別實(shí)體的小特征，如圓倒角、同心圓、圓孔等，得到其等效的幾何模型。

(2)從三維實(shí)體的等效模型中提取實(shí)體的特征邊(凹邊和反轉(zhuǎn)邊)，并將其加入特征邊鏈表。

(3)遍歷特征邊鏈表，找到相互連接的特征邊，構(gòu)成一個(gè)特征環(huán)，并將其加入特征環(huán)鏈表。注意：特征環(huán)有可能包含多連邊，也有可能只有一條邊，有可能是連續(xù)封閉的環(huán)，也有可能不是封閉的環(huán)。

(4)遍歷特征環(huán)鏈表。若取到的特征環(huán)是連續(xù)封閉的環(huán)，則直接通過(guò)該特征環(huán)及其支撐面構(gòu)造出分解面；若取到的特征環(huán)不是封閉的，則需要構(gòu)造合適的虛邊，將特征環(huán)補(bǔ)成連續(xù)封閉的環(huán)，然后再構(gòu)造分解面。

(5)通過(guò)分解面將三維實(shí)體分割為兩個(gè)部分，為了保證最終生成的體網(wǎng)格連續(xù)一致，需要建立這兩部分實(shí)體與分解面之間的拓?fù)潢P(guān)系。

(6)對(duì)分割后的實(shí)體進(jìn)行網(wǎng)格劃分時(shí)，為了保證空間網(wǎng)格的一致性，要先在分解面上生成面網(wǎng)格，再映射到被其分割的實(shí)體表面，之后用超限映射生成六面體網(wǎng)格。

(7)最后將各子體的六面體網(wǎng)格數(shù)據(jù)合并，即得到整個(gè)三維實(shí)體的網(wǎng)格數(shù)據(jù)。

2 實(shí)現(xiàn)過(guò)程

2.1 處理小特征

2.1.1 識(shí)別小特征

對(duì)于很多機(jī)械產(chǎn)品，為了避免應(yīng)力集中，都會(huì)引入倒圓角設(shè)計(jì)，這會(huì)導(dǎo)致捕獲實(shí)體模型的特征邊失敗。如圖1a所示的邊a、b，其二面角α=π，此時(shí)會(huì)將a、b兩條邊鑒別為自然邊，但是，很明顯圖1a所示的模型和圖1b是等效的，此時(shí)體分解失敗，所以需對(duì)倒圓角作特殊處理。要處理倒圓角，首先要識(shí)別倒圓角面，下面給出倒圓角面的三個(gè)條件：①面包含四條邊；②有兩條圓弧邊、兩條直線邊；③兩條圓弧邊是對(duì)邊，且其對(duì)應(yīng)的圓心角相同。若給定的實(shí)體表面滿足以上三個(gè)條件，則可以判定為倒圓角面。

(a)原始模型 (b)等效模型 圖1 倒圓角的處理Fig.1 The treatment of fillet

2.1.2 模板

機(jī)械產(chǎn)品中，經(jīng)常會(huì)遇見(jiàn)帶圓倒角、同心圓柱面、圓孔等小特征的實(shí)體，這些小特征的存在阻礙了六面體網(wǎng)格的生成，或使得生成的六面體網(wǎng)格質(zhì)量不好，為了解決此問(wèn)題，本文給出了如圖2～圖4所示的模板，圖中所標(biāo)識(shí)的點(diǎn)在三維空間是一條直線，引入模板以后，可以準(zhǔn)確地識(shí)別邊的類型，從而可以有效合理地進(jìn)行體分解流程。其中，圖2所示的模板可以準(zhǔn)確地提取三維實(shí)體模型的特征邊，從而對(duì)三維實(shí)體進(jìn)行有效體分隔，圖3所示的同心圓模板和圖4所示的圓孔模板，可以顯著地提高生成六面體網(wǎng)格的質(zhì)量。

圖2 倒圓角的模板Fig.2 The templates of fillet

圖3 同心圓的模板Fig.3 The templates of concentric circles

(a)幾何模型 (b)分解模型 (c)網(wǎng)格劃分圖4 圓孔的模板Fig.4 The templates of hole

2.2 形成合理的分解面

2.2.1 提取特征邊

特征邊是三維實(shí)體上形成分解環(huán)的邊，本文方法中所指的分解邊是指凹邊和反轉(zhuǎn)邊。凹邊和反轉(zhuǎn)邊是形成分解環(huán)，最終形成分解面的基礎(chǔ)，準(zhǔn)確獲取實(shí)體模型的分解邊是進(jìn)行體分解的關(guān)鍵步驟。獲取三維實(shí)體的分解邊主要包括以下兩個(gè)步驟：

(1)識(shí)別三維實(shí)體的小特征，根據(jù)圖2～圖4所示的模板，得到其等效模型。

(2)遍歷等效模型的邊界邊，計(jì)算其對(duì)應(yīng)的二面角，并標(biāo)記其類型，將凹邊和反轉(zhuǎn)邊存在鏈表中，如圖5所示的粗線邊。

圖5 圓倒角的特征邊提取示意圖Fig.5 The feature abstraction of fillet

2.2.2 形成分解環(huán)

獲取了三維實(shí)體的特征邊以后，遍歷得到的特征邊，相互連接的特征邊組成一個(gè)分解環(huán)。當(dāng)分解環(huán)是閉環(huán)時(shí)，可直接由分解環(huán)形成分解面，然后對(duì)三維實(shí)體進(jìn)行有效分割；當(dāng)分解環(huán)是開(kāi)環(huán)時(shí)，需通過(guò)搜索形成合適的“虛擬邊”，將開(kāi)環(huán)補(bǔ)成連續(xù)封閉的分解環(huán)。封閉的分解環(huán)才能有效地對(duì)三維實(shí)體進(jìn)行分割。如圖6a所示的實(shí)體模型，根據(jù)上述邊的定義提取到四條特征邊(凹邊)，如圖6b所示的粗線段。這4條特征邊首尾相連，構(gòu)成一個(gè)封閉的分解環(huán)，可直接形成分解面對(duì)三維實(shí)體進(jìn)行分割。

(a) (b) 圖6 分解環(huán)形成示意圖(一)Fig.6 The generation of parting loop

從圖5所示模型的等效模型中獲取的特征邊如圖7a所示的粗線段AB，通過(guò)該分解邊形成特征環(huán)時(shí)，沒(méi)有搜索到與其相連的分解邊，此時(shí)該特征邊形成的特征環(huán)即為一條孤立的線段，此時(shí)需要構(gòu)造虛邊。具體方法如下：在特征邊的端點(diǎn)處，作與其相連的兩條邊的延長(zhǎng)線，并與面上另一條邊相交，該端點(diǎn)與交點(diǎn)即為一條虛邊，如7b所示的粗線段AC和AD，重復(fù)上述過(guò)程，直到構(gòu)建的虛邊可以和分解邊構(gòu)成一個(gè)封閉的環(huán)。特征邊沿著兩個(gè)方向構(gòu)造虛邊，可以得到兩個(gè)分解環(huán)，如圖7c所示的ADEB和圖7d所示的ACFB，此時(shí)需要選擇一個(gè)相對(duì)合適的分解環(huán)。選取的準(zhǔn)則如下：設(shè)分解環(huán)上最短邊長(zhǎng)度為L(zhǎng)min，最長(zhǎng)邊長(zhǎng)度為L(zhǎng)max，比例β=Lmin/Lmax，取β較大的分解環(huán)，β越大，分解環(huán)的形狀越趨向正方形(當(dāng)分解環(huán)是四條邊的時(shí)候)，合理地選擇分解環(huán)可以有效地避免分割后的實(shí)體出現(xiàn)狹長(zhǎng)的片體，從而不利于生成高質(zhì)量的六面體網(wǎng)格。

(a) (b)

(c) (d)圖7 分解環(huán)形成示意圖(二)Fig.7 The generation of parting loop

2.2.3 形成分解面

三維實(shí)體之間是通過(guò)面連接的，所以需要由面來(lái)對(duì)實(shí)體進(jìn)行分割。形成分解環(huán)以后，如果分解環(huán)中存在虛邊，則根據(jù)分解環(huán)構(gòu)建一個(gè)虛平面；如果分割環(huán)中不存在虛邊，則先獲取分割環(huán)在三維實(shí)體上的支撐表面(該支撐表面不一定是平面)，得到支撐表面以后，根據(jù)支撐表面和分解環(huán)構(gòu)造出虛面，此時(shí)的構(gòu)造并不是真的生成一個(gè)面，而是在支撐面的基礎(chǔ)上加了一個(gè)約束邊界，這個(gè)虛面的幾何是基于支撐面，拓?fù)涫腔诜纸猸h(huán)。分解面形成如圖8所示。

圖8 分解面形成示意圖Fig.8 The generation of parting plane

2.3 體分解

對(duì)于復(fù)雜實(shí)體，可能會(huì)形成多個(gè)分解面，如何利用這些分解面對(duì)三維實(shí)體進(jìn)行全自動(dòng)分割是本文算法的關(guān)鍵。為了簡(jiǎn)化程序?qū)崿F(xiàn)過(guò)程，本文將所有分解面存儲(chǔ)在鏈表中，然后遍歷該鏈表。對(duì)于某一分解面，將目標(biāo)域分解為兩個(gè)子域，若這兩個(gè)子域中，至少有一個(gè)可以利用超限映射法進(jìn)行六面體網(wǎng)格劃分(子域包含6個(gè)實(shí)體面或?yàn)閳D3、圖4所示的模板)，則該分解過(guò)程合法，即可以利用該分解面對(duì)目標(biāo)域進(jìn)行分解，否則將該分解過(guò)程延后至下一次循環(huán)。上述操作可以形象地描述為，從復(fù)雜的目標(biāo)域上，逐步分割出一個(gè)可用映射法進(jìn)行網(wǎng)格劃分的子域。利用上述方法進(jìn)行有序的分解，可以避免遞歸算法，從而節(jié)省了內(nèi)存開(kāi)銷，且方便程序的調(diào)試與維護(hù)。分解面將目標(biāo)域分解為兩個(gè)子域以后，需要建立兩個(gè)子域與分解面的拓?fù)潢P(guān)系。目前，大多數(shù)CAD軟件采用基于邊界表征[12]的數(shù)據(jù)結(jié)構(gòu)(B-rep)。在邊界表征的數(shù)據(jù)結(jié)構(gòu)中，實(shí)體(body)是由面(face)表示的，面(face)是由環(huán)(loop)表示的，環(huán)(loop)是由邊(edge)表示的。本文依據(jù)邊界表征的數(shù)據(jù)結(jié)構(gòu)構(gòu)造分解面與兩個(gè)子域的拓?fù)潢P(guān)系，首先由分解邊構(gòu)造分解面，然后將分解面存儲(chǔ)于子域中用于存儲(chǔ)幾何面數(shù)組FaceArray。

2.4 子域網(wǎng)格劃分

2.4.1 子域邊界離散

用超限映射法生成六面體時(shí)，要求三維實(shí)體的幾何邊界離散線段數(shù)相等，因此對(duì)于每個(gè)分割后的塊體，在每個(gè)面上其對(duì)應(yīng)的離散線段數(shù)要相等。這實(shí)際上就是一個(gè)帶約束的線性規(guī)劃問(wèn)題，通過(guò)求解下列方程，即可得到所有邊界的離散線段數(shù)。

目標(biāo)函數(shù)為

約束條件為

式中，k為進(jìn)行體分解以后，三維實(shí)體包含的邊；Li為第i條邊的實(shí)際長(zhǎng)度；Ls為自適應(yīng)的網(wǎng)格尺寸值；iFX+為面上方向與X軸正向保持一致的所有邊界的集合；iFX-為面上方向與X軸負(fù)向保持一致的所有邊界的集合；iFY+為方向與Y軸正向保持一致的所有邊界的集合；iFY-為方向與Y軸負(fù)向保持一致的所有邊界的集合。

上述問(wèn)題可以采用整數(shù)規(guī)劃方法來(lái)求解[13]。

2.4.2 超限映射法

對(duì)三維實(shí)體進(jìn)行分割以后得到的可映射子體，可以用超限映射法[13]生成六面體網(wǎng)格。對(duì)于任意的包含6個(gè)實(shí)體表面，且是單連通域的三維實(shí)體，其對(duì)應(yīng)的參數(shù)空間為[0,1]×[0,1]×[0,1]，設(shè)(ξ,η,γ)為參數(shù)空間[0,1]×[0,1]×[0,1]中任意的一點(diǎn)，則其對(duì)應(yīng)在物理空間的坐標(biāo)是由(ξ,η,γ)在6個(gè)表面上對(duì)應(yīng)的坐標(biāo)和8個(gè)頂點(diǎn)決定的[14]。對(duì)于任意的包含4條邊界，且是單連通域的面，其對(duì)應(yīng)的參數(shù)空間為[0,1]×[0,1]，設(shè)(u,v)為參數(shù)空間[0,1]×[0,1]中任意的一點(diǎn)，則其對(duì)應(yīng)在物理空間的坐標(biāo)是由(u,v)在4條邊界上對(duì)應(yīng)的坐標(biāo)和“四邊形”面的4個(gè)頂點(diǎn)決定的[13]，現(xiàn)以面的超限映射法加以說(shuō)明。設(shè)“四邊形”面的參數(shù)域如圖9a所示，則(u,v)對(duì)應(yīng)的物理域坐標(biāo)可由下式得到：

式中，f1、f2、f3、f4為實(shí)體表面上四條邊界，如圖9b所示。

(a)參數(shù)域

(b)物理域圖9 “四邊形”面的物理域及參數(shù)域Fig.9 The physical domain and parameter domain of four sides face

2.5 子域網(wǎng)格合并

在各子域上使用超限映射法進(jìn)行網(wǎng)格劃分以后，需要將子域網(wǎng)格合并，從而得到整個(gè)目標(biāo)域的網(wǎng)格。本文首先對(duì)分解面進(jìn)行網(wǎng)格劃分，并將其作為包含該分解面的兩個(gè)子域在該表面上的網(wǎng)格數(shù)據(jù);其次，對(duì)子域上其他表面進(jìn)行網(wǎng)格劃分，并將生成的表面網(wǎng)格數(shù)據(jù)加入存儲(chǔ)目標(biāo)域網(wǎng)格數(shù)據(jù)的容器中；然后，依次在各子域上采用超限映射法生成域內(nèi)的網(wǎng)格數(shù)據(jù)，并將生成的域內(nèi)網(wǎng)格數(shù)據(jù)加入目標(biāo)域網(wǎng)格數(shù)據(jù)容器中。按照上述方式得到的網(wǎng)格數(shù)據(jù)在分解面上是連續(xù)一致的，且界面處的網(wǎng)格數(shù)據(jù)只存儲(chǔ)一遍。界面上的網(wǎng)格生成示意圖見(jiàn)圖10。

圖10 界面上的面網(wǎng)格生成Fig.10 The mesh generation of interface

3 數(shù)值算例

下面給出三個(gè)機(jī)械零部件的六面體網(wǎng)格劃分實(shí)例，如圖11～圖13所示。三個(gè)實(shí)體模型不能直接利用商業(yè)軟件對(duì)其進(jìn)行六面體劃分，必須要進(jìn)行人工分割，分解為簡(jiǎn)單子域才能有效生成六面體網(wǎng)格。本文利用特征識(shí)別技術(shù)可以全自動(dòng)地進(jìn)行體分解，得到一系列簡(jiǎn)單規(guī)則可映射的子域，然后利用超限映射法生成質(zhì)量較好的結(jié)構(gòu)化六面體網(wǎng)格，實(shí)現(xiàn)自動(dòng)劃分。

(a)原始模型

(b)分解模型

(c)生成的六面體網(wǎng)格圖11 六面體網(wǎng)格生成算例一Fig.11 The first example of hexahedral mesh generation

(a)原始模型 (b)分解體

(c)進(jìn)一步分解體 (d)生成的六面體網(wǎng)格圖12 六面體網(wǎng)格生成算例二Fig.12 The second example of hexahedral mesh generation

圖13 六面體網(wǎng)格生成算例三Fig.13 The third example of hexahedral mesh generation

4 總結(jié)

(1)由于子域網(wǎng)格是由超限映射法生成的，因此本文提出的方法能夠生成質(zhì)量較好的結(jié)構(gòu)化六面體網(wǎng)格，且算法效率較高。

(2)對(duì)于含有倒角、同心圓、小孔等小特征的三維實(shí)體，本文引入的模板能夠生成質(zhì)量較好的六面體網(wǎng)格。

(3)相對(duì)于傳統(tǒng)的商業(yè)軟件，本文的算法能夠全自動(dòng)地生成六面體網(wǎng)格，避免了大量人機(jī)交互操作，提高了網(wǎng)格生成的效率。

(4)數(shù)值實(shí)驗(yàn)表明，對(duì)于無(wú)法用商業(yè)軟件直接進(jìn)六面體劃分的復(fù)雜實(shí)體，本文提出的方法能夠全自動(dòng)地進(jìn)行六面體網(wǎng)格劃分。

[1] SCHONNING A, OOMMEN B, IONESCU I, et al. Hexahedral Mesh Development of Free-formed Geometry: the Human Femur Exemplified[J]. Computer-aided Design,2009,41(8):566-572.

[2] GORDON W J, THIEL L C. Transfinite Mappings and Their Application to Grid Generation[J]. Applied Mathematics & Computation,1982,10/11(4):171-233.

[3] MUKHERJEE N, PEDDI B, CABELLO J, et al. Automatic Hexahedral Sweep Mesh Generation of Open Volumes[C]//International Meshing Roundtable. Orlando,2013:333-347.

[4] 陳建軍, 肖周芳, 曹建,等. 多源掃掠體全六面體網(wǎng)格自動(dòng)生成算法[J].浙江大學(xué)學(xué)報(bào)工學(xué)版,2012,46(2):274-279. CHEN Jianjun, XIAO Fangzhou, CAO Jian, et al. Automatic Hexahedral Mesh Generation for Many-to-one Sweep Volume[J]. Journal of Zhejiang University,2012,46(2):274-279.

[5] ZHANG H, ZHAO G. Adaptive Hexahedral Mesh Generation Based on Local Domain Curvature and Thickness Using a Modified Grid-based Method[J]. Finite Elements in Analysis & Design,2007,43(9):691-704.

[6] 黃麗麗,趙國(guó)群. 基于柵格法的三維六面體網(wǎng)格質(zhì)量?jī)?yōu)化[J]. 中國(guó)機(jī)械工程,2009,20(21):2603-2608. HUANG Lili, ZHAO Guoqun. Optimization of Grid-based Three-dimensional Hexahedral Meshes[J]. China Mechanical Engineering,2009,20(21):2603-2608.

[7] SUN L, ZHAO G. Adaptive Hexahedral Mesh Generation and Quality Optimization for Solid Models with Thin Features Using a Grid-based Method[J]. Engineering with Computers,2016,32(1):61-84.

[8] KREMER M, BOMMES D, LIM I, et al. Advanced Automatic Hexahedral Mesh Generation from Surface Quad Meshes[C]// Proceedings of International Meshing Roundtable.Orlando,2013:147-164.

[9] KAWAMURA Y, ISLAM M S, SUMI Y. A Strategy of Automatic Hexahedral Mesh Generation by Using an Improved Whisker-weaving Method with a Surface Mesh Modification Procedure[J]. Engineering with Computers,2008,24(3):215-229.

[10] TAM T K H, ARMSTRONG C G. Finite Element Meshcontrol by Integer Programming[J]. International Journal for Numerical Methods in Engineering,1993,36(15):2581-2605.

[11] PRICE M A, ARMSTRONG C G, SABIN M A. Hexahedral Mesh Generation by Medial Surface Subdivision: Part I. Solids with Convex Edges[J]. International Journal for Numerical Methods in Engineering,1995,38(19):3335-3359.

[12] 孫家廣. 計(jì)算機(jī)圖形學(xué)[M]. 3版.北京: 清華大學(xué)出版社, 1998. SUN Jiaguang. Computer Graphics[M]. 3rd ed. Beijing: Tsinghua University Press,1998.

[13] WHITELEY M, WHITE D, BENZLEY S, et al. Two and Three-quarter Dimensional Meshing Facilitators[J]. Engineering with Computers,1996,12(3):144-154.

[14] LI T S, ARMSTRONG C G, MCKEAG R M. Quad Mesh Generation fork-sided Faces and Hex Mesh Generation for Trivalent Polyhedral[J]. Finite Elements in Analysis & Design,1997,26(4):279-301.