摘 要：本文介紹了現(xiàn)有的四面體網(wǎng)格質(zhì)量度量函數(shù)及網(wǎng)格優(yōu)化方法，本文采用錯(cuò)誤函數(shù)作為基于優(yōu)化光順?biāo)惴ǖ哪繕?biāo)函數(shù)，對(duì)四面體網(wǎng)格進(jìn)行光順，算法運(yùn)用多線(xiàn)程加速方式實(shí)現(xiàn)，并提出了一種預(yù)防死鎖的方法。實(shí)驗(yàn)結(jié)果表明：采用多線(xiàn)程加速方式，在頂點(diǎn)數(shù)量較大時(shí)，具有明顯優(yōu)越性。

關(guān)鍵詞：網(wǎng)格質(zhì)量度量準(zhǔn)則；網(wǎng)格優(yōu)化；多線(xiàn)程加速

中圖分類(lèi)號(hào)：TH123

虛擬手術(shù)仿真技術(shù)是近年來(lái)的研究熱點(diǎn)，作為組成三維仿真模型的基本單元——四面體網(wǎng)格的質(zhì)量好壞，則嚴(yán)重影響仿真的效果。由于所要分析的問(wèn)題越來(lái)越復(fù)雜，工程人員逐漸使用網(wǎng)格自動(dòng)生成工具來(lái)剖分網(wǎng)格，但得到的初始網(wǎng)格（三維情況）往往質(zhì)量很差，而網(wǎng)格單元質(zhì)量直接影響有限元計(jì)算中的收斂性，因此必須對(duì)得到的網(wǎng)格進(jìn)行優(yōu)化。本文將在總結(jié)已有四面體網(wǎng)格質(zhì)量的判定以及網(wǎng)格優(yōu)化的方法之后，運(yùn)用現(xiàn)有方法，使用多線(xiàn)程加速方式實(shí)現(xiàn)計(jì)算效率的提高。

1 單元質(zhì)量衡量準(zhǔn)則及網(wǎng)格優(yōu)化方法

長(zhǎng)期以來(lái)，對(duì)四面體單元質(zhì)量的度量或評(píng)判并沒(méi)有公認(rèn)的標(biāo)準(zhǔn)，但是對(duì)一個(gè)合理的單元質(zhì)量的度量準(zhǔn)則需要滿(mǎn)足的原則卻有著共識(shí)：?jiǎn)卧钠揭?、旋轉(zhuǎn)、反射、均勻縮放均應(yīng)不改變其度量值；當(dāng)且僅當(dāng)四面體為正四面體時(shí)其度量值取最大值1.0；當(dāng)四面體體積趨于零時(shí)，其度量值也接近于零[5]。以下為常用的度量準(zhǔn)則表達(dá)式如表1。

表1 幾種不同四面體網(wǎng)格質(zhì)量度量準(zhǔn)則

準(zhǔn)則表達(dá)式備注

ρ[1]

V是四面體的體積，Si為四面體各個(gè)面的面積，li為四面體的邊長(zhǎng)。

θ[2，3]θ=min（θ1，θ2，θ3，θ4） 可求得θ1，lij為連接頂點(diǎn)i和j的邊的長(zhǎng)度，θ2，θ3，θ4可通過(guò)指標(biāo)輪換得到。

Q[4]

Cd=1832.8208是為了使四面體的質(zhì)量度量值取最大值而采用的比例常數(shù)。

γ[2，5]

V是四面體的體積，lij為連接頂點(diǎn)i和j的邊的長(zhǎng)度

根據(jù)文獻(xiàn)11，本文最終選擇度量準(zhǔn)則ρ作為四面體單元質(zhì)量度量標(biāo)準(zhǔn)。

現(xiàn)有的網(wǎng)格優(yōu)化算法主要分為插入或刪除點(diǎn)[5]、邊交換或面交換[6]以及光順。傳統(tǒng)的Laplacian光順[7]做法是局部調(diào)整每一個(gè)網(wǎng)格點(diǎn)到與其相關(guān)聯(lián)的節(jié)點(diǎn)的中心，但不能保證一定能夠改善網(wǎng)格的質(zhì)量，且容易產(chǎn)生不合法單元。基于優(yōu)化的光順[8，9]是采用對(duì)目標(biāo)函數(shù)求最優(yōu)解的方法來(lái)確定節(jié)點(diǎn)的新位置，從而解決了Laplacian光順易產(chǎn)生不合法單元的問(wèn)題，提高網(wǎng)格質(zhì)量。

移動(dòng)模型表面的點(diǎn)往往會(huì)改變模型的外部輪廓，故本文采用光順的方法，只針對(duì)模型內(nèi)部的點(diǎn)。通過(guò)對(duì)度量準(zhǔn)則取倒數(shù)，得到錯(cuò)誤函數(shù)[11]。四面體網(wǎng)格中所有單元錯(cuò)誤函數(shù)之和作為基于優(yōu)化光順?biāo)惴ǖ哪繕?biāo)函數(shù)，通過(guò)對(duì)目標(biāo)函數(shù)最小值的求解，達(dá)到優(yōu)化四面體網(wǎng)格的目的。

