王小慶，金先龍

(1.上海交通大學(xué) 機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240;2.上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240)

網(wǎng)格生成是整個(gè)有限元計(jì)算的重要一環(huán)。在大型工程領(lǐng)域，比如抗震分析，抗爆炸沖擊分析等，數(shù)值計(jì)算方法的應(yīng)用越來(lái)越多，對(duì)求解精度的要求也越來(lái)越高，而這些對(duì)單元網(wǎng)格的需求規(guī)模也越來(lái)越大［1－4］，傳統(tǒng)的串行網(wǎng)格生成限于時(shí)間及內(nèi)存的限制［5－6］，很難滿足超大規(guī)模數(shù)值計(jì)算的精度需要，而超級(jí)計(jì)算機(jī)的快速發(fā)展，也為網(wǎng)格生成的并行化提供了必需的硬件環(huán)境，并行網(wǎng)格劃分已經(jīng)成為網(wǎng)格生成領(lǐng)域的研究熱點(diǎn)。

并行網(wǎng)格生成的研究方法主要集中在算法并行和問(wèn)題并行兩個(gè)方面。采用問(wèn)題并行方式的研究較多，Laemmer等［7－8］提出了直接對(duì)幾何體進(jìn)行區(qū)域分解的二維三角形和三維四面體網(wǎng)格的并行生成方法，其對(duì)幾何體的分區(qū)分別采用最大和最小極慣性軸遞歸對(duì)分的方法，提供了幾何體劃分的一種新思路，但是較難達(dá)到負(fù)載的均衡。Ivanov等［9－10］提出了對(duì)初始面網(wǎng)格分區(qū)，然后再并行生成體網(wǎng)格的方法。Ito等［11］提出了對(duì)初始體網(wǎng)格進(jìn)行分區(qū)，然后利用前沿推進(jìn)法進(jìn)行并行網(wǎng)格生成的方法。在算法并行方面，Blelloch等［12－13］闡述了并行 Delaunay網(wǎng)格生成方法。Lhner等［14］對(duì)前沿推進(jìn)法的并行化進(jìn)行了研究。國(guó)內(nèi)方面，關(guān)振群等［15］總結(jié)了有限元網(wǎng)格生成方法的研究進(jìn)展，指出網(wǎng)格并行生成是未來(lái)有限元網(wǎng)格的發(fā)展方向。李水鄉(xiāng)等［16］采用問(wèn)題并行方式，實(shí)現(xiàn)了基于局域網(wǎng)的小規(guī)模有限元網(wǎng)格的并行生成。齊龍等［17］對(duì)GPU在并行網(wǎng)格生成中的應(yīng)用進(jìn)行了研究。王磊等［18］總結(jié)了各類Delaunay四面體網(wǎng)格并行生成算法的優(yōu)缺點(diǎn)，并探討了其發(fā)展趨勢(shì)。

本文在上述研究基礎(chǔ)上，提出了一種大規(guī)模非結(jié)構(gòu)化四面體網(wǎng)格并行生成方法，采用問(wèn)題并行的方式，直接讀取CAD幾何體文件，然后進(jìn)行初始網(wǎng)格劃分，并通過(guò)相對(duì)體積比及最優(yōu)分區(qū)來(lái)控制初始網(wǎng)格數(shù)量，接著采用圖劃分理論對(duì)有限元網(wǎng)格進(jìn)行區(qū)域分解，然后各個(gè)分區(qū)采用分裂法進(jìn)行并行網(wǎng)格劃分。利用初始體網(wǎng)格進(jìn)行分區(qū)克服了利用初始面網(wǎng)格進(jìn)行區(qū)域分解所帶來(lái)的負(fù)載不均衡問(wèn)題。提出的初始網(wǎng)格數(shù)量控制方法，最大限度的降低了此串行部分對(duì)并行效率的影響。針對(duì)分區(qū)間節(jié)點(diǎn)匹配的難題，提出了一種基于共享單元的邊界處理方法，有效地解決了分區(qū)間邊界協(xié)調(diào)的問(wèn)題。最后通過(guò)算例證明了整個(gè)四面體并行生成方法可以獲得良好的加速比和并行效率，所產(chǎn)生的網(wǎng)格具有較高質(zhì)量可以直接用于后續(xù)并行計(jì)算。

