作者：咼嘉妮　胡久鄉　盧正 來源：CAD世界網

摘　要 　提出了一種基于網格生成遞歸法的并行區域劃分算法，該算法依據網格生成代價的估算分析，采用迭代分解法對區域進行并行劃分.在曙光1000A 系統上的運行結果表明，該網格算法的效率和加速比均優于串行遞歸算法.
關鍵詞　有限元網格；并行區域劃分算法；網格生成代價；迭代分解法

　　基于網絡生成遞歸法^［1～3］，本文提出了一種并行區域劃分算法，該算法滿足以下四個基本原則：a. 任務平衡性原則.能生成平衡的子區域集，即在各子區域中生成網格的時間大致相等.b. 邊界最簡原則.子區域的邊界結構簡單，邊界處理所需時間短，處理器間消息傳遞的費用低.c. 網格均勻原則.并行生成的最終網格形狀均勻，無奇異單元.d. 區域劃分代價最小原則.區域劃分算法本身的代價盡可能小.

1、基本思想及相關算法
　　在網格生成遞歸法中，如果每個子區域都包含相同的單元數，就比較容易實現任務平衡.因此，首先按照單元數估算待處理區域的網格生成代價，然后根據當前參與并行處理的處理器數N對區域進行分解，并對分解所得子區域進行邊界處理，最終獲得相互之間既平衡又獨立的N個并行子任務.
1.1　網格生成代價的估算算法
　　網格生成代價與分布于待處理區域中的單元數目緊密相關，而單元數目是由該區域的總面積S和區域內單元分布密度決定的.估算公式如下：

G=S/S_tri,?。?）

S_tri=［L²/(2M²)］sin60°,　?。?）

M=Σ2l_i/(d_i1+d_i2) (i=1, 2, …， n),　　　（3）

式中，G表示三維平面區域 (含有n條邊) 的網格生成代價；S表示該區域的面積；Stri 是對區域內單個三角形單元平均面積的粗略估算值；L表示區域的周長；M表示區域邊界上總的離散段數；l_i表示第i條邊的長度;d_i1和d_i2 分別表示第i條邊前端點和后端點處的單元邊長.
　　據以上公式，網格生成代價的估算算法步驟如下.
　　步驟 1　依次計算區域各邊界的長度和離散段數,根據式(2)和(3)計算S_tri；
　　步驟 2　應用三角累加算法計算S,根據式(1)求得G.三角累加算法步驟如下：a. 置變量S為0；b. 選取區域的任意一個頂點作為V₁,按順時針 (或逆時針) 方向取與V₁ 依次相鄰的兩個點為V₂和V₃;c. 計算由V₁, V₂ 和 V₃ 組成的三角形的面積，累加入變量S;d. 以V₃為新的V₂，按順時針 (或逆時針) 方向取與V₂相鄰的第一個點為新的V₃，若V₃=V₁，則算法停止；否則轉 c.
1.2　區域劃分算法
　　估算出整個區域的網格生成代價后，區域劃分算法的任務就是：尋找N-1條分割線，將區域分割為N個子區域，使得各個子區域的網格生成代價大致相等.由于無法精確確定子區域中的單元數目，因此允許各子區域的網格生成代價在G/N附近有一個δ誤差，亦即所產生的子區域的網格生成代價都屬于［G/N-δ, G/N+δ］區間.
　　根據原則 b，任何一條分割線的兩端點都定位在區域的邊界上，而不落在區域中.在圖 1 （a） 中，任何一個子域都只通過一條完整的分割線與其他子域相鄰；而在圖1（b） 中，子域3與子域1和2之間的相鄰關系不便于邊界處理.

　　　　　　(a)　　　　　　　　　(b)
圖 1　區域劃分效果圖