2 數(shù)據(jù)結(jié)構(gòu)

頂點(diǎn)、邊、面和四面體是四面體網(wǎng)格的基本組成單元[12]。頂點(diǎn)包含位置信息。邊包含頂點(diǎn)索引及相鄰四面體索引。面包含頂點(diǎn)索引、面積等。四面體包含頂點(diǎn)索引、邊索引、體積及單元質(zhì)量等。頂點(diǎn)類(lèi)、邊類(lèi)、面類(lèi)及四面體類(lèi)分別存放在四個(gè)線(xiàn)性鏈表。

模型僅給出四面體的頂點(diǎn)索引以及頂點(diǎn)的位置信息，邊類(lèi)、面類(lèi)及四面體類(lèi)都需要構(gòu)造。邊類(lèi)構(gòu)造較為簡(jiǎn)單，通過(guò)遍歷四面體索引，找到該條邊，并檢測(cè)該邊是否存在，若已經(jīng)存在，則繼續(xù)下一條邊，若不存在則加入線(xiàn)性表中，并將該四面體索引添加到該邊的相鄰四面體成員中。面類(lèi)構(gòu)造與邊類(lèi)構(gòu)造類(lèi)似，不作贅述。通過(guò)四面體面類(lèi)的三個(gè)頂點(diǎn)得到兩個(gè)向量，由兩向量叉乘取模除以2易得四面體各面的面積。四面體體積可由頂點(diǎn)構(gòu)成的三個(gè)向量點(diǎn)乘除以6得到。此時(shí)，四面體單元質(zhì)量很容易算出。

3 多線(xiàn)程加速

由于模型的頂點(diǎn)數(shù)量龐大，使用并行計(jì)算的算法可充分利用多核CPU的計(jì)算能力。為了使系統(tǒng)中的多線(xiàn)程能有條不紊地運(yùn)行，須提供用于實(shí)現(xiàn)線(xiàn)程間同步和通信的機(jī)制，本文采用互斥鎖（mutex）[10]實(shí)現(xiàn)線(xiàn)程間對(duì)資源互斥訪(fǎng)問(wèn)。

所用多線(xiàn)程需要多次循環(huán)，需要設(shè)置一個(gè)計(jì)數(shù)器counter1初值為線(xiàn)程數(shù)，每個(gè)線(xiàn)程執(zhí)行完計(jì)數(shù)器減1，當(dāng)計(jì)數(shù)器為1時(shí)，主線(xiàn)程開(kāi)始執(zhí)行，并將完成對(duì)運(yùn)算結(jié)果的復(fù)制，為下次循環(huán)提供更新的數(shù)據(jù)。若此時(shí)設(shè)置counter1為初值，由于各線(xiàn)程的響應(yīng)速度不同，必將引起死鎖。因而需要設(shè)置一個(gè)障礙（barrier）來(lái)使各個(gè)線(xiàn)程在一次循環(huán)結(jié)束時(shí)做一次同步。為此加入另一個(gè)計(jì)數(shù)器counter2初值為線(xiàn)程數(shù)，當(dāng)各線(xiàn)程執(zhí)行完成后，由主線(xiàn)程設(shè)置counter1為線(xiàn)程數(shù)，開(kāi)啟下一次循環(huán)。將主線(xiàn)程視為0號(hào)線(xiàn)程，其余線(xiàn)程編號(hào)分別為1，2，3…，關(guān)鎖lock操作用于將mutex關(guān)上，開(kāi)鎖unlock操作用于打開(kāi)mutex，則主線(xiàn)程與其他線(xiàn)程所執(zhí)行的操作偽代碼如下表所示：

表2 各線(xiàn)程不同階段所執(zhí)行的偽代碼對(duì)應(yīng)表

Thread 0Thread 1，2，3…

While（counter1！=1）；Lock mutex；

--counter1；

Unlock mutex；

Copy result；

counter2=num_threads；

Lock mutex；

counter1=0；

Unlock mutex；While（counter1！=0）；

While（counter2！=1）；

Counter1=num_threads；Lock mutex；

--counter2；

Unlock mutex；

Lock mutex；

counter2=0；

Unlock mutex；While（counter2！=0）；

4 實(shí)驗(yàn)結(jié)果與分析

實(shí)驗(yàn)環(huán)境為技嘉臺(tái)式電腦，處理器是Intel（R）Core（TM）i7-3770 CPU @3.40GHz 四核，安裝內(nèi)存為8.00GB，操作系統(tǒng)是windows 7旗艦版64位，編譯器為Microsoft Visual Studio 2008。為了測(cè)試多線(xiàn)程加速情況，采用了頂點(diǎn)數(shù)量不同的模型，表3為不同模型的執(zhí)行移動(dòng)點(diǎn)操作所用時(shí)間。

