邵鐵政，李世森

（天津大學(xué)，天津300072）

隨著計(jì)算機(jī)能力的不斷強(qiáng)大，三維空間內(nèi)的有限元計(jì)算被越來(lái)越廣泛的使用，因此三維空間內(nèi)的網(wǎng)格劃分也隨之不斷地被研究和改進(jìn)。三維問(wèn)題所采用的網(wǎng)格劃分單元一般是四面體、六面體。而四面體網(wǎng)格具有靈活性較高及能夠更好的逼近邊界的特點(diǎn)。四面體網(wǎng)格生成的算法已經(jīng)有許多重要的進(jìn)展，如局部變換法、Delaunay 算法、四叉樹(shù)/八叉樹(shù)算法、單元生長(zhǎng)法、網(wǎng)格前沿法等，其中以Delaunay 算法應(yīng)用最廣，因?yàn)槠渚哂袊?yán)謹(jǐn)?shù)臄?shù)學(xué)理論證明作為基礎(chǔ)，能夠由初始散亂點(diǎn)生成形態(tài)優(yōu)化的網(wǎng)格，故稱為Delaunay 網(wǎng)格［1］。

相對(duì)于二維領(lǐng)域而言，目前Delaunay 算法在三維領(lǐng)域內(nèi)的研究還不夠成熟，因?yàn)槿S空間內(nèi)，點(diǎn)、邊、面的管理及單元之間的關(guān)系更加復(fù)雜，單元重疊的可能性較大，而且不易檢驗(yàn)。經(jīng)研究大部分Delaunay 算法［1-3］，散亂點(diǎn)集的外邊界進(jìn)行Delaunay 三角剖分,而且需要先在散亂點(diǎn)內(nèi)部形成一個(gè)初始的四面體剖分，然后通過(guò)交換或者尋找?guī)缀侮P(guān)系將內(nèi)部四面體的剖分更新為新的Delaunay 四面體剖分，這樣使得網(wǎng)格劃分過(guò)程的前期工作量較大。本文適用于外包面為凸的散亂點(diǎn)四面體剖分，不需要內(nèi)部初始四面體劃分，提高了四面體網(wǎng)格劃分的效率。

1 三維散亂點(diǎn)Delaunay 四面體網(wǎng)格生成方法

1.1 相關(guān)概念

Delaunay 算法源于1934 年B.Delaunay 提出的Delaunay 準(zhǔn)則［4］。

三維Delaunay 準(zhǔn)則的定義如下：

（1）定義1:對(duì)于三維空間內(nèi)由一堆散點(diǎn)組成的四面體網(wǎng)格，若每一個(gè)四面體的外接球均不包含除此四面體頂點(diǎn)以外的網(wǎng)格點(diǎn)，則此四面體稱為Delaunay 四面體，由這些四面體形成的網(wǎng)格稱為Delaunay 網(wǎng)格。在已知大量初始散點(diǎn)的情況下，若不存在多點(diǎn)共球的情況，則由這些散點(diǎn)形成的Delaunay 網(wǎng)格是唯一的［5-6］。

為了方便尋找Delaunay 四面體，本文引入了一個(gè)新的角度概念——球缺角。由二維Delaunay 三角形圓心角判定方法推出三維空間下的部分相關(guān)定義如下。

為了闡述方便本文引入了平面域內(nèi)、域外定義：

（2）定義2:如圖1-a 以平面ABC 為例，右手展開(kāi)拇指與其余四指垂直，四指沿A-B-C 方向握拳，拇指所指方向的空間，即平面ABC 上部O 點(diǎn)所在空間為域內(nèi)，反方向的空間及平面ABC 則為域外。

四面體球缺角定義如下：

（3）定義3：如圖1-a，ΔABC 為在三維空間內(nèi)一個(gè)三角形，D 為在空間內(nèi)與之構(gòu)成四面體的點(diǎn)，O 為A、B、C、D 四點(diǎn)的外接球球心，半徑為R。OA 為O 與ΔABC 任意頂點(diǎn)形成向量，則OA與ΔABC 的負(fù)法向量n的夾角就是點(diǎn)D 對(duì)ΔABC 的球缺角δ(0＜δ＜π)。

