超大型三維有限元模型的快速可視化算法

鄒 靜， 紀(jì)洪廣

（1. 清華大學(xué)水利水電工程系，北京 100084；2. 北京科技大學(xué)土木與環(huán)境學(xué)院，北京 100083）

超大型三維有限元模型的快速可視化算法

鄒 靜1,2， 紀(jì)洪廣2

（1. 清華大學(xué)水利水電工程系，北京 100084；2. 北京科技大學(xué)土木與環(huán)境學(xué)院，北京 100083）

針對(duì)超大型三維有限元模型的特點(diǎn)，總結(jié)了其可視化任務(wù)的基本需求和關(guān)鍵算法。由此分別為網(wǎng)格、單元、剖面可視化，以及繪制彩色云圖等不同任務(wù)，編制了結(jié)構(gòu)簡(jiǎn)單、運(yùn)行可靠的快速顯示算法。通過具體工程實(shí)例與典型商業(yè)程序的對(duì)比驗(yàn)證，初步證明了該算法的可行性和優(yōu)越性，有助于簡(jiǎn)化超大型有限元程序的編寫過程。

三維有限元; 可視化; 彩色云圖; 實(shí)體切割

隨著礦山開采實(shí)踐的不斷深入，三維數(shù)值分析以強(qiáng)大的分析能力逐漸成為主流。由于實(shí)際工程尺度的加大，三維計(jì)算模型的規(guī)?？焖僭鲩L(zhǎng)。近年來(lái)，單元數(shù)從十年前的數(shù)千增長(zhǎng)到數(shù)萬(wàn)，進(jìn)而增加到數(shù)十萬(wàn)甚至超越一百萬(wàn)。即使不考慮本構(gòu)模型復(fù)雜化的情況，礦山巖土分析也因此面臨很大壓力。如何快速、可靠地求解這些模型，是很有實(shí)際意義的問題。而在這些問題中，基于“所見即所得”的分析與測(cè)試需要，計(jì)算模型構(gòu)造與成果的快速可視化，是一個(gè)無(wú)法逾越的重要環(huán)節(jié)。盡管基本的算法和思路已不鮮見，不少軟件與程序包也已實(shí)現(xiàn)相應(yīng)功能。但要么程序過于龐大，要么源代碼受保護(hù)，編寫簡(jiǎn)單高效的可視化模塊始終存在不小挑戰(zhàn)。原有算法或程序雖然在原理上依然正確，但由于先前缺乏相應(yīng)的針對(duì)性，難免暴露問題出現(xiàn)反應(yīng)遲鈍的現(xiàn)象。在各種分析方法中，以有限元方法算法最成熟，同時(shí)其三維計(jì)算模型結(jié)構(gòu)簡(jiǎn)單，并且可與很多分析方法通用網(wǎng)格。因此，本文以有限元計(jì)算模型為研究對(duì)象，在實(shí)踐基礎(chǔ)上針對(duì)模型規(guī)模超大的情況，總結(jié)出關(guān)鍵、核心的算法需求，進(jìn)而構(gòu)造出簡(jiǎn)潔、高效的可視化實(shí)現(xiàn)流程。本項(xiàng)研究不僅有助于數(shù)值計(jì)算研究的進(jìn)展，對(duì)工程圖形領(lǐng)域的研究也有一定的借鑒意義。

1 核心算法

1.1 算法需求分析

有限元乃至其他算法三維模型的顯示，無(wú)外乎包括顯示要素與結(jié)果兩個(gè)方面。要素即幾何要素，包括：點(diǎn)、網(wǎng)格、單元等。結(jié)果即計(jì)算結(jié)果，顯示方式等：云圖、切片、等值線圖等幾種?？梢哉J(rèn)為可視化任務(wù)的本質(zhì)是：表示基本元素及其連接關(guān)系。在超大規(guī)模有限元模型的背景下，最緊要的莫過于：省略不可見元素，抽取關(guān)鍵要素，從而有效降低繪圖工作量。據(jù)此，可將算法需求分為以下兩大類：

1） 去除重復(fù)元素

在計(jì)算模型規(guī)模較小的情況下，可以根據(jù)元素的拓?fù)錁?gòu)造完整地匯出圖形。即四面體單元繪成6條邊、4個(gè)面，或者塊體單元顯示為12條邊、6個(gè)面，繪制元素?cái)?shù)量為單元數(shù)的數(shù)倍之多。在這些元素當(dāng)中位置重疊而不可見的占大部分，必須顯示的元素實(shí)際常不到總數(shù)的1/10。因此，當(dāng)計(jì)算模型規(guī)模十分龐大時(shí)，在元素序列中快速刪除重復(fù)元素，是加快繪圖速度一項(xiàng)十分重要的算法需求。

