基于有效并行求解策略的顯式有限元分析并行算法

付朝江,王天奇,林悅榮

(1.福建工程學院 土木工程學院,福州 350108; 2.福建省土木工程新技術與信息化重點實驗室,福州 350108)(*通信作者電子郵箱cjfu@163.com)

0 引言

有限元法進行結構非線性動力反應的計算非常耗時,并行計算可提供行之有效的方法。顯式時間離散可避免大量矩陣的存儲和重復求解線性方程系統,但受條件穩定性限制,較小的時間步長導致模擬需要很多時間步數,采用顯式有限元結構動力分析的并行計算方法可以在較短的時間內求解大規模問題。

利用有限元法求解結構動力問題時,最廣泛采用的顯式時間積分技術為中心差分算法,通常為顯式Newmark-β法,可表示為:

(1)

(2)

(3)

(4)

(5)

算法1 中心差分算法。

時間步循環

節點/自由度循環

1)更新位移和速度向量(式(1)和(2))

2)計算有效節點力(式(3)和(4))

3)計算加速度向量(式(5))

4)計算穩定時間步

結束節點循環

5)更新時間t=t+Δt

結束時間步循環

因為第2)步涉及材料應力更新的單元應力計算，所以第1)和3)步的計算成本相對第2)步小。采取區域分解,各個計算步均可進行并行化。

算法中的每個時間步構成一個自然的同步點,為了在t+Δt時繼續進行,所有節點和單元在t時的量需為已知。在下一個全局同步點,負載平衡使所有的處理器同時完成各自的工作任務。在顯式結構動力有限元分析中,CPU時間與單元數成比例。顯然,顯式非線性動力分析的并行實現主要挑戰在于負載平衡。

有限元計算的并行數值算法大都采取區域分解技術,將有限元網格劃分為子網格,每個子網格由單個處理器獨立處理。通過采取有效的靜態網格劃分技術來實現處理器間的負載平衡和通信最小化[1-2]。

動態負載平衡成為近幾年來并行計算研究的重要領域,已有的動態負載平衡算法主要目的是設計各種方法來減小計算過程中出現的負載不平衡問題,即應用自適應重網格化技術[3-4]。這些算法的研究主要是在同構并行計算機上開展，針對網絡機群環境的分布式多指令多數據流(Multiple Instruction Multiple Data stream, MIMD)并行計算系統開展的研究甚少。大部分動態負載平衡算法應用于以下兩步[5-7]:第1步,算法確定工作負載的數量,在鄰近的處理器中進行分配以維持各處理器中負載均衡,如Cybenko[5]提出的擴散算法和Boillat[7]的確定工作流算法。第2步是任務選定,即算法選擇有限元網格的單元或節點，該步利用處理器之間單元或節點的交換來改進現有的分解。隨后,將所選定的單元或節點在鄰近的處理器中進行任務遷移。動態負載平衡算法確定的分解質量很大程度取決于所采用的選擇過程[9]。

已有的動態負載平衡算法是針對具體的應用來處理計算過程中出現的計算不平衡而開展研究,并在專用的同構并行計算機上進行，很少有算法是針對網絡機群環境的分布式MIMD并行計算系統開展[10]。本文在網絡機群并行計算環境下，采取區域分解、動態任務分配及計算與通信重疊的優化通信方法，開展非重疊通信區域分解并行算法、重疊通信區域分解并行算法、群動態任務分配算法、動態任務分配算法及動態負載平衡算法的研究，探討基于有效并行求解策略的有限元并行算法的有效性，并對算法的性能進行了評價分析。

1 并行求解策略

上述中心差分算法可采用重疊區域分解或非重疊區域分解技術實現。采用重疊區域分解技術時,在每個處理器中位于網格劃分的界面單元被復制,節點僅屬于各自區域。這樣,這些共享單元的內力計算需要在每個處理器中復制,該計算工作非常耗時。相反,采用非重疊區域分解技術時,單元僅屬于各自區域,界面節點為各區域共享,如圖1所示。因而,界面單元力計算不需復制,而只是涉及界面節點位移、速度和加速度向量的計算需要復制。在這兩種方法中,處理器間的通信階數(即重疊區域分解中的加速度向量和非重疊區域分解中的內力向量)是相同的。因此,針對顯式非線性動力分析的并行實現,采取非重疊區域分解技術的方法較為優越。該非重疊區域分解中心差分時間步并行算法簡稱為非重疊通信區域分解并行算法(non-overl apped算法)可用算法2描述。

