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

鄒 靜， 紀洪廣

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

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

鄒 靜1,2， 紀洪廣2

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

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

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

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

1 核心算法

1.1 算法需求分析

有限元乃至其他算法三維模型的顯示，無外乎包括顯示要素與結果兩個方面。要素即幾何要素，包括：點、網格、單元等。結果即計算結果，顯示方式等：云圖、切片、等值線圖等幾種。可以認為可視化任務的本質是：表示基本元素及其連接關系。在超大規模有限元模型的背景下，最緊要的莫過于：省略不可見元素，抽取關鍵要素，從而有效降低繪圖工作量。據此，可將算法需求分為以下兩大類：

1） 去除重復元素

在計算模型規模較小的情況下，可以根據元素的拓撲構造完整地匯出圖形。即四面體單元繪成6條邊、4個面，或者塊體單元顯示為12條邊、6個面，繪制元素數量為單元數的數倍之多。在這些元素當中位置重疊而不可見的占大部分，必須顯示的元素實際常不到總數的1/10。因此，當計算模型規模十分龐大時，在元素序列中快速刪除重復元素，是加快繪圖速度一項十分重要的算法需求。

2） 明確元素間連接關系

邊的繪制取決于點的連接關系，面的繪制取決于邊的連接方式，而顯示實體則要涉及面的鄰接關系。各種關系當中以邊的鄰接關系，也就是頂點的連接順序最為重要，切片圖和等值線圖都要依此為基礎。具體明確點素間連接關系方法，用數學語言描述實際如圖1所示，即為無向圖的深度優先遍歷。

圖1 無向圖的深度優先遍歷

1.2 關鍵算法選擇

根據上面總結的算法需求，結合C/C++語言與 Opengl功能的特點，幾個關鍵算法可以歸納如下。

1） 快速排序

使用排序可將相同元素移到相鄰位置，然后在只要從頭至尾比較前后位置，即可檢測到所有內容重復的元素。如果操作大尺度元素序列計算量必然巨大，此時可采用快速排序算法[1]來處理。雖然采用查找再插入的常規方法可在空間復雜度上占優勢，但其時間復雜度大致為(1/800～1/8)O(n2)，當元素很多時遠大于快速排序的O(n)。在各種現有的快速排序實現當中，以C語言的qsort庫函數速度最快且性能最為穩定。該函數由stdlib.h文件引用，具體調用方式為：

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

參數1為待排序數據塊的首地址，參數2為數據塊單元個數，參數3為數據塊的單元尺寸，參數4為自定義的單元比較函數，輸入參數為單元首地址，輸出變量為表示單元大小關系的整數。

假設數列的元素長度為n，且已經為每個元素的組成要素排序，則元素E1、E2的比較函數算法可表述如下：

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） 深度優先遍歷

由于鄰接表（見圖2）主要用于繪制切片，本文無向圖的節點數不會超越 12個，一個節點所連的邊數也不會超過3個。因此，選擇12×12的整數數組來表示，其具體的構造方式如圖所示，深度優先遍歷全圖的流程如圖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 鄰接表構造

圖3 遍歷無向圖流程

3） 圖形顯示

本文圖形顯示過程如圖4所示，先生成矢量模型再調用 Opengl繪制圖形，以免單元、節點數據過多造成繪圖緩慢。其中用到的Opengl線、面繪制語句，以及若干繪圖模式設定語句，可以歸納匯總如下：

圖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();

繪制云圖時采用平滑模式：

glShadeModel(GL_SMOOTH);

為了顯示正確的三維效果，調用以下兩條語句啟用深度測試：

glEnable(GL_DEPTH_TEST);

glDepthFunc(GL_LESS);

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

glEnable(GL_LINE_SMOOTH);

glHint(GL_LINE_SMOOTH_HINT,

GL_FASTEST);

2 典型可視化任務的實現算法

有限元和其它數值計算方法中，最典型的可視化任務有：節點、網格、單元與切片的可視化，以及數據云圖繪制等4種類型。相應的實現算法簡述如下，單元具體構造參見有關經典算法[3-4]。

2.1 網格可視化