2） 明確元素間連接關(guān)系

邊的繪制取決于點(diǎn)的連接關(guān)系，面的繪制取決于邊的連接方式，而顯示實(shí)體則要涉及面的鄰接關(guān)系。各種關(guān)系當(dāng)中以邊的鄰接關(guān)系，也就是頂點(diǎn)的連接順序最為重要，切片圖和等值線圖都要依此為基礎(chǔ)。具體明確點(diǎn)素間連接關(guān)系方法，用數(shù)學(xué)語(yǔ)言描述實(shí)際如圖1所示，即為無(wú)向圖的深度優(yōu)先遍歷。

圖1 無(wú)向圖的深度優(yōu)先遍歷

1.2 關(guān)鍵算法選擇

根據(jù)上面總結(jié)的算法需求，結(jié)合C/C++語(yǔ)言與 Opengl功能的特點(diǎn)，幾個(gè)關(guān)鍵算法可以歸納如下。

1） 快速排序

使用排序可將相同元素移到相鄰位置，然后在只要從頭至尾比較前后位置，即可檢測(cè)到所有內(nèi)容重復(fù)的元素。如果操作大尺度元素序列計(jì)算量必然巨大，此時(shí)可采用快速排序算法[1]來(lái)處理。雖然采用查找再插入的常規(guī)方法可在空間復(fù)雜度上占優(yōu)勢(shì)，但其時(shí)間復(fù)雜度大致為(1/800～1/8)O(n2)，當(dāng)元素很多時(shí)遠(yuǎn)大于快速排序的O(n)。在各種現(xiàn)有的快速排序?qū)崿F(xiàn)當(dāng)中，以C語(yǔ)言的qsort庫(kù)函數(shù)速度最快且性能最為穩(wěn)定。該函數(shù)由stdlib.h文件引用，具體調(diào)用方式為：

_CRTIMP void __cdecl qsort(void *, size_t, size_t, int ( __cdecl *) (const void *, const void *));

參數(shù)1為待排序數(shù)據(jù)塊的首地址，參數(shù)2為數(shù)據(jù)塊單元個(gè)數(shù)，參數(shù)3為數(shù)據(jù)塊的單元尺寸，參數(shù)4為自定義的單元比較函數(shù)，輸入?yún)?shù)為單元首地址，輸出變量為表示單元大小關(guān)系的整數(shù)。

假設(shè)數(shù)列的元素長(zhǎng)度為n，且已經(jīng)為每個(gè)元素的組成要素排序，則元素E1、E2的比較函數(shù)算法可表述如下：

for i 1 to n

if E1[i]> E2[i]

return 1

else if E1[i]< E2[i]

return -1

end

end

return 0

其中，返回1表示E1> E2，返回-1表示E1< E2，而0表示E1= E2。

2） 深度優(yōu)先遍歷

由于鄰接表（見圖2）主要用于繪制切片，本文無(wú)向圖的節(jié)點(diǎn)數(shù)不會(huì)超越 12個(gè)，一個(gè)節(jié)點(diǎn)所連的邊數(shù)也不會(huì)超過3個(gè)。因此，選擇12×12的整數(shù)數(shù)組來(lái)表示，其具體的構(gòu)造方式如圖所示，深度優(yōu)先遍歷全圖的流程如圖3所示。

添加邊(i, j)的算法為：

l ←K(2, i)

K(l, i)←j

K(2, i)←K(2, i)+1

K(3, i)←K(3, i)+1

l ←K(2, j)

K(l, j)←i

K(2, j)←K(2, j)+1

K(3, j)←K(3, j)+1

刪除邊(i, j)的算法為：

l

←K(2,i)

K(l, i)←-1

K(3, i)←K(3, i)+1

l ←K(2,j)

K(l, j)←-1

K(3, j)←K(3, j)+1

圖2 鄰接表構(gòu)造

圖3 遍歷無(wú)向圖流程

3） 圖形顯示

本文圖形顯示過程如圖4所示，先生成矢量模型再調(diào)用 Opengl繪制圖形，以免單元、節(jié)點(diǎn)數(shù)據(jù)過多造成繪圖緩慢。其中用到的Opengl線、面繪制語(yǔ)句，以及若干繪圖模式設(shè)定語(yǔ)句，可以歸納匯總?cè)缦拢?/p>

圖4 繪圖過程

繪制線條：

glBegin (GL_LINES);

glColor3f(1.0f,1.0f,1.0f);

glVertex3f (x1,y1,z1);

glColor3f(1.0f,1.0f,1.0f);