圖1 非重疊區域分解的網格劃分

算法2 非重疊通信區域分解并行算法。

時間步循環

節點/自由度循環

1)更新位移和速度向量(式(1)和(2))

2)計算有效節點力(式(3)和(4))

3)通過消息傳遞交換共享節點的力

4)計算加速度向量(式(5))

5)計算穩定時間步

結束節點循環

6)更新時間t=t+Δt

結束時間步循環

1.1 計算與通信重疊

在算法2描述的并行算法中, 通過將通信和計算重疊，實現通信延遲的最小化，以進一步改善算法的性能。該并行算法中,在子區域中每個單元的內力向量計算完成之后,通信才開始。然而,如果由共享節點構成的共享單元的內力計算完成,得到所需的信息(共享節點的內力向量)同相鄰節點進行交換,此時開始通信而不是等待內力計算結束才開始通信,其余單元的計算可以獨立進行,從而實現通信和計算重疊。通信和計算重疊可通過將單元力計算的任務劃分為兩部分：第一部分計算共享節點的內力向量,該階段通信發送和第二部分的計算同時進行；第二部分計算內部單元的內力向量。該優化通信的非重疊區域分解中心差分時間步并行算法簡稱為重疊通信區域分解并行算法(overlapped算法)可用算法3描述。

算法3 重疊通信區域分解并行算法。

時間步循環

節點/自由度循環

1)更新位移和速度向量(式(1)和(2))

2)計算有效節點力(式(3)和(4))

3)設立非阻塞通信以交換共享節點的力,同時進行內部節點的有效節點力的計算

4)等待處理器間通信完成

5)計算加速度向量(式(5))

6)計算穩定時間步

結束節點循環

7)更新時間t=t+Δt

結束時間步循環

1.2 動態任務分配

現有的顯式非線性動力分析的并行算法是假定整個計算過程中計算狀態保持不變,據此假設,計算負載按單元或節點數相等簡單地分配給各個處理器，利用圖劃分算法即可實現[1]；然而,這種劃分不適應材料非線性。通常塑性局限于某些單元,而這些單元需要按彈塑性情況處理,可能需要對單元每個的積分點按很復雜的材料準則進行處理；而其他不考慮塑性的單元,涉及的計算相對簡單、消耗的計算時間較少,因此，并行求解中會出現負載不平衡。本節對非線性動力分析的動態任務分配方法進行探討。

在中心差分法中,動態任務分配的并行化只在計算階段考慮[11],可利用主-從(Master-Slave)模式實現。Master(即主處理器)作為整個時間積分過程的協調者并分配任務給Slave(從處理器);Slave從Master接收任務,進行必要的計算并將計算結果發送給Master,等待Master的下一步數據。這里無需進行有限元網格劃分,位移計算在Master中進行,每個單元的內力計算分配給Slave。所有的Slave被分配任務后,Master等待Slave返回的結果。一旦接收到某個Slave的結果數據包,一個新的數據包(包含涉及另一個網格單元的數據)就會發送給Slave。這樣,任務動態分配給所有處理器,計算負載將會很好地平衡，然而在動態任務分配中的同步通信會較高。動態任務分配(Dynamic Task Allocation, DTA)算法如算法4。

算法4 動態任務分配(DTA)算法。

主處理器:

時間步循環

1)更新位移和速度向量(式(1)和(2))

2)除去所有單元標記

3)發送無標記的單元數據給空閑的從處理器

4)標記信息包已發送的單元

5)對單元進行循環

(a)接收來自從處理器的單元內力向量

(b)發送相應于無標記的單元信息包給從處理器

(c)標記相應的單元

(d)組集內部節點力向量

結束單元循環

6)計算有效節點力向量(式(3))

7)計算加速度向量(式(5))

8)計算穩定時間步

9)更新時間t=t+Δt

結束時間步循環

從處理器:

當循環沒有結束時

接收來自主處理器的單元信息包

計算單元內力向量

發送單元內力向量給主處理器

結束循環

1.3 群動態任務分配

由DTA算法可知,從處理器需要頻繁地和主處理器進行通信,從而增加了算法的通信復雜度并隨網格中單元數線性增加。為減小處理器間的通信量,提出群算法。群算法形成多個群,每個群由一組單元構成,每個群中的單元數不必相同。然而,群數比處理器數多。DTA中群形成算法(DTA with clustering算法)如算法5。