　　根據原則c，多條分割線不交于同一點，因為這種分割在N較大的情況下會導致極點現象——在多條分割線的交點處出現小內角的長薄單元，影響網格形狀.
　　本文采用迭代分解法劃分區域，步驟如下.
　　步驟 1　記待分解的區域為 D;
　　步驟 2　循環執行以下步驟N-2次:a. 在D中找一條分割線將區域分為兩個子域D₁和D₂，使D₁的網格生成代價約等于G/N;b. 將D₁從D中去掉，以D₂為新的D，轉步驟 a;
　　步驟 3　在D中找一條分割線將區域分為兩個子域D₁和D₂，使D₁和D₂的網格生成代價均屬于 ［G/N-δ,G/N+δ］.
　　本算法中，尋找分割線的方法是：先在區域中找到一個合適的頂點作為分割線的前端點X，然后從這個頂點出發，沿順時針（或逆時針)方向按對半查找法在區域的邊界上尋找分割線的后端點Y.對半查找算法步驟如下：
　　步驟 1　從X點出發，按順時針（或逆時針）方向取與X依次相鄰的兩個點為R和Y；
　　步驟 2　計算分割線XY左側 （或右側） 子域的網格生成代價g;
　　步驟 3　若g∈［G/N-δ,G/N+δ］，則算法停止；若g>G/N+δ，則令T=Y，以T與R的中點為新的Y，轉步驟2；若g<G/N-δ,則令R=Y.若Y是原區域的頂點，則沿順時針 （或逆時針） 方向取與Y相鄰的第一個點為新的Y，否則，以R與T的中點為新的Y，轉步驟2.
　　根據以下原則確定分割線的前端點X：對第一次分割，選取區域中最大頂角所在頂點為X；對第 i （i=2, 3, …, N-1） 次分割，選取分割線兩端點中較大頂角所在的頂點為X.根據這一策略，每次分割都選取頂角較大的點作為前端點，一定程度地避免了極點的形成，但還不能完全避免極點的產生.對此，引入分割線共點控制機制：首先定義參數R，由它指明能交于一點的分割線的最大數目,然后，每次分割都用變量r記錄X處的共點率，一旦發現r>R，就變換X.
　　圖 2 是對區域ABCDEF作6次劃分的示意圖.在圖2（a）中，R取值為2，尋找第5條分割線時，因為∠CHI比∠HIE大,X原應定位在H點，但由于已有兩條分割線經過H點，因此如果這次還選取H作為前端點， 那么共點率r增為3,超過預定義的R值,于是，改選I點為X；而在圖2(b)中，因為R取值為3，所以，H仍可作為第5條分割線的前端點.

　　　　　　(a) R=2　　　　　　　　　(b) R=3
圖 2　共點控制

　　該區域劃分算法的特點是：并行任務的負載平衡程度由δ調節；所有分割點都定位在區域的邊界上，因為每次分割剩下的子域只有一條邊是在分割過程中新加入的分割線，無論選取這條邊的哪個端點作為X，X和Y總是落在區域的邊界上，有效地滿足了邊界最簡原則,；通過R靈活控制分割線共點率，消除極點，保證網格均勻度.

2、并行實現
　　本文研究的并行區域劃分算法已在曙光1000A系統上實現.根據曙光1000A系統目前的條件和特點，選取PVM作為并行編程環境，分三個步驟實現有限元網格的并行生成：首先，由node₁上的master任務程序檢測當前PVM中的結點數N，運用區域劃分算法分解待處理的區域為N個子區域，將各子區域分配到各node上；然后，并行執行各處理器上的slave任務程序，在各子區域中生成均勻網格；最后，由node₁上的master任務程序收集、組裝在各處理器上生成的子區域網格.這里，采用了動態負載平衡策略——farm模式，但有一點與通常情況下不同：為了提高處理器利用率，node₁在子區域分配過程中也分得一個子區域，所以，在數據分配之后，子區域網格傳回之前，主處理器也不空閑.
　　由于在并行劃分階段已經對邊界進行了預處理，因此，一旦子區域分配給各個處理器，所有處理器將同時獨立地在自己的數據子塊上根據預定義的單元分布密度生成網格.經過大致相等的時間段后，它們將完成各自的工作，并將生成的子區域網格傳送給master任務程序.整個并行算法的通信很小，僅包括從node₁向其他各node廣播子區域數據，以及各node向node₁回送生成的子區域網格單元.

3、實驗結果
　　表 1 以實例說明上述并行算法在曙光1000A系統上的運行效果 （表中括號內外的數據分別對應例 1 和例 2）.在例 1 中，δ＝G/(5N(N-1))，G=15　987 個,實際生成的單元總數為12　531個；例 2 中，δ＝G/(5N(N-1)),G＝33　452 個，實際生成的單元總數為28　769個.由于區域劃分和組裝數據所耗費的時間是μs級的，因此將并行算法的執行時間 (t) 分為兩部分：子區域網格生成時間 (t₁) 和通信時間 (t₂).t₁ 指各node在其分配到的子區域中生成網格的時間；t₂包括node₁向其他各node發送子區域數據的時間，以及各node向node₁回送子區域網格數據的時間.表 1 給出了使用2,3,4,5,6個node進行并行網格生成的時間開銷及加速比 (ν).使用一個處理器時，既不需要做區域劃分，也不用通信，因此，在區域中生成網格的時間即為全部運行時間.

表 1　并行算法執行時間及加速比