（4）四面體球缺角的特性1:由以上定義易知D 對(duì)ΔABC 的球缺角δ 唯一。

（5）四面體球缺角的特性2:在空間內(nèi)當(dāng)所有散亂點(diǎn)位于ΔABC 域內(nèi)，一點(diǎn)與已知三角形構(gòu)成Delaunay 四面體的條件是：這個(gè)點(diǎn)對(duì)ΔABC 的球缺角最大。

1.2 相關(guān)證明

（1）特性1 證明:根據(jù)定義3 可知，四面體球缺角即為球心與三角形所構(gòu)成的圓錐的圓錐角的一半，所以得證。

（2）特性2 證明:圖1-b 為圖1-a 的正視圖，球O 過(guò)A、B、C、D 四點(diǎn)，r 為已知ΔABC 外接圓半徑，R 為球O 的半徑，h 為O 到ΔABC 的距離，H 為ΔABC 域內(nèi)球缺高，則球缺的體積

若D-ABC 為Delaunay 四面體，則根據(jù)[定義1]只需證明球缺內(nèi)不包含任何點(diǎn)。若證明[特性2]，只需證明球內(nèi)任意點(diǎn)D1對(duì)ΔABC 的球缺角大于D 對(duì)ΔABC 的球缺角即可。

∵h(yuǎn)=r×cotδ，R=r/sinδ，H=h+R

∴球缺體積V（δ）=（r/sinδ+r×cotδ）2×[3×（r/sinδ）-（r×cotδ+r/sinδ）] ×π/3

=（-cos3δ+3×cosδ+2）×πr3/（3×sin3δ）

由以上推導(dǎo)公式可知:

對(duì)于已經(jīng)固定的ΔABC 外接圓半徑r 不變，即球缺體積為以δ 為自變量的函數(shù)，則可求出并化簡(jiǎn)得到V-δ 函數(shù)的一階導(dǎo)數(shù)

V′（δ）=-（1+cosδ）2×πr3/sin4δ

∴由于(0＜δ＜π)，V′＜0

得證V 與δ 反相關(guān)

如圖1-b 易知在域內(nèi)，D1點(diǎn)對(duì)ΔABC 的球缺體積必然小于點(diǎn)D 對(duì)ΔABC 的球缺體積。即D1點(diǎn)對(duì)ΔABC球冠角δ 更大。

當(dāng)域內(nèi)球缺小于半球時(shí)，δ 的范圍是（π/2，π），即cotδ 為負(fù)，證明過(guò)程相同。

所以［特性2］得證。

2 編程實(shí)現(xiàn)散亂點(diǎn)四面體劃分的步驟

2.1 散亂點(diǎn)數(shù)據(jù)讀入

數(shù)據(jù)文件要求：散亂點(diǎn)按統(tǒng)一格式存入inputA.txt 文件中，順序?yàn)椤俺跏键c(diǎn)號(hào)，x 坐標(biāo)，y 坐標(biāo)，z 坐標(biāo)”。

2.2 程序運(yùn)行思維圖

如圖2 所示，程序運(yùn)行分為以下幾個(gè)步驟：

（1）三維空間點(diǎn)集獲取：讀取輸入文件中的所有點(diǎn)的坐標(biāo)，讀取內(nèi)容包括“初始點(diǎn)號(hào)，x 坐標(biāo)，y 坐標(biāo)，z 坐標(biāo)”。

（2）初始化三角形：讀取文件中的第一、二、三點(diǎn)為初始三角形。

（3）點(diǎn)定位搜索：由初始三角形根據(jù)四面體球缺角的特性2 在散亂點(diǎn)中尋找與之構(gòu)成Delaunay 四面體的點(diǎn)，并生成Delaunay 四面體。

（4）循環(huán)運(yùn)行搜索：在生成每一個(gè)四面體的同時(shí)，有3 個(gè)新的三角形生成，分別以新生成的三角形為初始三角形繼續(xù)尋找對(duì)應(yīng)的點(diǎn)。由此循環(huán)運(yùn)行直至將所有的散亂點(diǎn)進(jìn)行四面體劃分。