算法5 DTA中群形成算法。

1)除去所有的單元標記

2)計算網格中每個單元的度

//度表示連接到該單元的無標記單元的數目

3)對無標記的單元循環計算

選擇具有最高度的無標記單元

循環計算

如果單元的度為零,則

(a)將單元分配給相鄰的群

否則

(a)設置群數為:cluster_conter=cluster_conter+1

(b)標記單元并選擇與該單元相鄰的所有單元來組成一個群

(c)標記所有選擇的單元

(d)重算所有無標記的單元的度

結束循環

結束循環

DTA算法是對群操作而不是單元。每次涉及每個群的信息發送給空閑的從處理器,從處理器計算群中所有單元內力并組集以形成該群的整體內力向量,將其發送給主處理器并等待下一組群信息;主處理器依次組集由從處理器接收的內力向量以形成整個網格的整體內力向量,并進一步形成有效的節點力向量和計算該時間步的節點加速度。群DTA的優點就是通信復雜度減小了NE/NC,其中：NC為形成的群數、NE為網格中的單元數。附加的開銷是群的形成,該部分可在預處理階段完成。

1.4 動態負載平衡

工作站機群環境下有限元計算的負載平衡算法可按下述任一方法[12]設計：

方法一 將網格劃分成多個子網格,子網格數少于處理器數。首先,每個子網格分配給各個處理器。隨著計算進行和負載不平衡出現,調整處理器數以減小負載不平衡。

方法二 將網格劃分成子網格,子網格數大于處理器數。將一組子網格分配給每個處理器。隨著計算進行和負載不平衡出現,在處理器間調動子網格以減小負載不平衡。

方法三 將網格劃分成子網格,子網格數等于處理器數。每個子網格分配給一個處理器。隨著計算進行和負載不平衡出現,通過在子網格間轉移單元來改進初始的劃分以減小負載不平衡。

方法一為一個多級方法,會導致多個處理器在同一時間通信,導致網絡工作站的性能大大降低。方法二不會產生這種通信方式；然而,單元塊在處理器間的移動不能保證切割邊緣的最小化。方法三不會產生復雜的通信方式，因為單元遷移代替了單元塊移動,切割邊緣的最小化可以實施。由于網絡工作站并行計算要求處理器間的通信最小化以充分利用處理器資源,設計動態負載平衡算法時,方法三最優。

大多數動態負載平衡算法采用擴散算法[4]，但這些算法需要較小的迭代次數收斂,明顯地增加涉及動態負載平衡的開銷。本文采用方法三實現負載動態平衡。動態負載平衡算法構架的設計可分為下述5個階段。

1.4.1 評價處理器能力

該階段動態負載平衡器收集數據。將處理器能力CPi定義為每秒內處理器能完成的節點/單元數的計算指令。在并行算法執行期間,如果一個處理器同時正在執行各種其他工作,處理器能力隨時間變化。處理器能力確定是通過將工作負載除以處理器利用時間。工作負載Wi為處理器i的單元/節點數,處理器利用時間Ti為處理器i的運行真實時間，處理器能力CPi可表示為CPi=Wi/Ti。

1.4.2 檢查負載不平衡

該階段負載平衡器檢查是否需要重新平衡計算負載。為便于判斷,采用負載不平衡系數來評估[13],該系數定義為:LIF=Tmax/Tideal-1。其中：Tmax為處理器利用時間的最大值；Tideal為總工作負載和總處理能力的比值,即表示當工作負載完全平衡和處理器是有效利用時的理想時間。

(6)

其中:NP為處理器數,如果總工作負載是完全平衡,LIF為零;否則將會大于零。當負載不平衡系數LIF超過容許值時,將進行負載平衡初始化。

1.4.3 形成負載轉移矩陣

這個階段基于計算的負載不平衡,計算在處理器間將要轉移的工作量,即轉移的單元數，可通過形成負載轉移矩陣來實現。

依據LIF,要求負載平衡,需要形成負載轉移矩陣F,其大小為NP×NP[9]。負載轉移矩陣的每個元素Fij確定在處理器i和j間進行轉移的工作量,以實現計算負載平衡。本文負載轉移矩陣F不需顯式形成,而是計算負載轉移參數,該參數表示將要發送或從相鄰處理器接收的單元數。負載轉移參數表示為:

