陳禮杰,吳慧,李鐵瑞,高博青
?
基于曲面擬合的復(fù)雜自由曲面網(wǎng)格劃分
陳禮杰1,吳慧2,李鐵瑞1,高博青1
(1. 浙江大學(xué) 建筑工程學(xué)院,浙江 杭州,310058;2. 浙江財經(jīng)大學(xué) 公共管理學(xué)院,浙江 杭州,310018)
為了實現(xiàn)復(fù)雜自由曲面的建筑網(wǎng)格劃分,將復(fù)雜自由曲面擬合成1個曲面,并利用映射法進行網(wǎng)格劃分。首先提出基于非均勻有理B樣條(NURBS)的線面求交和基于NURBS的面面求交2種方法,獲得復(fù)雜自由曲面的×點云;將×點云作為控制點構(gòu)造NURBS曲面,然后把修剪后的NURBS曲面邊界按給定數(shù)目或給定長度劃分,連接參數(shù)域上的邊界節(jié)點后得到參數(shù)域網(wǎng)格,將其映射回空間曲面。對邊界處網(wǎng)格進行投影處理以減少曲面擬合出現(xiàn)的邊界褶皺及內(nèi)縮問題。采用加密點陣技術(shù)減小曲面擬合的誤差。研究結(jié)果表明:與多曲面網(wǎng)格劃分方法相比,采用曲面擬合及映射法能夠得到線條流暢、大小基本一致的網(wǎng)格,具有良好的建筑美學(xué)效果。
非均勻有理B樣條(NURBS);復(fù)雜自由曲面;曲面擬合;映射法;網(wǎng)格劃分
隨著計算機科學(xué)技術(shù)的進步、建筑施工工藝的發(fā)展和人們建筑審美水平的提高,自由曲面結(jié)構(gòu)以其優(yōu)美的造型越來越多地出現(xiàn)在人們的視野當中?!白杂汕妗敝笩o法用單個或幾個解析函數(shù)表達的曲面,也可理解為那些明顯區(qū)別于傳統(tǒng)建筑造型的曲面[1]。由于解析曲面具有一定的周期性或規(guī)律性,其網(wǎng)格劃分比較方便。而自由曲面的不規(guī)則性使傳統(tǒng)的網(wǎng)格劃分方法難以適用于自由曲面結(jié)構(gòu)。因此,如何在自由曲面上生成線條流暢、大小基本一致的網(wǎng)格,是空間結(jié)構(gòu)領(lǐng)域的熱點和難題[2]。網(wǎng)格劃分起源于有限元領(lǐng)域[3],常用的方法有Delaunay法[4]、波前法[5]、映射法[6]等。但有限元領(lǐng)域的網(wǎng)格劃分方法難以滿足建筑美學(xué)線條流暢、大小基本一致的要求。因此,近年來有研究者針對自由曲面的建筑網(wǎng)格劃分進行了相關(guān)研究,如:WINSLOW等[7]以桿件走向為優(yōu)化變量,以結(jié)構(gòu)性能為目標函數(shù),用遺傳算法實現(xiàn)自由曲面的多目標優(yōu)化;SU等[8]繪出自由曲面的主應(yīng)力跡線,采用改進的波前法自動生成網(wǎng)格;潘煒等[9]按照展開前、后面積變化最小的原則展開曲面,并利用線推進法生成平面網(wǎng)格,最后映射回空間曲面;丁慧等[10]提取等參線并以桿件長度相等為條件分割等參線,連接等參線上的節(jié)點生成網(wǎng)格。復(fù)雜自由曲面由多個非均勻有理B樣條(NURBS)裁剪曲面組成,且存在內(nèi)部空洞等復(fù)雜情況。針對該類自由曲面,現(xiàn)有的網(wǎng)格劃分方法是在每個裁剪曲面上分別劃分網(wǎng)格,然后對交界附近網(wǎng)格進行調(diào)整。該方法的缺點是網(wǎng)格劃分不易進行,交界處網(wǎng)格不好調(diào)整,難以達到建筑美學(xué)線條流暢的要求。本文作者提出基于NURBS的線面求交和面面求交2種方法,獲得復(fù)雜自由曲面的×點云。第1種方法是構(gòu)造矩形×平面點陣,過每個點作平行于Z坐標軸的射線與NURBS曲面求交,得到×點陣;第2種方法是根據(jù)NURBS曲面的特征,構(gòu)造個平面并與NURBS曲面求交,對交線均勻采樣個點,修補邊界后得到點陣。利用向心參數(shù)化方法得到點陣的和向節(jié)點矢量,構(gòu)造以為控制點陣的NURBS曲面。將修剪后的NURBS曲面邊界按給定數(shù)目或給定長度劃分,連接參數(shù)域上的邊界節(jié)點后形成參數(shù)域網(wǎng)格,將其映射回空間曲面,對邊界處網(wǎng)格作投影處理,即可得到線條流暢、大小基本一致的四邊形網(wǎng)格。
曲面擬合技術(shù)在逆向工程和計算機輔助設(shè)計等領(lǐng)域中得到了廣泛應(yīng)用[11]。利用計算機視覺技術(shù)(如激光掃描)對自由曲面進行數(shù)字化,可得到大量散亂數(shù)據(jù)點,根據(jù)數(shù)據(jù)點集可得到實物的CAD模型[12]。
在建筑工程領(lǐng)域,逆向工程用于古建測繪和古建保護、快速建立城市模型、空間網(wǎng)格結(jié)構(gòu)的逆向建模等方面。李娜[13]借鑒逆向工程的操作原理,利用激光三維掃描系統(tǒng)進行實物樣件的點云數(shù)據(jù)采集,采用蒙皮法進行曲面造型,然后進行網(wǎng)格劃分。該技術(shù)基于對實物進行激光掃描得到大規(guī)模散亂數(shù)據(jù)點。在復(fù)雜自由曲面網(wǎng)格劃分時,一般建筑師已經(jīng)提供符合其設(shè)計意圖的曲面CAD模型;為了獲得線條流暢、大小基本一致的網(wǎng)格,本文作者提出基于NURBS的線面求交和面面求交2種方法,獲得復(fù)雜自由曲面的×點云,構(gòu)造以×點云為控制點陣的NURBS曲面,為映射法劃分網(wǎng)格提供基礎(chǔ)。
圖1所示為某工程CAD模型。下面以圖1為例,介紹基于NURBS的線面求交法和基于NURBS的面面求交法。
圖1 某工程CAD模型
根據(jù)NURBS曲面的強凸包性[14]和控制點陣為矩形拓撲[12]的特點,構(gòu)造1個×矩形平面點陣作射線與NURBS曲面求交,提取原曲面信息。算法流程如圖2所示。
圖2 基于NURBS的線面求交法流程圖
Fig. 2 Flow chart of NURBS based ray-surface intersection method
圖3 延伸曲面上邊界
圖4 N×M矩形點陣及曲面包圍盒
NURBS曲面與射線求交是一種基礎(chǔ)的幾何算法,常用于雷達散射截面計算(RCS)的射線追蹤計算[15]、逆向工程中曲面的制造誤差計算[16]等方面。
首先利用Krawczyk算子判斷式(1)在某一區(qū)間是否存在解,然后對存在唯一解的區(qū)間給定1個初值,利用牛頓迭代法求解。詳細流程參見文獻[16]。
圖5 裁剪NURBS曲面的有效交點與無效交點
前面提出一種提取原曲面信息的方法,現(xiàn)通過構(gòu)造個平面與NURBS曲面求交提取原曲面信息。算法流程如圖6所示。
1) 補齊曲面空洞。當復(fù)雜自由曲面內(nèi)部存在空洞時,為了保證平面與曲面的交線在平面上能夠擬合為一條曲線,用平面補齊空洞。補齊空洞及平面切割方式如圖7所示。
圖6 基于NURBS的面面求交法流程圖
圖7 補齊空洞及平面切割方式
2) 選擇平面切割方式。復(fù)雜自由曲面由多個NURBS裁剪曲面組成,應(yīng)根據(jù)曲面各部分的走向和特征,選擇合適的切割方式,如直線式和圓弧式切割(見圖7)。原曲面由11個NURBS曲面組成(不包括補齊空洞的平面),分為5部分切割。第1部分由平面上點1和點2及平面數(shù)1確定。首先過點1和點2作直線,然后以為正法線,等距離地作1個平面與曲面求交。第2部分由平面上點3、點4、點5及平面數(shù)2確定。首先過點4以直線為正法線作平面,然后以點4為旋轉(zhuǎn)中心,軸為旋轉(zhuǎn)軸,將平面旋轉(zhuǎn)2?1次,每次旋轉(zhuǎn)角度為∠345/(2?1),即可得到2個平面與曲面相交。其余3個部分依此類推。
3) NURBS曲面與平面求交。用平面與CAD模型的個曲面求交,如圖8所示。平面與4個曲面相交,共有5條交線,其中與曲面2有2條交線即交線2 和4。
4) 處理同一平面上的交線(見圖8)。首先處理裁剪曲面與平面求交的情況。曲面1是裁剪曲面,其刪除部分不顯示,因此,把位于曲面1刪除部分且與平面相交的交線刪除。對曲面4進行同樣處理。平面與NURBS曲面共有5條交線,每條交線均勻采樣一定數(shù)目點,然后構(gòu)造1條插值于所有采樣點的NURBS曲線。每個平面都得到1條與原曲面相交的曲線。
圖8 同一平面上交線處理
5) 擬合曲線采樣。曲面擬合中NURBS曲面的型值點或控制點必須為矩形拓撲,故得到平面和原曲面擬合后的交線后,對該交線均勻采樣個點,由此保證NURBS曲面與個平面求交后得到的點陣為矩形拓撲。
6) 補齊邊界。對于與NURBS曲面相交的第1個和最后1個平面(下面稱邊界平面),總存在一部分曲面在邊界平面之外,其大小不能忽略。取這部分曲面的邊界線,并均勻采樣個點,可以保證完整地提取原曲面的信息。補齊邊界方式如圖9所示。
圖9 補齊邊界
經(jīng)過以上步驟,可以得到矩形拓撲點陣。
利用曲面擬合中常用的向心參數(shù)化[14]方法構(gòu)造矩形拓撲點陣的節(jié)點矢量和。獲得點陣和節(jié)點矢量和之后,進行曲面擬合。
對于矩形拓撲點陣,傳統(tǒng)的曲面擬合方法是將其作為型值點陣,創(chuàng)建1個非有理的(,)次B樣條曲面(,),使其插值于這些點。
式中:Q為型值點;u和v分別為其對應(yīng)的和向參數(shù)值;R(u, v)為其對應(yīng)的有理基函數(shù);P為待求的控制點。
另一種方法用1個非有理的(,)次B樣條曲面(,)逼近這些點,以控制點為優(yōu)化變量,滿足如下條件:
為了準確地表示復(fù)雜自由曲面,其矩形拓撲點陣的點數(shù)×不能太少。若使用插值或逼近方法,則至少需要解1個(×)×(×)的線性方程組,這將使計算時間大大增加。
針對復(fù)雜自由曲面擬合問題,采用控制點逼近方法。控制多面體是對B樣條曲面的平面片逼近,次數(shù)越低,逼近效果越好;當節(jié)點矢量無限加細時,控制多面體將收斂于曲面[17]。因此,以為控制點陣,在節(jié)點矢量和上構(gòu)造非有理的雙二次B樣條曲面(,),根據(jù)CAD模型的原有邊界情況對(,)進行裁剪?;贜URBS的線面求交法和面面求交法的擬合曲面與裁剪曲面分別如圖10和圖11所示,2種求交方法的擬合曲面雖不盡相同,但裁剪曲面都能反應(yīng)復(fù)雜自由曲面的特征,可為映射法網(wǎng)格劃分提供基礎(chǔ)。
(a) 擬合曲面;(b) 裁剪曲面
(a) 擬合曲面;(b) 裁剪曲面
圖12所示為映射法劃分網(wǎng)格示意圖。建筑師提供的CAD模型由3個曲面構(gòu)成(見圖12(a)),利用基于NURBS的線面求交法,將3個NURBS裁剪曲面擬合為1個曲面(見圖12(b)),利用映射法劃分網(wǎng)格。
(a) 某工程CAD模型;(b) 邊界等長度分段;(c) 空間曲面上網(wǎng)格;(d) 處理后網(wǎng)格
1) 將空間曲面邊界分段。以二分法為基礎(chǔ),將空間曲面的邊界按弦長分段。分段方式有給定長度、給定段數(shù)、等比例或指定長度數(shù)組等(見圖12(b))。
2) 將分段點投影回參數(shù)域,連接參數(shù)域上的分段點后形成參數(shù)域網(wǎng)格。分段后的邊界點反向映射回參數(shù)域中,在參數(shù)域中將邊界點按照建筑師期望的走向依次連直線。
3) 將參數(shù)域網(wǎng)格映射回空間曲面,得到空間網(wǎng)格(見圖12(c))。
4) 處理邊界處網(wǎng)格。由圖12(c)可知:用映射法在擬合曲面上畫出的網(wǎng)格在原有曲面交界處線條流暢,但在擬合曲面的邊界處質(zhì)量不佳,為此需處理邊界處網(wǎng)格。下面介紹該具體方法。
前面以圖1為例詳細介紹了基于NURBS的線面求交法(以下簡稱方法一)及基于NURBS的面面求交法(以下簡稱方法二)的算法流程。表1所示為2種方法的比較結(jié)果(其中,為將復(fù)雜自由曲面擬合為1個曲面所需的時間,U為擬合曲面上邊界最大誤差,L為擬合曲面下邊界最大誤差)。
方法一耗時最久的步驟是曲面與直線求交,需進行49萬(700×700)次求交運算;方法二耗時最久的步驟是曲面與平面求交,只需進行900次求交運算。故方法二的速度遠快于方法一的速度。
表1 2種方法的比較
圖13所示為方法一邊界處褶皺。方法一中,當曲面邊界走向不嚴格為平面點陣走向時,邊界附近的控制點不能很好地還原曲面的邊界信息,出現(xiàn)“褶皺”問題,導(dǎo)致上、下邊界最大誤差較大。此時,可以通過加密平面點陣的方法來減小褶皺,但會使時間大大 增加。
圖14所示為方法二下邊界誤差。原曲面下邊界出現(xiàn)直角轉(zhuǎn)折,在下邊界附近的3個采樣點由此也產(chǎn)生較大轉(zhuǎn)折(見方法二步驟5))。因為NURBS曲面的強凸包性,擬合曲面在此處出現(xiàn)最大誤差為250.5 mm,可以通過加密采樣點的方法減小誤差。而上邊界處轉(zhuǎn)折較小,最大誤差僅為2.1 mm。
圖13 方法一邊界處褶皺
圖14 方法二下邊界誤差
2種方法除邊界處出現(xiàn)誤差外,其余部分誤差均小于5 mm。用映射法網(wǎng)格劃分時需將邊界處網(wǎng)格投影至原曲面以消除誤差,方法一還需手動處理因“褶皺”而畸形的網(wǎng)格。
通過2種方法的比較可以得出以下結(jié)論:1) 方法二比方法一復(fù)雜;2) 方法二的速度遠快于方法一的速度,最大誤差也比方法一的??;3) 利用映射法劃分網(wǎng)格時,2種方法均需將邊界處網(wǎng)格投影回原曲面,方法一還需手動處理因“褶皺”而畸形的網(wǎng)格。
以某工程CAD模型(見圖12(a))為例,用分曲面的方法劃分網(wǎng)格,如圖15所示。首先將曲面1的邊界1等長度分段,根據(jù)欲劃分的網(wǎng)格形式確定邊界2的分段點,最后連接邊界1分段點或同時連接邊界1、邊界2分段點(見圖15(a))。劃分好曲面1的網(wǎng)格后,憑經(jīng)驗確定曲面2上邊界3的分段點,連接邊界2分段點或邊界2、邊界3分段點(見圖15(b))。同理畫出曲面3的網(wǎng)格(見圖15(c))。但上述過程依賴于工程師的經(jīng)驗,且網(wǎng)格劃分容易失敗。
在邊界2附近,曲線1為邊界1分段點相連所得曲線,曲線2為邊界1和2分段點相連所得曲線,因此,導(dǎo)致網(wǎng)格大小不均勻(見圖15(a))。由于分別在3個曲面上劃分網(wǎng)格,導(dǎo)致邊界處網(wǎng)格不均勻,線條不流暢(見圖15(c))。對邊界處質(zhì)量不佳的曲線進行一定調(diào)整,調(diào)整后網(wǎng)格均勻性有一定提高,但交界處曲線流暢性仍然不佳(見圖15(d))。
(a) 曲面1;(b) 曲面1和2;(c) 曲面1,2和3;(d) 調(diào)整后網(wǎng)格劃分
現(xiàn)用曲面擬合的方法對圖12(a)進行網(wǎng)格劃分。該工程平面投影為橢圓形,平面面積為146.0 m× 116.0 m。用方法一獲得600×600點陣,擬合為1個曲面,等分外邊界,在參數(shù)域上將外邊界分段點依次連為直線后將參數(shù)域網(wǎng)格映射回空間曲面(見圖12(c))。將邊界處網(wǎng)格投影回原曲面,并處理極少數(shù)因“褶皺”而畸形的網(wǎng)格,得到最終的網(wǎng)格劃分(見圖12(d))。
圖15中在原有的3個曲面上分別劃分網(wǎng)格,因此,在交界附近網(wǎng)格不均勻,線條過渡不流暢,且難以調(diào)整,網(wǎng)格劃分過程依賴工程師個人經(jīng)驗,且容易失??;圖12(d)中線條流暢,網(wǎng)格大小基本一致,且過渡均勻。雖然在邊界處出現(xiàn)一定誤差,但可以通過加密點陣減小誤差,或?qū)⑦吔缣幘W(wǎng)格投影回原曲面,以消除邊界處誤差。
1) 提出2種點云生成方法:方法1中構(gòu)造1個×矩形平面點陣并作射線與復(fù)雜自由曲面求交,從而獲得矩形拓撲點陣;方法2需合理選擇平面切割方式,將復(fù)雜自由曲面分為幾個部分,構(gòu)造個平面與復(fù)雜自由曲面求交,并將多條交線擬合后均勻采樣個點,修補邊界后獲得矩形拓撲點陣。
2) 將二者進行比較,令點云為控制點陣,把復(fù)雜自由曲面擬合為1個NRUBS曲面,以映射法劃分空間網(wǎng)格。
3) 與分曲面劃分方法相比,曲面擬合和映射法能夠得到線條流暢、大小基本一致的四邊形網(wǎng)格。
[1] 沈世釗, 武岳. 結(jié)構(gòu)形態(tài)學(xué)與現(xiàn)代空間結(jié)構(gòu)[J]. 建筑結(jié)構(gòu)學(xué)報, 2014, 35(4): 1?10. SHEN Shizhao, WU Yue. Structural morphology and modern space structures [J]. Journal of Building Structures, 2014, 35(4): 1?10.
[2] 李承銘, 盧旦. 自由曲面單層網(wǎng)格的智能布局設(shè)計研究[J]. 土木工程學(xué)報, 2011, 44(3): 1?7. LI Chengming, LU Dan. Study of intelligent layout design of single-layer lattice shell of free form surface[J]. China Civil Engineering Journal, 2011, 44(3): 1?7.
[3] LO S H. Finite element mesh generation and adaptive meshing[J]. Progress in Structural Engineering & Materials, 2002, 4(4): 381?399.
[4] SI Hang. TetGen, a delaunay-based quality tetrahedral mesh generator[J]. ACM Transactions on Mathematical Software, 2015, 41(2): 11.
[5] PIRZADEH S. Structured background grids for generation of unstructured grids by advancing-front method[J]. AIAA Journal, 2012, 31(2): 285?289.
[6] LISEIKIN V. Grid generation methods[M]. Amsterdam, Netherlands: Springer, 2010: 10?14.
[7] WINSLOW P, PELLEGRINO S, SHARMA S B. Multi-objective optimization of free-form grid structures[J]. Structural & Multidisciplinary Optimization, 2010, 40(1/2/3/4/5/6): 257?269.
[8] SU Liang, ZHU Shunlai, XIAO Nan, et al. An automatic grid generation approach over free-form surface for architectural design [J]. Journal of Central South University, 2014, 21(6): 2444?2453.
[9] 潘煒, 吳慧, 李鐵瑞, 等. 基于曲面展開的自由曲面網(wǎng)格劃分[J]. 浙江大學(xué)學(xué)報(工學(xué)版), 2016, 50(10): 1973?1979. PAN Wei, WU Hui, LI Tierui, et al. Grid generation on free-form surface based on surface flattening[J]. Journal of Zhejiang University: Engineering Science, 2016, 50(10): 1973?1979.
[10] 丁慧, 羅堯治. 自由形態(tài)網(wǎng)殼結(jié)構(gòu)網(wǎng)格生成的等參線分割法[J]. 浙江大學(xué)學(xué)報(工學(xué)版), 2014, 48(10): 1795?1801. DING Hui, LUO Yaozhi. Isoparametric line dividing grid generation method for free-form lattice shells[J]. Journal of Zhejiang University: Engineering Science, 2014, 48(10): 1795?1801.
[11] BRUJIC D, AINSWORTH I. Fast and accurate NURBS fitting for reverse engineering[J]. The International Journal of Advanced Manufacturing Technology, 2011, 54(5): 691?700.
[12] 來新民, 黃田. 基于NURBS的散亂數(shù)據(jù)點自由曲面重構(gòu)[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 1999, 11(5): 433?436. LAI Xinmin, HUANG Tian. NURBS based approach of surface reconstruction from scattered data[J]. Journal of Computer Aided Design and Computer Graphics, 1999, 11(5): 433?436.
[13] 李娜. 空間網(wǎng)格結(jié)構(gòu)幾何形態(tài)研究與實現(xiàn)[D]. 浙江大學(xué)建筑工程學(xué)院, 2009: 77?92. LI Na. Research and implementation of geometric configuration of space grid structures[D]. Hangzhou: Zhejiang University. College of Civil Engineering and Architecture, 2009: 77?92.
[14] PIEGL L, TILLER W. The NURBS book[M]. New York, USA: Springer, 1995: 94?95.
[15] 官火梁, 吳強, 席平. RCS計算中NURBS曲面和射線求交的快速計算[J]. 圖學(xué)學(xué)報, 2006, 27(1): 87?91. GUAN Huoliang, WU Qiang, XI Ping. A fast algorithm for intersection calculation of ray and NURBS surface in predicting radar cross section calculation[J]. Journal of Engineering Graphics, 2006, 27(1): 87?91.
[16] 王保慶, 張俐, 李東升. 逆向工程中NURBS曲面與直線交點快速計算[J]. 圖學(xué)學(xué)報, 2010, 31(2): 155?158. WANG Baoqing, ZHANG Li, LI Dongsheng. Rapid calculation of intersection points between NURBS surface and line in Reverse Engineering[J]. Journal of Engineering Graphics, 2010, 31(2): 155?158.
[17] BOOR C D. Cutting corners always works[J]. Computer Aided Geometric Design, 1987, 4(S1/2): 125?131.
Grid generation on complex free-form surface based on surface fitting
CHEN Lijie1, WU Hui2, LI Tierui1, GAO Boqing1
(1. College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China;2.Public Administration College, Zhejiang University of Finance & Economics, Hangzhou 310018, China)
In order to obtain architectural mesh of complex free-form surface, the complex free-form surface was fitted into one surface and then grids were generated using mapping method. First, two methods called NURBS based on ray-surface intersection and NURBS based on plane-surface intersection were proposed to obtain the×point cloud of complex free-form surface. Secondly, a NURBS surface was constructed using the×point cloud as control points, then the boundaries of trimmed NURBS surface were divided by a given number or a given length, and the boundary nodes of the parametric domain were connected to get parametric domain grids. After that, the grids were mapped onto the free-form surface. Lastly, the grids near the boundaries were projected onto original surfaces to reduce boundary fold and shrinkage caused by surface fitting. To reduce the error of surface fitting, the number of point cloud was increased. The results show that compared with the method which respectively generates grids on multiple trimmed surfaces, the proposed methods produce fluent and uniform grids, with good architectural aesthetics.
non-uniform rational B-splines (NURBS); complex free-form surface; surface fitting; mapping method; grid generation
10.11817/j.issn.1672-7207.2018.07.019
TU393.3
A
1672?7207(2018)07?1718?08
2017?07?02;
2017?08?20
國家自然科學(xué)基金資助項目(51378457,51678521);浙江省自然科學(xué)基金資助項目(LY15E080017) (Projects(51378457, 51678521) supported by the National Natural Science Foundation of China; Project(LY15E080017) supported by the Natural Science Foundation of Zhejiang Province)
高博青,博士,教授,從事空間結(jié)構(gòu)研究;E-mail: bqgao@zju.edu.cn
(編輯 伍錦花)