2.3 搜索方法

步驟（3）中搜索方法分為以下幾個(gè)部分：

（1）確定搜索區(qū)域:為減少程序的運(yùn)算量，將所有散亂點(diǎn)分為域內(nèi)、域外2 個(gè)部分（方法見(jiàn)3-(3)）,只對(duì)域內(nèi)部分進(jìn)行搜索。

（2）計(jì)算球心坐標(biāo)：分別將域內(nèi)的所有點(diǎn)按順序分別與初始三角形的3 個(gè)頂點(diǎn)組合，由子程序計(jì)算此空間四點(diǎn)的外接球的球心坐標(biāo)O。

（3）計(jì)算保存向量：根據(jù)三角形頂點(diǎn)坐標(biāo)求出三角形外接圓圓心坐標(biāo)O1，向量OO1即為初始三角形負(fù)法向量。保存球心O 與任意三角形頂點(diǎn)形成的向量OA (根據(jù)[特性一]，O 與任意頂點(diǎn)形成的向量與OO1夾角大小相同)。

（4）計(jì)算OA 與OO1夾角

計(jì)算出OA 與OO1夾角，并保存，以便經(jīng)歷一次搜索后找出最大的夾角，并確定為最大球缺角，該角對(duì)應(yīng)的點(diǎn)與此初始三角形構(gòu)成Delaunay 四面體。

3 程序的準(zhǔn)確性和高效性

為提高程序運(yùn)行的準(zhǔn)確性和效率，本程序進(jìn)行了以下幾方面的處理:

（1）順序記錄：在每次計(jì)算所得的新三角形點(diǎn)數(shù)據(jù)，都按照順時(shí)針?lè)较蜻M(jìn)行記錄,保存的新生成的數(shù)據(jù)格式符合初始運(yùn)算格式，使程序可以循環(huán)進(jìn)行計(jì)算。

（2）查重：由于初始三角形不同，對(duì)于每次生成的新的四面體來(lái)說(shuō)有可能與前面生成的四面體重復(fù)，因此，對(duì)新生成的四面體進(jìn)行查重，假設(shè)生成的新四面體4 個(gè)點(diǎn)的序號(hào)分別為（a，b，c，d），將(a，b，c，d)分別與之前生成的各數(shù)組對(duì)比，如果出現(xiàn)與（a，b，c，d）數(shù)值構(gòu)成相同的數(shù)組，則取消該四面體，不存入，由此循環(huán)下去，并不斷檢驗(yàn)是否重復(fù)，直到令所有點(diǎn)都找到Delaunay 四面體則計(jì)算過(guò)程結(jié)束。

（3）域內(nèi)域外判斷：對(duì)于新生成的三角形來(lái)說(shuō)一部分點(diǎn)位于該三角形的域內(nèi)，另一部分位于域外。本程序?qū)λ悬c(diǎn)與初始三角形的位置關(guān)系進(jìn)行了預(yù)判，域外的三角形不必要進(jìn)行復(fù)雜的計(jì)算，因此提高了程序的運(yùn)行。四面體體積判定法:根據(jù)三維行列式正負(fù)的意義，A（x1，y1，z1）、B(x2，y2，z2)、C（x3，y3，z3)為空間的三點(diǎn)，D（x0，y0，z0）與三點(diǎn)的行列式表示的是四面體的體積（當(dāng)D 位于平面ABC 的域內(nèi)時(shí)V 為正，否則V 為負(fù)）。

四面體體積

