丁道紅,章 青

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室,江蘇南京 210098;2.河海大學(xué)力學(xué)與材料學(xué)院,江蘇南京 210098)

Voronoi單胞是以自然相鄰插值法進(jìn)行形函數(shù)構(gòu)造的自然單元法的基礎(chǔ)[1-2]。自然單元法形函數(shù)的構(gòu)造方法有Sibson法[3]和Laplace法[4]2種,其中Sibson法是以插值點對應(yīng)于某一結(jié)點的二階Voronoi單胞的勒貝格測度與插值點的一階Voronoi單胞的勒貝格測度的比值作為形函數(shù)。三維Voronoi單胞的勒貝格測度即為Voronoi單胞的體積,目前關(guān)于Voronoi單胞體積的計算最常用的方法是Lasserre方法[5-8]。Lasserre方法是在已知凸多面體控制方程的前提下,通過遞推公式求得凸多面體的體積。該方法因具有簡單可行、操作方便等優(yōu)點而得到廣泛的應(yīng)用,但該方法在應(yīng)用過程中要求控制方程不能存在冗余的方程。本文在與Voronoi單胞相對應(yīng)的Delaunay三角網(wǎng)的啟示下,通過Delaunay三角網(wǎng)將三維Voronoi單胞劃分成若干個不重疊且充滿Voronoi單胞的四面體,從而可簡單地通過四面體的體積計算得到Voronoi單胞的體積。

1 Voronoi單胞

對于Voronoi單胞的概念[9],可通過數(shù)學(xué)方式進(jìn)行描述。記Rn空間上任意2點xI和xJ的歐氏距離為d(xI,xJ),設(shè)P={x1,…,xn}為Rn空間上任意n個互異的點,則與其中任意一點xI對應(yīng)的Voronoi單胞的定義為

二階Voronoi單胞是在Rn空間上任意n個互異的點P={x1,…,xn}的基礎(chǔ)上,增加一新的待插值點x′,形成新的一階Voronoi單胞,其中結(jié)點xI為插值點x′的自然鄰點,結(jié)點xI的二階Voronoi單胞為結(jié)點xI的一階Voronoi單胞與插值點x′的一階Voronoi單胞重疊的部分,即

將具有公共邊界的Voronoi單胞所對應(yīng)的結(jié)點稱為自然鄰點,其公共邊界為該2點的中垂面(二維問題為2點的垂直平分線),即2自然鄰點xI和xJ的Voronoi單胞公共邊界所在的平面(或直線)L為

其法線方向為結(jié)點xI指向xJ。與結(jié)點xI對應(yīng)的Voronoi單胞在平面(或直線)L的結(jié)點xI所在的一側(cè),即結(jié)點xI對應(yīng)的Voronoi單胞滿足

假設(shè)結(jié)點xI(xI,yI,zI)為插值點xp(xp,yp,zp)的1個自然鄰點,與結(jié)點xI(xI,yI,zI)對應(yīng)的自然鄰點共有m個,分別為(x1,y1,z1),…,(xm,ym,zm),則與結(jié)點xI對應(yīng)的二階Voronoi單胞可寫成

令

則式(5)可寫為

式(6)描述了二階Voronoi單胞的控制域,若能得到二階Voronoi單胞的頂點,則可通過Delaunay三角網(wǎng)將控制域劃分成若干個小域進(jìn)行分析。因此若要得到二階Voronoi單胞的勒貝格測度,首先需要得到二階Voronoi單胞的頂點。二階Voronoi單胞的頂點是邊界的交點,即式(6)取等號時所得的滿足式(6)的點。對于n維空間問題,通過聯(lián)立式(6)中n個方程,建立1個新的方程組

若|A|≠0,則式(7)必有解

若x=A-1b滿足式(6),則為二階Voronoi單胞的1個頂點。

2 Delaunay三角網(wǎng)

Delaunay三角網(wǎng)的概念是建立在Voronoi單胞基礎(chǔ)之上的。將具有公共邊界的Voronoi單胞所對應(yīng)的結(jié)點連接所得到的三角形,稱作Delaunay三角形,將Rn空間上任意n個互異的點通過Delaunay三角形劃分,可形成Delaunay三角網(wǎng)。

文獻(xiàn)[10]證明了Delaunay三角網(wǎng)具有如下性質(zhì)。

a.唯一性：不論從區(qū)域何處開始構(gòu)建,最終都得到一致的結(jié)果。

b.最優(yōu)性：任意2個相鄰三角形形成的凸四邊形的對角線如果可以互換,那么2個三角形6個內(nèi)角中最小的角度不會變大。

c.最規(guī)則：如果將三角網(wǎng)中每個三角形的最小角進(jìn)行升序排列,則Delaunay三角網(wǎng)的排列得到的數(shù)值最大。

d.具有凸多邊形的外殼：三角網(wǎng)最外層的邊界形成1個凸多邊形的外殼。

以上這些性質(zhì)保證了Delaunay三角網(wǎng)所劃分的區(qū)域不具有重疊部分,且能覆蓋結(jié)點所占的凸區(qū)域。