ψi=(Ti-Tideal) ×CPi

(7)

ψi為正表明單元數從子圖i需要發送給相鄰子網格;類似地,ψi為負表明從相鄰子網格接收單元數。根據ψi值,可容易地構成矩陣Fij。為識別和標記每個子網格i中遷移的單元,本文采用ψ設計遷移算法。

1.4.4 選擇遷移單元

根據計算的ψi值,選擇相應子網格i的單元遷移到子網格j;反之亦然。采用啟發算法實現這種遷移[14],如算法6所述。將遷移算法應用于相應子圖并改變子網格i和子網格j間的頂點。當轉移子圖的邊界頂點給它相鄰的指定子圖時,遷移算法可最小化邊緣切割。一旦遷移頂點確定,它對應的單元也就標定。

算法6 遷移算法。

1)按ψi值,確定主動的(ψi為正)、被動的(ψi為負)和中性的(ψi為零)子圖。

2)從主動子圖收集一組頂點,該頂點為界面頂點或子圖不相連部分的頂點或按相同優先次序的其他頂點。

3)標記所有頂點。

4)從主動子圖列表的組中提取一個標記的頂點。

5)通過從主動子圖i到相鄰被動子圖j交換頂點,計算邊緣切割的減小ECij。

6)計算所有相鄰被動子圖的最大ECij和標記相應的頂點及被動子圖。

7)重復第4)～7)步，直到一組中所有頂點標記完并更新最大ECij值和相應的數據結構。

8)交換標記的頂點給被動子圖j和更新界面頂點數據結構。

9)更新負載平衡參數ψi和ψj。

10)從主動頂點組中移動標記頂點。

11)標記組中所有頂點。

12)重復第4)～11)步,直到平衡參數ψi為零。

13)重復第1)～12)步,直到所有子圖的平衡參數ψ為零。

1.4.5 遷移數據

一旦指定為遷移到接收處理器的單元被確定,這些單元連同它們的數據一起轉移。為了最小化遷移開銷,所有數據結構須變換為一個流。這個流發送到接收處理器,它負責將接收的數據結構恢復到初始狀態。

算法7 動態負載平衡算法。

對時間步循環

對節點/自由度循環

1)更新位移和速度向量(式(1)和式(2))

2)計算共享節點的有效節點力(式(3))

3)建立非阻塞通信調用來交換共享節點力和同時進行計算內部點的有效節點力

4)等待處理器間通信結束

5)計算加速度向量(式(4))

6)計算穩定的時間步

結束節點循環

If(需要檢查負載不平衡)

檢查負載不平衡

If(負載不平衡為真)

計算數據遷移信息

選取將要遷移的單元

遷移選擇的單元

Endif

Endif

7)更新時間:t=t+Δt

結束時間步循環

本文提出的動態負載平衡算法應用于結構非線性動態分析的并行計算，在動態分析中，顯式動態分析算法是進行時間積分，當負載平衡在通信后進行初始化,負載評估發生在計算階段。通常,在動態分析中,在每個時間步終點之后負載平衡可初始化[15]。然而,在每個時間步終點重新平衡常常導致太多的開銷并可能引起整體計算過程速度下降。因此,負載平衡啟動通常每隔n個時間步進行一次,n的最佳值取決于負載平衡的成本和負載平衡進展的方式。

2 數值算例

本文提出的非線性動力分析并行算法的數值計算在DELL工作站機群上進行,該機群由4臺雙CPU的DELL工作站通過100.0 Mb/s以太網連接而成,共有8個CPU(2.4 GHz Xeon chips,512 KB cache),每個節點內存1.0 GB。每臺工作站都有真實的IP地址,使用消息傳遞接口(MPI)[16-17]。操作系統為Red Hat Linux9.0。

圖2(a)為1/4的圓筒形殼,L=1 000 mm,兩直邊固定。受垂直均布荷載p=2.0 kN/m2,施加方式如圖2(b)。殼體厚度2 mm,E=2.06×105MPa,泊松比v=0.3,屈服應力fy=235 MPa,質量密度為7.8×103kg/m3。進行圓筒形殼的非線性分析時，考慮其幾何非線性。

圖2 計算模型