glVertex3f (x2,y2,z2);

glEnd ();

繪制面：

glBegin (GL_POLYGON);

glColor3f(1.0f,1.0f,1.0f);

glVertex3f (x1,y1,z1);

…

glEnd();

繪制云圖時(shí)采用平滑模式：

glShadeModel(GL_SMOOTH);

為了顯示正確的三維效果，調(diào)用以下兩條語(yǔ)句啟用深度測(cè)試：

glEnable(GL_DEPTH_TEST);

glDepthFunc(GL_LESS);

由于線消隱與面消隱機(jī)制存在差異，同時(shí)繪制面與其邊線時(shí)常出現(xiàn)線被面掩蓋的現(xiàn)象。文獻(xiàn)[2]中采用面模擬線的方法回避這個(gè)問題，而筆者經(jīng)嘗試后發(fā)現(xiàn)：如果先繪線再繪面并啟用線平滑模式，也能較好地解決問題。調(diào)用線平滑模式的語(yǔ)句如下：

glEnable(GL_LINE_SMOOTH);

glHint(GL_LINE_SMOOTH_HINT,

GL_FASTEST);

2 典型可視化任務(wù)的實(shí)現(xiàn)算法

有限元和其它數(shù)值計(jì)算方法中，最典型的可視化任務(wù)有：節(jié)點(diǎn)、網(wǎng)格、單元與切片的可視化，以及數(shù)據(jù)云圖繪制等4種類型。相應(yīng)的實(shí)現(xiàn)算法簡(jiǎn)述如下，單元具體構(gòu)造參見有關(guān)經(jīng)典算法[3-4]。

2.1 網(wǎng)格可視化

網(wǎng)格可視化的關(guān)鍵如圖5所示。先根據(jù)單元構(gòu)造繪出所有邊線，再合并位置重疊（實(shí)為端點(diǎn)編號(hào)重合）的邊線，從而生成可以透視的骨架網(wǎng)格。令數(shù)組e[1…ne]代表邊線集合，則網(wǎng)格可視化偽代碼可表述如下。為了準(zhǔn)確比較e[i]和e[l]大小，已先為邊線節(jié)點(diǎn)排好順序。

l←0

for i 1 to ne

if e[i]≠e[l]

wend： 繪制邊線e[l]

l←i

if i=ne

goto wend

end

end

end

圖5 邊線合并

2.2 單元可視化

單元可視化重點(diǎn)在于外表面的可視化，實(shí)際上就是顯示單元群外殼的過程，所以單元可視化可歸結(jié)為圖6所示的外殼抽取[5-6]過程。具體操作方法為：先根據(jù)單元構(gòu)造生成面元列表f [1…nf ]（其元素為排序后的面節(jié)點(diǎn)編號(hào)），然后由快速排序?qū)⑾嗤嬉频较噜徫恢茫儆上旅鎮(zhèn)未a如下表述的算法，刪除重復(fù)面元抽取外殼繪成圖形。

圖6 外殼抽取

l←0

for i 1 to fe

if f [i]≠f [l]

wend：按照存儲(chǔ)的節(jié)點(diǎn)順序，重排f [l]的節(jié)點(diǎn)次序

繪制邊界面f [l]

l←i

if i=nf

goto wend

end

else

i←i+1

l←i

end

end

需要注意的是，使用 Opengl同時(shí)繪制大量線與面時(shí)，繪制線比繪制面會(huì)消耗更多時(shí)間，因而往往會(huì)明顯拖慢繪圖的整體速度。因此，還需要按照上一小節(jié)的方法，刪除多余邊線抽取表面網(wǎng)格單獨(dú)繪圖。

2.3 云圖繪制

云圖繪制核心算法與上一小節(jié)基本相同，區(qū)別在于后者只給單個(gè)面繪上單一顏色，而前者則需要根據(jù)數(shù)值繪成漸變顏色。超大型有限元模型的單元已經(jīng)足夠密集，不需要像文獻(xiàn)[7]一樣布置精確的顏色方案。只需為每個(gè)頂點(diǎn)設(shè)定單獨(dú)顏色值即可，數(shù)值選擇取則決于節(jié)點(diǎn)變量值的相對(duì)大小。從而，首先確定繪制的變量是否為節(jié)點(diǎn)變量，如果為單元變量則還需要轉(zhuǎn)化類型。然后，按照一定規(guī)律確定顏色值集合——也即色譜，根據(jù)變量相對(duì)值選定顏色數(shù)值即可繪制云圖，簡(jiǎn)要的操作流程如圖7所示。

圖7 云圖繪制流程