（4）避免計(jì)算誤差：本文采取了一系列方法避免計(jì)算過(guò)程中產(chǎn)生的誤差，其中以下2 個(gè)方法用于整個(gè)程序中相關(guān)的計(jì)算步驟。①為避免計(jì)算機(jī)在計(jì)算時(shí)出現(xiàn)除以0 的數(shù)學(xué)錯(cuò)誤，本程序在編譯時(shí)避免了除法運(yùn)算，全部使用正負(fù)號(hào)的判斷，避免了計(jì)算機(jī)儲(chǔ)存或計(jì)算中出現(xiàn)誤差的情況。②在步驟（2）、（3）中計(jì)算球心、外接圓圓心時(shí)，本程序采用克拉默法則對(duì)方程組進(jìn)行求解，但是在程序編譯過(guò)程中發(fā)現(xiàn)，編程平臺(tái)自帶的IMSL 函數(shù)庫(kù)在計(jì)算行列式過(guò)程中DET（）函數(shù)使用lin_sol_lsq（）來(lái)解方程［7］，對(duì)于行列式值為0 或非常接近于0 的矩陣會(huì)失效，所以本程序中使用了重新編譯的行列式計(jì)算函數(shù)進(jìn)行相關(guān)計(jì)算，消除了此類計(jì)算誤差。

4 時(shí)間復(fù)雜度分析

為了對(duì)本程序的計(jì)算效率進(jìn)行比較，本文進(jìn)行了程序運(yùn)行時(shí)間復(fù)雜度的分析［8］。圖3 為不同點(diǎn)數(shù)的情況下程序運(yùn)行時(shí)間的趨勢(shì)線情況。

根據(jù)趨勢(shì)線可以看出計(jì)算的點(diǎn)的數(shù)量y 與程序運(yùn)行時(shí)間x 的關(guān)系擬合曲線為y=8E-06x2，即時(shí)間與點(diǎn)數(shù)總量呈二次關(guān)系，該計(jì)算方法有良好的應(yīng)用特性。

5 結(jié)束語(yǔ)

球缺角定義的提出，使Delaunay 四面體的剖分有了更加量化的方法，區(qū)別于以往的定性判斷，量化的方法誤差更小，更容易控制和比較。雖然在本方法的實(shí)踐過(guò)程中仍有一些方面不完善，但是相信在不斷地改進(jìn)后，將對(duì)空間散亂點(diǎn)集Delaunay 四面體剖分方法的研究做出貢獻(xiàn)，可用于各個(gè)領(lǐng)域的工程應(yīng)用軟件的開(kāi)發(fā)。

［1］關(guān)文革，武強(qiáng)，賈麗萍，等.約束數(shù)據(jù)域Delaunay 四面體網(wǎng)格生成算法［J］.華中科技大學(xué)學(xué)報(bào)：自然科學(xué)版，2005（5）:67-69.GUAN W G，WU Q，JIA L P，et al. Algorithm of mesh generation of Delaunay tetrahedral in constrained domain［J］. Journal of Huazhong University of Science and Technology，2005（5）: 67-69.

［2］WU Q，XU H.An approach to computer modeling and visualization of geological faults in 3d［J］.International Journal of Geographical Information Systems，1993，7（6）:501-524.

［3］XU H，WU Q.Design & implementation of visualization for 3d samdwich geological bodies［J］. Computer Applictations，2001，21（12）:56-60.

［4］Lawson C L. Software for C1 surface interpolation，Mathematical Software III［M］. New York: Academic Press，1977: 161-194.

［5］陳學(xué)工，潘懋.空間散點(diǎn)集Delaunay 四面體剖分切割算法［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)報(bào)，2002（1）:93-95.CHEN X G，PAN M. Delaunay Triangulation Cutting Algorithm for A Set of Irregularly Located Spatial Points［J］. Journal of Computer Aided Design & Computer Graphics，2002（1）: 93-95.

［6］王建華，徐強(qiáng)尋，張銳.任意形狀三維物體的Delaunay 網(wǎng)格生成算法［J］.巖石力學(xué)與工程學(xué)報(bào)，2003，22（5）:717-720.WANG J H，XU Q X，ZHANG R. Delaunay algorithm and related procedure to generate the tetrahedron mesh for an object with arbitrary boundary［J］.Chinese Journal of Rock Mechanics and Engineering，2003，22（5）:717-720.

［7］李麗.三維空間Delaunay 三角剖分算法的研究及應(yīng)用［D］.大連：大連海事大學(xué)，2010.［8］彭國(guó)倫.Fortran95 程序設(shè)計(jì)［M］.北京：中國(guó)電力出版社，2010.