將殼劃分為12×24結構網格,采用8節點的殼單元。通過采用時間步為0.2 ms來驗證算法的穩定性和精度。用隱式積分算法(Newmark算法)[18]和本文基于區域分解的顯式非線性動力時間積分算法(算法4)計算節點A位移時間歷程反應,計算結果如圖3所示?？梢钥闯霰疚牟⑿酗@式時間積分算法的結果與Newmark算法的結果非常吻合,表明本文并行顯式時間積分算法是有效的,即具有很好的穩定性和精度。

圖3 節點A位移時程反應

為確定動態負載平衡(DLB)算法中負載不平衡系數LIF的容許值,設定每隔25個時間步啟動一次負載平衡,LIF的容許值分別設為0.05、0.25、0.5時,進行算例計算(自由度為4 710),DLB算法計算時間如表1所示。可以看出,容許值為0.05、0.5時計算所需的時間明顯增加,算法的性能下降。由于LIF取為0.05時,該取值小,計算過程中出現負載不平衡的次數多,一旦出現負載不平衡需要進行負載平衡初始化,數據重分配,從而增加了通信開銷;而LIF取為0.5時,由于該取值較大,計算過程中出現負載不平衡次數減少,但由于LIF取值較大,一旦出現負載不平衡時,負載不平衡程度高,需要進行大量的數據重分配,數據遷移產生的通信開銷也較大。計算表明LIF的容許值為0.25時較優。

取LIF的容許值為0.25,負載平衡啟動每隔時間步n分別設為15、25、35,進行算例計算(自由度為4 710),DLB算法計算時間如表2所示。由表2可見,啟動間隔時間步為15時,因啟動負載不平衡檢查次數較多而增加時間；而啟動間隔時間步為35時,因步數較大,出現負載不平衡時,負載不平衡程度較高,進行數據重分配時數據遷移產生的通信開銷較大,致使計算時間比啟動間隔時間步為25時稍有增加。故本文算例采用啟動間隔時間步為25,LIF的容許值為0.25。

表1 不同 LIF 容許值的動態負載平衡計算時間

表2 不同啟動間隔時間步 n 的動態負載平衡計算時間

為評估算法的性能加密有限元網格以增加問題大小(自由度)，為保證解的精度采用嚴格的收斂準則，分別進行非重疊通信區域分解并行算法、重疊通信區域分解并行算法、群動態任務分配算法、動態任務分配算法及動態負載平衡算法的并行實現。為評價五種算法的并行性能,分別用1、2、4、6、8個處理器對所求問題進行求解,用并行加速比對結果進行性能分析,并行加速比定義為:speedup=用1個處理器的運行時間/用n個處理器的運行時間。

五種并行算法的加速比如圖4～5所示?？梢钥闯?隨著問題大小的增加,相應算法的加速比提高。基于動態任務分配的并行算法的性能并不理想，盡管動態任務分配算法中考慮了動態負載平衡,其加速比較小，結果并非滿意。而群動態任務分配算法性能優于動態任務分配算法,但次于區域分解算法。這是由于動態任務分配算法中處理器間的通信開銷很高,因而性能下降。在每個時間積分步,動態任務分配算法中主和從處理器間的通信總數等于網格中的單元數,因此,通信復雜度隨網格中的單元數而增加。再者,所有從處理器和主處理器不斷地進行通信,導致較高的通信延遲，因而,這些通信開銷隨著處理器數的增加而增加。

群動態任務分配算法減少了動態任務分配時處理器間的通信開銷;然而,相對于優化的區域分解算法,通信開銷仍然較高。群動態任務分配算法的性能比動態任務分配算法有所提高,這是由于處理器間的通信減少。由于每個從處理器處理一組單元,從處理器和主處理器間的通信會錯開,這樣導致較小的通信延遲。區域分解算法中將部分通信和計算重疊較非重疊通信和計算時,其性能得以提高。

動態負載平衡算法不會產生復雜的通信方式,由于單元遷移代替單元塊移動,切割邊緣的最小化可以實施,處理器間的通信開銷較少,動態負載平衡算法最優。

本文算法與Newmark算法的CPU時間(自由度為4 710)如表3所示?？梢钥闯?對相同規模的問題本文算法比隱式算法快,相應的加速比高;并且算法的計算性能隨著問題大小的增加而提高,顯示本文并行顯式時間積分算法具有良好的性能。

圖4 算法的加速比(自由度為4 710)

圖5 算法的加速比(自由度為17 322)

處理器數CPU耗時/s本文算法Newmark算法1845.90987.652563.93705.464313.30395.068179.98235.15