表3 不同模型執(zhí)行移動(dòng)點(diǎn)操作所用時(shí)間對(duì)比表

模型（頂點(diǎn)個(gè)數(shù)）單核（秒）2核（秒）4核（秒）

7880.0034920.0026340.002330

37860.0085090.0053970.004376

51780.014710.0090120.006113

263420.136990.085610.06017

若定義單元質(zhì)量在0～0.4之間為質(zhì)量較差單元，采用本文的優(yōu)化算法，質(zhì)量較差單元所占比重由15%降低至7%，四面體質(zhì)量整體得到很大改善。從表2可以看出，各模型在多個(gè)核參與并行計(jì)算的情況下，程序運(yùn)行時(shí)間比單核都得到了不同程度的降低，這是因?yàn)椴⑿兴惴▌?chuàng)建了多個(gè)線(xiàn)程同時(shí)執(zhí)行，提高計(jì)算效率。上表中不難看出：對(duì)于頂點(diǎn)數(shù)量較少的第一個(gè)模型，多核情況下運(yùn)行時(shí)間并沒(méi)有明顯少于單核情況。這是因?yàn)槟Ｐ偷囊?guī)模不夠大，且線(xiàn)程切換與同步會(huì)有開(kāi)銷(xiāo)[10]。為了避免各個(gè)線(xiàn)程負(fù)載沉重，盡量做到計(jì)算負(fù)載平衡，筆者在執(zhí)行并行運(yùn)算前已用排序法將模型進(jìn)行優(yōu)化。

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

本文通過(guò)選用適當(dāng)?shù)亩攘繕?biāo)準(zhǔn)及優(yōu)化函數(shù)，實(shí)現(xiàn)了采用多線(xiàn)程加速方法進(jìn)行四面體網(wǎng)格優(yōu)化的目的。實(shí)驗(yàn)結(jié)果表明，在將網(wǎng)格質(zhì)量提高相同情況下，多線(xiàn)程能夠很好的提高運(yùn)算效率。

參考文獻(xiàn)：

[1]Freitag L A，Ollivier-Gooch C.A comparison of tetrahedral mesh improvement techniques[J].1996：87-100.

[2]Lo S H.Optimization of tetrahedral meshes based on element shape measures[J].Computers structures，1997（05）：951-961.

[3]Liu A，Joe B.Relationship between tetrahedron shape measures[J].BIT Numerical Mathematics，1994（02）：268-287.

[4]Zavattieri P D，Dari E A，Buscaglia G C.Optimization strategies in unstructured mesh generation[J].International Journal for Numerical Methods in Engineering，1996（12）：2055-2071.

[5]Ollivier-Gooch C F.An unstructured mesh improvement toolkit with application to mesh improvement，generation and refinement[C].AIAA 36th Aerospace Sciences Meeting.AIAA.1998（0218）：2.

[6]Freitag L A，Ollivier-Gooch C.Tetrahedral mesh improvement using swapping and smoothing[J].International Journal for Numerical Methods in Engineering，1997（21）：3979-4002.

[7]Field D A.Laplacian smoothing and Delaunay triangulations[J].Communications in applied numerical methods，1988（06）：709-712.

[8]Canann S A，Tristano J R，Staten M L.An Approach to Combined Laplacian and Optimization-Based Smoothing for Triangular，Quadrilateral，and Quad-Dominant Meshes[C]//IMR.1998：479-494.

[9]Escobar J M，Montenegro R，Montero G，et al.Smoothing and local refinement techniques for improving tetrahedral mesh quality[J].Computers structures，2005（28）：2423-2430.

[10]湯小丹，梁紅兵，哲鳳屏，湯子瀛，計(jì)算機(jī)操作系統(tǒng)（第三版），2007：71-81.

[11]董亮，劉厚林，談明高.離心泵四面體網(wǎng)格質(zhì)量衡量準(zhǔn)則及優(yōu)化算法[J].西安交通大學(xué)學(xué)報(bào)，2011（11）：100-105.

[12]賈世宇，潘振寬.虛擬手術(shù)中基于最少單元分裂的切割仿真技術(shù)[J].系統(tǒng)仿真學(xué)報(bào)，2008（06）：1487-1492.

作者簡(jiǎn)介：李密（1988-），女，山東濱州人，2011級(jí)碩士研究生，研究方向：計(jì)算機(jī)圖形學(xué)；賈世宇（1972-），男，遼寧營(yíng)口人，博士，副教授，研究方向：醫(yī)學(xué)仿真技術(shù)、計(jì)算機(jī)圖形學(xué)、計(jì)算機(jī)游戲設(shè)計(jì)；張肖醞（1990-），女，河南安陽(yáng)人，2013級(jí)碩士研究生，研究方向：計(jì)算機(jī)圖形學(xué)。

作者單位：青島大學(xué) 信息工程學(xué)院，山東青島 266071