圖1 并行網(wǎng)格生成流程圖Fig.1 Flow Chart of parallel mesh generation

1 并行網(wǎng)格生成框架

此處所闡述的并行網(wǎng)格劃分方法，作為后續(xù)自主開發(fā)的并行求解器的前處理器的核心算法，克服了傳統(tǒng)串行網(wǎng)格劃分時(shí)間及內(nèi)存耗費(fèi)的瓶頸［19］，程序采用C++和Fortran語(yǔ)言編寫，利用MVAPICH并行庫(kù)實(shí)現(xiàn)并行通信。本文采用問(wèn)題并行方式，串行網(wǎng)格劃分部分采用前沿推進(jìn)法。

1.1 并行網(wǎng)格劃分算法

如圖1為并行網(wǎng)格劃分的流程圖。輸入文件為通用CAD幾何文件，采用前沿推進(jìn)法［20］進(jìn)行初始網(wǎng)格劃分，并控制其產(chǎn)生合理的初始網(wǎng)格數(shù);接著利用圖論理論的圖劃分方法，進(jìn)行區(qū)域分解，將各區(qū)域的數(shù)據(jù)發(fā)送到其余進(jìn)程，然后采用分裂法進(jìn)行并行網(wǎng)格生成。后續(xù)通過(guò)邊界條件處理提交并行求解器進(jìn)行求解計(jì)算，最后形成可視化后處理文件。在整個(gè)文件系統(tǒng)中，除了初始幾何文件為單一文件傳遞外，其余文件均采用多文件傳輸處理，每個(gè)核負(fù)責(zé)處理自己的文件，并行多文件的利用，在很大程度上提高了程序的執(zhí)行效率。

1.2 最優(yōu)初始網(wǎng)格生成方法

初始網(wǎng)格劃分的數(shù)量決定整個(gè)算法中串行部分的時(shí)間，劃分?jǐn)?shù)量過(guò)多，串行時(shí)間會(huì)過(guò)長(zhǎng)，因而降低并行效率，劃分?jǐn)?shù)量過(guò)少，會(huì)導(dǎo)致模型失真，網(wǎng)格質(zhì)量下降，同時(shí)所能進(jìn)行的分區(qū)數(shù)也受到限制，影響并行規(guī)模的擴(kuò)大。為合理解決此問(wèn)題，此處采用雙重約束對(duì)初始網(wǎng)格數(shù)量進(jìn)行控制和優(yōu)化。首先，保證初始網(wǎng)格的總體積和原幾何體體積相對(duì)差值不大于5%，從而保證初始網(wǎng)格的幾何貼切度，使有限元模型不失真，如式(1)為其約束條件。不同相對(duì)體積比初始網(wǎng)格效果如圖2所示，可以看到對(duì)于單位球，其初始網(wǎng)格的體積和原幾何體體積相對(duì)差值控制在4.8%時(shí)基本能保證模型不失真。其次，在保證最佳分區(qū)效果即保證分區(qū)均衡度的情況下，確定每個(gè)分區(qū)最少單元數(shù)，由于圖剖分和圖本身的拓?fù)溆泻艽笙嚓P(guān)性，因此只能通過(guò)多次試驗(yàn)的方法確定其在一定分區(qū)數(shù)時(shí)對(duì)應(yīng)的最佳單元數(shù)。

式中 a，b，c是共頂點(diǎn)的三條棱邊的長(zhǎng)度，α，β，γ 為相鄰棱邊組成的面角，ω為這三個(gè)面角和的一半。

圖2 不同體積比初始網(wǎng)格效果示意圖Fig.2 Effect of different volume ration of initial mesh

1.3 區(qū)域分解技術(shù)

有限單元網(wǎng)格的區(qū)域分解，采用圖論圖劃分方法。首先將單元信息轉(zhuǎn)換為其對(duì)應(yīng)的對(duì)偶圖，然后將對(duì)偶圖利用圖論的方法進(jìn)行區(qū)域分解，其過(guò)程分為粗化圖階段，圖分區(qū)階段及圖還原階段，最后再將分區(qū)結(jié)果映射還原為有限元網(wǎng)格的區(qū)域分解信息。圖3所示為有限元網(wǎng)格與其相應(yīng)對(duì)偶圖的相互對(duì)應(yīng)關(guān)系，每個(gè)有限元單元對(duì)應(yīng)對(duì)偶圖中的一個(gè)頂點(diǎn)，具有共享邊(面)的單元，在對(duì)偶圖中其相應(yīng)的頂點(diǎn)之間，形成一條邊。圖中頂點(diǎn)和邊上的數(shù)字代表權(quán)重，此處所有節(jié)點(diǎn)和邊按照等權(quán)重處理，大小均置為1。