研究matlab工具箱的.m源文件，可以得到常見色譜的計(jì)算方法。其中，灰度色譜計(jì)算式如式(1)所示，中n表示色譜的長(zhǎng)度，而r、g、b分別表示三基色數(shù)值（取值范圍為[0,1]）。

hsv色譜的計(jì)算式如下式(2)與式(3)所示

2.4 剖面可視化

制作切片的方法大致有：文獻(xiàn)[8]所述的直接切割單元法，以及文獻(xiàn)[9]調(diào)用Opengl的圖像剖切法。后者算法較為復(fù)雜且難以輸出剖面數(shù)據(jù)，本文采取前一種方法的思路，但對(duì)各種類型單元使用同一切割方案，不再區(qū)分不同的相交情況。主要思路如圖8所示：查找與剖面相交的單元，然后逐一剖切再組合成完整的剖面。判斷一個(gè)單元是否與剖面相交，根據(jù)其在剖面上下側(cè)是否均有其節(jié)點(diǎn)。給定節(jié)點(diǎn)(x, y, z)與剖面A(x-x0)+B(y-y0) +C(z-z0)+D=0，相應(yīng)位置關(guān)系判計(jì)算方法見式(4)

其中，d>0表示點(diǎn)在剖面上側(cè)，d=0表示點(diǎn)在剖面中間，d<0表示點(diǎn)在剖面下側(cè)。

圖8 相交單元檢測(cè)

在作數(shù)據(jù)切片時(shí)需要插取交點(diǎn)數(shù)值，僅作單元切片時(shí)則不需要插值，具體操作步驟簡(jiǎn)述如下：

Step 1 計(jì)算節(jié)點(diǎn)與剖面的位置關(guān)系，據(jù)此抽取所有與剖面相交的單元，記錄單元編號(hào)集合成數(shù)組。

Step 2 釋放所有單元邊線形成交線數(shù)組，同時(shí)在單元數(shù)組中記錄邊的編號(hào)，以及其在單元內(nèi)部的編號(hào)。

Step 3 刪除重復(fù)元素緊縮邊線數(shù)組，并更新相交單元數(shù)組對(duì)應(yīng)編號(hào)。

Step 4 計(jì)算每條邊與剖面的交點(diǎn)，同時(shí)插取交點(diǎn)處數(shù)值，并記錄成交點(diǎn)數(shù)組。

Step 5 逐個(gè)單元?dú)w納交線，構(gòu)造并填充交點(diǎn)鄰接表。

Step 6 遍歷交點(diǎn)鄰接表形成單元剖面，然后組合成最終的模型剖面。

3 應(yīng)用實(shí)例

為了驗(yàn)證本文算法的實(shí)際性能，選取以下3個(gè)實(shí)例與現(xiàn)有同類軟件作對(duì)比試驗(yàn)。實(shí)現(xiàn)語(yǔ)言為C++且混合使用devcpp 5/vc++ 6.0編程，編譯環(huán)境是32位Windows Vista Business系統(tǒng)。機(jī)器配置為：Intel 酷睿 2.4GHz雙核處理器, 頻率1066MHz的2G容量DDR3內(nèi)存，頻率475MHz、1G容量集成顯卡。

1） 規(guī)則網(wǎng)格

本算例幾何尺寸為 1×1×1m，單元?jiǎng)澐殖叽?00×100×100，共計(jì)100萬(wàn)個(gè)八節(jié)點(diǎn)規(guī)則單元。最終運(yùn)算所得單元、網(wǎng)格圖形如圖9所示，為了保持圖形紙面清晰程度網(wǎng)格只展示平面圖形。當(dāng)采用大型商業(yè)有限元程序 ANSYS對(duì)比時(shí)，不但生成網(wǎng)格明顯吃力，而且啟動(dòng)單元圖形也十分遲鈍（15s以上），而運(yùn)行本文算法時(shí)則沒有此類問題。

2） 自由網(wǎng)格

本算例為起伏的地表模型，使用網(wǎng)格自由劃分工具tetgen（網(wǎng)址為：http://tetgen.berlios.de/）劃分網(wǎng)格，得到單元數(shù) 842245、節(jié)點(diǎn)數(shù) 149885的全四面體網(wǎng)格模型。根據(jù)本文算法得到的圖形如圖10所示，其中剖切面為過模型中心點(diǎn)且垂直于X軸的平面，繪圖視線方向平行于X軸正方向。本文算法顯示這些圖形基本流暢，當(dāng)將其導(dǎo)入巖土工程常用數(shù)值計(jì)算工具Flac3D，發(fā)現(xiàn)繪制單元或云圖大致需要3秒的反應(yīng)時(shí)間。