目前關(guān)于Delaunay三角網(wǎng)劃分的方法很多[11-12],如Lawson算法、Bowyer-Watson算法等。此外,一些商業(yè)軟件具有成熟的Delaunay三角網(wǎng)劃分的功能,如MATLAB軟件中,可通過Delaunay命令進(jìn)行相關(guān)工作。

在三維問題中,通過Delaunay三角網(wǎng)劃分所得的基本體為四面體,四面體的體積可通過簡單的幾何方法得到,如已知某四面體的4個頂點分別為(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4),則該四面體的體積為

3 算 例

已知三維空間8結(jié)點1(0,0,0),2(1,0,0),3(1,0,1),4(0,0,1),5(0,1,0),6(1,1,0),7(1,1,1),8(0,1,1),插值點xp(0.7236,0.1382,0.1382),則結(jié)點6為插值點的自然鄰點。結(jié)點6的自然鄰點有3個,分別為2,5,7。結(jié)點6的二階Voronoi單胞為

其中

通過式(10)中任取3個方程并取等號進(jìn)行聯(lián)立方程,代入上式進(jìn)行判斷是否符合要求,可得結(jié)點6的二階Voronoi單胞的頂點為(0.5,0.5,-1.0854),(1.2927,0.5,0.5),(0.5,0.7542,0.5),(0.5,0.5,0.5)共4個結(jié)點,通過Delaunay三角網(wǎng),可劃分成1個四面體,通過式(8)可知其體積為0.0533。

4 結(jié) 語

在已知Voronoi單胞頂點的前提下,利用Delaunay三角網(wǎng)將Voronoi單胞劃分成若干四面體,通過求解四面體的體積得到Voronoi單胞的體積。該方法具有原理簡單、可行性強(qiáng)、理解性好等優(yōu)點。

此外,該方法不僅可以應(yīng)用于三維Voronoi單胞的體積計算,對于任意1個已知邊界頂點的凸多面體均可應(yīng)用該方法進(jìn)行體積的計算。

[1]BRAUN J,SAMBRIDGEM.A numerical method for solving partial differential equations on highly irregular evolving grids[J].Nature,1995,376：655-660.

[2]SUKUMAR N,MORAN B,BELYTSCHKO T.The natural element method in solid mechanics[J].International Journal for Numerical Methods in Engineering,1998,43(5)：839-887.

[3]SIBSON R.A brief description of natural neighbour interpolation[G]//BARNETT V.Interpreting Multivariate Data.New York：Wiley,1981：21-36.

[4]BELIKOV V V,SEMENOV A Y.Non-Sibsonian interpolation on arbitrary system of points in Euclidean space and adaptive isolines generation[J].Applied Numerical Mathematics,2000,32(4)：371-387.

[5]LASSERRE JB.An analytical expression and an algorithm for the volume of a convex polyhedron inRn[J].Journal of Optimization Theory and Applications,1983,39(3)：363-377.

[6]LASSERREJB.A laplace transform algorithm for the volume of a convex polytope[J].Journal of the ACM,2001,48(6)：1126-1140.

[7]王建華,張英新,高紹武.三維彈塑性自然單元法算法實現(xiàn)[J].計算力學(xué)學(xué)報,2006,23(5)：594-598.(WANG Jianhua,ZHANG Yingxin,GAO Shaowu.The computational methods of natural element method in three dimensional elasto-plastic analysis[J].Chinese Journal of Computational Mechanics,2006,23(5)：594-598.(in Chinese))

[8]江濤,章青.基于Lasserre算法的自然單元法形函數(shù)計算[J].力學(xué)與實踐,2008,30(4)：1-5.(JIANG Tao,ZHANGQing.The shape functions in the natural element method based on Lasserre's algorithm[J].Mechanics in Engineering,2008,30(4)：1-5.(in Chinese))

[9]楊永清,馮鈞,王志堅.基于Voronoi圖的復(fù)雜對象空間方位關(guān)系的推理計算[J].河海大學(xué)學(xué)報：自然科學(xué)版,2008,36(3)：414-417.(YANG Yongqing,FENG Jun,WANG Zhijian.Reasoning and calculation of spatial direction relationship between complex objects based on Voronoi diagram[J].Journal of Hohai University：Natural Sciences,2008,36(3)：414-417.(in Chinese))

[10]SEIDEL R.The upper bound theorem for polytopes：an easy proof of its asymptotic version[J].Computational Geometry,1995,5(2)：115-116.

[11]王建華,徐強(qiáng)勛,張銳.任意形狀三維物體的Delaunay網(wǎng)格生成算法[J].巖石力學(xué)與工程學(xué)報,2003,22(5)：717-722.(WANG Jianhua,XUQiangxun,ZHANG Rui.Delaunay algorithm and related procedure to generate the tetrahedron mesh for an object with arbitrary boundary[J].Chinese Journal of RocKMechanics and Engineering,2003,22(5)：717-722.(in Chinese))

[12]WATSON D F.Computing then-dimensional delaunay tessellation with application to Voronoi polytopes[J].The Computer Journal,1981,24(2)：167-172.