圖3 有限元網(wǎng)格及其對(duì)偶圖Fig.3 Finite element model and its dual graph

有限元模型轉(zhuǎn)換形成的對(duì)偶圖，利用CSR(compressed storage)格式進(jìn)行存儲(chǔ)。CSR格式在存儲(chǔ)稀疏圖領(lǐng)域有廣泛應(yīng)用。這種格式利用兩個(gè)數(shù)組來(lái)表達(dá)圖的鄰接信息，第一個(gè)數(shù)組存儲(chǔ)第i個(gè)頂點(diǎn)相鄰頂點(diǎn)的開始和終了地址索引，第二個(gè)數(shù)組存儲(chǔ)著每個(gè)頂點(diǎn)的鄰接頂點(diǎn)，根據(jù)第一個(gè)數(shù)組中第i個(gè)頂點(diǎn)的索引可以在第二個(gè)數(shù)組中找到與其鄰接的頂點(diǎn)。程序中具體實(shí)現(xiàn)采用metis的k路圖劃分方法［21］?；趫D論理論的區(qū)域分解能夠保證負(fù)載的均衡，提高并行效率。

1.4 基于分裂法的并行網(wǎng)格生成

初始網(wǎng)格進(jìn)行區(qū)域分解后，由主進(jìn)程發(fā)往各個(gè)進(jìn)程，各從進(jìn)程接收到數(shù)據(jù)后，利用分裂法對(duì)網(wǎng)格進(jìn)行并行加密，實(shí)現(xiàn)網(wǎng)格的并行生成，直到滿足網(wǎng)格數(shù)量及精度要求［20，22］。分裂法是在原單元各個(gè)邊的中點(diǎn)生成新的節(jié)點(diǎn)，作為新增網(wǎng)格的頂點(diǎn)，然后分裂產(chǎn)生新的單元，一個(gè)三角形單元分裂為四個(gè)三角形單元，一個(gè)四面體單元分裂為8個(gè)四面體單元，其過(guò)程如圖4所示。

圖4 二維、三維單元分裂示意圖Fig.4 Schematic of 2D and 3D element split

采用分裂法很大程度上減少了各個(gè)分區(qū)之間的通信量，但是為了保證網(wǎng)格的整體性，各個(gè)分區(qū)邊界上的節(jié)點(diǎn)應(yīng)該保持一致的編號(hào)，否則無(wú)法進(jìn)行后續(xù)計(jì)算或者網(wǎng)格的合并。為此，各個(gè)分區(qū)邊界之間需要保持一定的通信，這里提出了一種基于共享單元的邊界判定方法，首先根據(jù)對(duì)偶圖區(qū)域分解結(jié)果，確定邊界四面體單元;通過(guò)遍歷所有邊界單元的邊，并在其所有共享進(jìn)程之間進(jìn)行通信，確定邊界共享邊;采用隨機(jī)算法在共享進(jìn)程間分配共享邊;劃分每個(gè)分區(qū)的邊為內(nèi)部邊、屬于自己的邊界邊、不屬于自己的邊界邊三種類型，分別插入新的節(jié)點(diǎn)，并在相應(yīng)邊界之間交換節(jié)點(diǎn)信息;產(chǎn)生新的單元;進(jìn)入下一個(gè)循環(huán)?；诠蚕韱卧姆椒?，借鑒了文獻(xiàn)［10］基于共享節(jié)點(diǎn)的方法，同時(shí)克服了其無(wú)法處理圖5情況的一個(gè)缺陷，圖示該邊的兩個(gè)節(jié)點(diǎn)96和113被兩個(gè)進(jìn)程共享，采用基于共享節(jié)點(diǎn)的方法會(huì)判定其為邊界邊，但是事實(shí)上它是內(nèi)部邊;而采用基于單元的邊界判定方法，通過(guò)對(duì)該邊界單元四條邊的遍歷，及與其共享進(jìn)程的通信，發(fā)現(xiàn)其共享進(jìn)程Core1并沒有該邊，這樣可以判定該邊為內(nèi)部邊，而非邊界邊?；诠蚕韱卧至逊ǖ木唧w實(shí)現(xiàn)流程如圖6所示。