網格可視化的關鍵如圖5所示。先根據單元構造繪出所有邊線，再合并位置重疊（實為端點編號重合）的邊線，從而生成可以透視的骨架網格。令數組e[1…ne]代表邊線集合，則網格可視化偽代碼可表述如下。為了準確比較e[i]和e[l]大小，已先為邊線節點排好順序。

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 單元可視化

單元可視化重點在于外表面的可視化，實際上就是顯示單元群外殼的過程，所以單元可視化可歸結為圖6所示的外殼抽取[5-6]過程。具體操作方法為：先根據單元構造生成面元列表f [1…nf ]（其元素為排序后的面節點編號），然后由快速排序將相同面移到相鄰位置，再由下面偽代碼如下表述的算法，刪除重復面元抽取外殼繪成圖形。

圖6 外殼抽取

l←0

for i 1 to fe

if f [i]≠f [l]

wend：按照存儲的節點順序，重排f [l]的節點次序

繪制邊界面f [l]

l←i

if i=nf

goto wend

end

else

i←i+1

l←i

end

end

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

2.3 云圖繪制

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

圖7 云圖繪制流程

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

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

2.4 剖面可視化

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

其中，d>0表示點在剖面上側，d=0表示點在剖面中間，d<0表示點在剖面下側。

圖8 相交單元檢測

在作數據切片時需要插取交點數值，僅作單元切片時則不需要插值，具體操作步驟簡述如下：

Step 1 計算節點與剖面的位置關系，據此抽取所有與剖面相交的單元，記錄單元編號集合成數組。

Step 2 釋放所有單元邊線形成交線數組，同時在單元數組中記錄邊的編號，以及其在單元內部的編號。

Step 3 刪除重復元素緊縮邊線數組，并更新相交單元數組對應編號。

Step 4 計算每條邊與剖面的交點，同時插取交點處數值，并記錄成交點數組。

Step 5 逐個單元歸納交線，構造并填充交點鄰接表。

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

3 應用實例

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

1） 規則網格

本算例幾何尺寸為 1×1×1m，單元劃分尺寸100×100×100，共計100萬個八節點規則單元。最終運算所得單元、網格圖形如圖9所示，為了保持圖形紙面清晰程度網格只展示平面圖形。當采用大型商業有限元程序 ANSYS對比時，不但生成網格明顯吃力，而且啟動單元圖形也十分遲鈍（15s以上），而運行本文算法時則沒有此類問題。

2） 自由網格

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

3） 復雜網格

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

4 結 論

從上面的對比實驗可以看出：本文算法不僅能正確完成指定功能，而且在運行速度、穩定性上優于一些現有商業軟件或免費工具。如果再充分調用 Opengl的圖形加速功能，算法的運行速度還能進一步提升。因此，本文算法能更好滿單元數量不斷增長的需求，并且算法結構簡單容易實現，可在面向超大型三維數值模型的有限元（或他方法）程序編制時使用。

圖9 簡單網格

圖10 自由網格

圖11 復雜網格

[1] 朱經緯. 三維數據場的三維重建與模型的虛擬切割研究[D]. 武漢: 華中科技大學, 2007.

[2] 周國祥, 王春艷. 基于 OpenGL 三維非均勻 FDTD網格圖形的消隱處理[J]. 計算機應用研究, 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] 陳良玉, 楊文通. 一種大型三維有限元網格的顯示和高效消隱方法[J]. 計算機輔助設計與圖形學學報, 1997, 9(2): 164-167.

[6] 余衛平. 有限元網格圖形處理技術及計算結果的可視化[J]. 計算機輔助設計與圖形學學報, 2003, 15(12): 1561-1565.

[7] 李建波, 陳健云, 林 皋. 針對三維有限元數據場的精確后處理算法[J]. 計算機輔助設計與圖形學學報, 2004, 16(8): 1169-1175.

[8] 梁振光, 唐任遠. 電磁場有限元結果的剖切顯示[J].電機與控制學報, 2004, 8(3): 242-246.

[9] 楊曉松, 顧元憲, 李云鵬, 等．有限元網格體繪制中的剖切算法[J]. 中國圖象圖形學報, 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

國家“863”重點課題資助項目（2008AA062104）；“十一五”支撐課題資助項目（2008BAB33B03）；國家“973”重點課題資助項目（2010CB226803，2010CB731501）

鄒 靜（1982-），男，湖南新化人，博士，主要研究方向為力學數值計算程序的可視化算法設計。