3 結語

對顯式非線性動力有限元分析的并行求解算法進行研究。算法采用了非重疊區域分解技術、動態任務分配和動態負載平衡,通過將通信和計算重疊,對區域分解技術中的通信進行優化。數值算例表明優化通信的區域分解算法的性能優于傳統的區域分解算法;同時也表明區域分解技術遠優于基于動態任務分配的算法。基于動態任務分配算法的性能主要受處理器間通信復雜性的影響,通過采取群方法來降低動態任務分配的通信復雜度。動態負載平衡算法由于單元遷移代替單元塊移動,可以實施切割邊緣的最小化,處理器間的通信開銷較少,動態負載平衡算法最優。數值研究表明群動態任務分配算法的性能高于動態任務分配算法,群動態任務分配算法的性能低于區域分解算法的性能,動態負載平衡算法最優。

參考文獻(References)

[1] WALSHAW C, CROSS M. Parallel optimisation algorithms for multilevel mesh partitioning[J]. Parallel Computing, 2000, 26(12): 1635-1660.

[2] DEVINE K, FLAHERTY J. Parallel adaptive hp-refinement techniques for conservation laws[J]. Applied Numerical Mathematics, 1996, 20(4): 367-386.

[3] ZHAO B, LIU Y, GOH S H. Parallel finite element analysis of seismic soil structure interaction using a PC cluster[J]. Computers and Geotechnics, 2016, 80: 167-177.

[4] HOU Y R, DU G Z. An expandable local and parallel two-grid finite element scheme[J]. Computers and Mathematics with Applications, 2016, 71(12): 2541-2556.

[5] CYBENKO G. Dynamic load balancing for distributed memory multiprocessors[J]. Journal of Parallel and Distributed Computing, 1989, 7(2): 279-301.

[7] BOILLAT J E. Load balancing and Poisson equation in a graph[J]. Concurrency and Computation: Practice and Experience, 1990, 2(4): 289-313.

[8] DIEKMANN R, FROMMER A, MONIEN B. Efficient schemes for nearest neighbor load balancing[J]. Parallel Computing, 1999, 25(7): 789-812.

[9] SZIVERI J, TOPPING BHV. Transient dynamic nonlinear analysis using MIMD computer architectures[J]. Journal of Computing in Civil Engineering, 2000, 14(2): 79-91.

[10] HU Y F, BLAKE R J. An improved diffusion algorithm for dynamic load balancing[J]. Parallel Computing, 1999, 25(4): 417-444.

[11] TOUHEED N, SELWOOD P, JIMACK P K. A comparison of some dynamic load-balancing algorithms for a parallel adaptive flow solver[J]. Parallel Computing, 2000, 26(12): 1535-1554.

[12] KRYSL P, BITTNAR Z. Parallel explicit finite element solid dynamics with domain decomposition and message passing: dual partitioning scalability[J]. Computers and Structures, 2001, 79(3): 345-360.

[13] RAO A R M. Explicit nonlinear dynamic finite element analysis on homogeneous/heterogeneous parallel computing environment[J]. Advances in Engineering Software, 2006, 37(11): 701-720.

[14] RAO A R M, RAO T, DATTAGURU B. Generating optimised partitions for parallel finite element computations employing float-encoded genetic algorithms[J]. Computer Modeling in Engineering and Sciences, 2004, 5(3): 213-234.

[15] KWAK J Y, CHUN T Y, SHIN S J, et al. Domain decomposition approach to flexible multibody dynamics simulation[J]. Computational Mechanics, 2014, 53(1): 147-158.

[16] 付朝江. 集群MPI環境下有限元結構分析并行計算研究[D]. 上海: 上海大學, 2006.(FU C J. The research on parallel computation of finite element structural analysis based on MPI cluster[D]. Shanghai: Shanghai University, 2006.)

[17] CAPODAGLIO G, AULISA E. A particle tracking algorithm for parallel finite element applications[J]. Computers and Fluids, 2017, 159: 338-355.

[18] 付朝江. 隱式非線性動力分析有限元并行求解格式[J]. 工程力學, 2010, 27(10): 27-33.(FU C J. Parallel algorithm for implicit nonlinear dynamic finite element analysis[J]. Engineering Mechanics, 2010, 27(10): 27-33.)

This work is partially supported by the National Natural Science Foundation of Chian (51378124).