圖5 特殊邊界棱邊示意圖Fig.5 Schematic of special edge

2 算例分析

本節(jié)通過(guò)算例分析該算法的有效性和實(shí)用性。算例測(cè)試環(huán)境均為上海超級(jí)計(jì)算機(jī)中心“魔方”曙光5000A超級(jí)計(jì)算機(jī)，整機(jī)理論計(jì)算峰值為200 Tflops，具體計(jì)算區(qū)域?yàn)镃區(qū)，該區(qū)由800個(gè)刀片式服務(wù)器組成，每個(gè)服務(wù)器包含四顆完全相同的AMD Barcelona 1.9 GHz四核處理器，每個(gè)刀片共享內(nèi)存為64 G。所有服務(wù)器操作系統(tǒng)均為64位SuSE Linux Enterprise Server 10 SP2，通過(guò)LSF作業(yè)系統(tǒng)進(jìn)行作業(yè)管理。

2.1 算例1

算例1選擇在電力、冶金、石油、化工、水利等行業(yè)有廣泛用途的滑動(dòng)式馬鞍支座，通過(guò)該算例驗(yàn)證算法的并行性能及可擴(kuò)展性。首先通過(guò)圖7直觀的展現(xiàn)采用八個(gè)計(jì)算核心時(shí)，并行網(wǎng)格劃分的具體實(shí)現(xiàn)過(guò)程，即包括幾何體讀入、初始網(wǎng)格劃分、區(qū)域分解以及并行網(wǎng)格生成等。

圖6 基于分裂法的并行網(wǎng)格生成流程圖Fig.6 Flow chart of parallel mesh generation based on the split method

通過(guò)不同規(guī)模的網(wǎng)格及不同規(guī)模的核數(shù)，對(duì)算例進(jìn)行分析，評(píng)估其時(shí)間性能，基于初始網(wǎng)格數(shù)量的控制方法，對(duì)于128核及其以下核數(shù)時(shí)，串行部分初始網(wǎng)格取877個(gè)，256核及其以上核數(shù)時(shí)，取初始網(wǎng)格取7 016個(gè)。如表1所示為不同規(guī)模不同核數(shù)時(shí)，網(wǎng)格劃總分時(shí)間?？梢钥吹剑瑢?duì)于千萬(wàn)數(shù)量網(wǎng)格，可以在幾秒內(nèi)完成網(wǎng)格的劃分;而對(duì)串行很難完成的億以及十億量級(jí)的網(wǎng)格劃分，利用該并行網(wǎng)格生成方法可以很快完成，其中對(duì)將近二十億的網(wǎng)格采用1024個(gè)處理核時(shí)，可以在大約0.28 min內(nèi)完成網(wǎng)格劃分。

圖7 八分區(qū)并行網(wǎng)格生成示意圖Fig.7 Diagram of parallel mesh generation with 8 partitions

表1 生成不同數(shù)量網(wǎng)格總時(shí)間(s)Tab.1 Time of different amounts of element generation

并行算法的可擴(kuò)展性是指隨問(wèn)題規(guī)模及計(jì)算資源的同比擴(kuò)大，并行算法效率的穩(wěn)定性。加速比是衡量可擴(kuò)展性的常用方法。如圖8為在問(wèn)題規(guī)模固定為億量級(jí)網(wǎng)格時(shí)，加速比隨核數(shù)的變化。由圖可以看到，在不考慮串行部分即初始網(wǎng)格劃分和分區(qū)時(shí)間的情況下，該算法可以獲得相當(dāng)好的加速比。在考慮串行部分后，整體的加速比有所下降，128核、256核和512核均可以達(dá)到70%的并行效率，在1 024核情況下仍然達(dá)到大約56%的并行效率。

2.2 算例2