3） 復(fù)雜網(wǎng)格

圖11是以某礦的實(shí)際地下巷道為基礎(chǔ)，采用tetgen劃分出的復(fù)雜四面體網(wǎng)格，總計(jì)節(jié)點(diǎn)數(shù)164681而單元數(shù)為588129。其中，圖11(a)圖為整體單元圖形，圖11(b)圖為前者方框處放大結(jié)果。該模型的特點(diǎn)為節(jié)點(diǎn)數(shù)相對(duì)較多，所以單元邊線數(shù)量十分巨大。使用tetgen自帶可視化工具tetview查看時(shí)，移動(dòng)、轉(zhuǎn)動(dòng)模型有些遲鈍（有時(shí)停頓 10s以上）。而如果使用基于本文算法編制的工具，則生成、防縮、移動(dòng)以及轉(zhuǎn)動(dòng)模型都基本流暢。

4 結(jié) 論

從上面的對(duì)比實(shí)驗(yàn)可以看出：本文算法不僅能正確完成指定功能，而且在運(yùn)行速度、穩(wěn)定性上優(yōu)于一些現(xiàn)有商業(yè)軟件或免費(fèi)工具。如果再充分調(diào)用 Opengl的圖形加速功能，算法的運(yùn)行速度還能進(jìn)一步提升。因此，本文算法能更好滿單元數(shù)量不斷增長(zhǎng)的需求，并且算法結(jié)構(gòu)簡(jiǎn)單容易實(shí)現(xiàn)，可在面向超大型三維數(shù)值模型的有限元（或他方法）程序編制時(shí)使用。

圖9 簡(jiǎn)單網(wǎng)格

圖10 自由網(wǎng)格

圖11 復(fù)雜網(wǎng)格

[1] 朱經(jīng)緯. 三維數(shù)據(jù)場(chǎng)的三維重建與模型的虛擬切割研究[D]. 武漢: 華中科技大學(xué), 2007.

[2] 周國(guó)祥, 王春艷. 基于 OpenGL 三維非均勻 FDTD網(wǎng)格圖形的消隱處理[J]. 計(jì)算機(jī)應(yīng)用研究, 2008, 25(1): 285-287.

[3] Smith I M, Griffiths D V. Programming the finite element method [M]. Chichester, UK: John Wiley & Sons, 2004: 587.

[4] Liu G R, Quek S S. Finite element method: a practical course [M]. London UK: Butterworth-Heinemann Publications, 2003: 184.

[5] 陳良玉, 楊文通. 一種大型三維有限元網(wǎng)格的顯示和高效消隱方法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 1997, 9(2): 164-167.

[6] 余衛(wèi)平. 有限元網(wǎng)格圖形處理技術(shù)及計(jì)算結(jié)果的可視化[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2003, 15(12): 1561-1565.

[7] 李建波, 陳健云, 林 皋. 針對(duì)三維有限元數(shù)據(jù)場(chǎng)的精確后處理算法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2004, 16(8): 1169-1175.

[8] 梁振光, 唐任遠(yuǎn). 電磁場(chǎng)有限元結(jié)果的剖切顯示[J].電機(jī)與控制學(xué)報(bào), 2004, 8(3): 242-246.

[9] 楊曉松, 顧元憲, 李云鵬, 等．有限元網(wǎng)格體繪制中的剖切算法[J]. 中國(guó)圖象圖形學(xué)報(bào), 2002, 7(1): 55-62.

Fast visualization algorithm for huge 3D-finite element models

Zou Jing1,2, Ji Hongguang2
( 1. Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China; 2. School of Civil & Environment, Beijing University of Sciences & Technology, Beijing 100083, China )

According to the features of huge 3D-FEM models, basic requirements of visualization and key algorithms are summarized. Then, a simple and reliable algorithm of fast rendering is prepared for different tasks of grid, unit, slice visualization and color contour respectively. Through contrastive testification of specific examples and typical business software by contrast, the feasibility and superiority of the algorithm is initially demonstrated, which can contribute to simplification of preparation for huge FEM program.

3D-FEM; visualization; color contour; solid cutting

TP 391.72

A

2095-302X (2012)04-0013-07

2010-07-02

國(guó)家“863”重點(diǎn)課題資助項(xiàng)目（2008AA062104）；“十一五”支撐課題資助項(xiàng)目（2008BAB33B03）；國(guó)家“973”重點(diǎn)課題資助項(xiàng)目（2010CB226803，2010CB731501）

鄒 靜（1982-），男，湖南新化人，博士，主要研究方向?yàn)榱W(xué)數(shù)值計(jì)算程序的可視化算法設(shè)計(jì)。