算例2選擇具有實(shí)際工程應(yīng)用背景的某核電站防浪堤抗震分析模型進(jìn)行大規(guī)模非結(jié)構(gòu)化網(wǎng)格的并行生成，主要考察生成網(wǎng)格的質(zhì)量。整個(gè)模型長(zhǎng)度2 km、寬度950 m、高度210 m。該算例采用1024個(gè)核進(jìn)行網(wǎng)格的并行劃分，總運(yùn)行時(shí)間為100.3 s，最終的網(wǎng)格單元數(shù)2 750 414 848，節(jié)點(diǎn)數(shù)480 617 472，如圖9所示為防浪堤抗震分析模型1024分區(qū)整體網(wǎng)格生成結(jié)果及局部網(wǎng)格詳圖。

圖8 并行網(wǎng)格生成加速比Fig.8 Speedup of parallel mesh generation

圖9 防浪堤抗震模型1024分區(qū)并行網(wǎng)格劃分結(jié)果Fig.9 Parallel mesh generation result of anti-seismic breakwater model with 1024 partitions

并行網(wǎng)格劃分方法克服了串行劃分方法的時(shí)間及內(nèi)存的瓶頸，但是其網(wǎng)格質(zhì)量必須滿足后續(xù)串行或并行求解器的需要，為此，針對(duì)該算例分析其網(wǎng)格質(zhì)量，這里將網(wǎng)格質(zhì)量Q定義如式:

式中，VolIdeal=8r3/9，為與該四面體同一個(gè)外接球的理想四面體的體積，r為該四面體的外接球半徑。VolActual為每個(gè)四面體的體積，通過(guò)式(2)可求得。理想的Q值為0，即該四面體為正四面體時(shí)，網(wǎng)格質(zhì)量最高。如圖10為防浪堤抗震模型網(wǎng)格質(zhì)量的統(tǒng)計(jì)結(jié)果圖，可以看到，Q值集中在0－0.3之間，和商業(yè)軟件串行產(chǎn)生的網(wǎng)格質(zhì)量分布情況類似，可以滿足后續(xù)計(jì)算的需要。

圖10 防浪堤抗震模型網(wǎng)格質(zhì)量分布圖Fig.10 Mesh quality distribution of anti-seismic breakwater model

3 結(jié)論

(1)所提出的非結(jié)構(gòu)化四面體網(wǎng)格并行生成方法，克服了串行網(wǎng)格劃分存在的瓶頸問(wèn)題，通過(guò)算例證明了其高效性和所生成網(wǎng)格的實(shí)用性。

(2)提出的初始網(wǎng)格數(shù)量控制方法，減少了串行部分影響，有效地提高了并行效率;所述的基于共享單元的邊界處理方式，有效地解決了各個(gè)分區(qū)之間邊界節(jié)點(diǎn)的匹配，對(duì)類似邊界之間的協(xié)調(diào)問(wèn)題具有借鑒意義。

［1］鄧榮兵，金先龍，陳峻.爆炸流場(chǎng)與玻璃幕墻動(dòng)力響應(yīng)的仿真計(jì)算方法［J］.振動(dòng)與沖擊，2011，30(3):14－17.DENG Rong-bing，JIN Xian-long，CHEN Jun.Numerical simulation method for blast flow and dynamic of glass curtain wall［J］.Journal of Vibration and Shock，2011，30(3):14 －17.

［2］ Zhou M，Xie T，Seol S，et al.Tools to support mesh adaptation on massively parallel computers［J］.Engineering with Computers，2012，28(3):287－301.

［3］Jou W.Comments on the feasibility of LES for commercial airplane wings［R］.AIAA －98－2801，1998.

［4］ Baum J D，Luo H，Lhner R.Numerical simulation of blast in the world trade center［C］.AIAA，Aerospace Sciences Meeting and Exhibit，33rd，Reno，NV.1995.

［5］馬新武，王芳，趙國(guó)群.具有內(nèi)部特征約束的四邊形網(wǎng)格生成方法［J］.計(jì)算力學(xué)學(xué)報(bào)，2012，29(6):855－900.MA Xin-wu，WANGFang，ZHAO Guo-qun.An algorithm for quadrilateral mesh generation with internal feature constraints［J］.Chinese Journal of Computational Mechanics，2012，29(6):855－900.

［6］田素壘，張志毅，陳敏，等.四面體網(wǎng)格生成方法的研究與實(shí)現(xiàn)［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2012，33(11):4416－4421.TIAN Su-lei，ZHANG Zhi-yi，CHEN Min，et al.Research and implementation of methods for tetrahedral mesh generation［J］.Computer Engineering and Design，2012，33(11):4416－4421.

［7］Laemmer I L.Parallel mesh generation［M］.Solving Irregularly Structured Problems in Parallel.Springer Berlin Heidelberg，1997:1 －12.

［8］ AndrH，Gluchshenko O N，Ivanov E G，et al.Automatic parallel generation of tetrahedral grids by using a domain decomposition approach［J］.Computational Mathematics and Mathematical Physics，2008，48(8):1367 －1375.

［9］ Ivanov E G， Andr H， Kudryavtsev A N. Domain decomposition approach for automatic parallel generation of tetrahedral grids［J］.Comput.Methods Appl.Math.，2006，6(2):178－193.

［10］ Ylmaz Y，zturan C，Tosun O，et al. Parallel mesh generation、migration and partitioning for the elmer application［R］.PRACE，2010:1－12.

［11］ Ito Y，Shih A M，Erukala A K，et al.Parallel unstructured mesh generation by an advancing front method［J］.Mathematics and Computers in Simulation，2007，75(5 －6):200－209.

［12］ Blelloch G E，Miller G L，Hardwick J C，et al.Design and implementation of a practical parallel Delaunay algorithm［J］.Algorithmica，1999，24(3－4):243－269.

［13］Chrisochoides N，Nave D.Parallel Delaunay mesh generation kernel［J］.International Journal for Numerical Methods in Engineering，2003，58(2):161－176.

［14］ Lhner R.A parallel advancing front grid generation scheme［J］. International Journal for Numerical Methods in Engineering，2001，51(6):663－678.

［15］關(guān)振群，宋超，顧元憲，等.有限元網(wǎng)格生成方法研究的新進(jìn)展［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2003，15(1):1－14.GUAN Zhen-qun， SONG Chao， GU Yuan-xian， et al.Recent advances of research on finite element mesh generation methods［J］.Journal of Computer-Aided Design ＆ Computer Graphics，2003，15(1):1 －14.

［16］李水鄉(xiāng)，王云鵬，陳永強(qiáng).基于局域網(wǎng)的有限元網(wǎng)格分布式并行生成［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2005，26(12):3165－3166.LI Shui-xiang， WANG Yun-peng， CHEN Yong-qiang.Distributed parallel FEM mesh generation on LAN［J］.Computer Engineering and Design，2005，26(12):3165－3166.

［17］齊龍，肖素梅，劉云楚，等.基于GPU的并行非結(jié)構(gòu)網(wǎng)格生成技術(shù)研究［J］.機(jī)械設(shè)計(jì)與制造，2013(2):184－186.QI Long，XIAO Su-mei，LIU Yun-chu，et al.Research on parallel unstructured mesh generation technology based on GPU［J］.Machinery Design ＆ Manufacturer，2013(2):184－186.

［18］王磊，聶玉峰，李義強(qiáng).Delaunay四面體網(wǎng)格并行生成算法研究進(jìn)展［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2011，23(6):923－932.WANGLei，NIE Yu-feng，LI Yi-qiang，Advances of research on parallel Delaunay tetrahedral mesh generation ［J］.Journal of Computer-Aided Design ＆ Computer Graphics，2011，23(6):923 －932.

［19］ Larwood B G，Weatherill N P，Hassan O，et al.Domain decomposition approach for parallel unstructured meshgeneration［J］. International journal for numerical methods in engineering，2003，58(2):177－188.

［20］ Schberl J. NETGEN:an advancing front 2D/3D-mesh generator based on abstract rules［J］. Computing and visualization in science，1997，1(1):41－52.

［21］ Karypis G， Kumar V. Metis:a software package for partitioning unstructured graphs，partitioning meshes，and computing fill-reducing orderings of sparse matrices，version 5.1.0［M］.Minneapolis:Department of Computer Science and Engineering University of Minnesota，2013:26－30.

［22］Garatani K，Nakamura H，Okuda H，et al.GeoFEM:High performance parallel FEM for solid earth［C］.High-Performance Computing and Networking.Springer Berlin Heidelberg，1999:133－140.