記得大三那一年有一門(mén)課叫做高等有限元，最后的作業(yè)就是網(wǎng)格剖分算法的實(shí)現(xiàn)，我和同學(xué)一起花了些時(shí)間做了一個(gè)Qt程序，他寫(xiě)算法，我寫(xiě)界面，最后成績(jī)竟然出奇的拿了90多...
創(chuàng)新互聯(lián)建站服務(wù)項(xiàng)目包括周寧網(wǎng)站建設(shè)、周寧網(wǎng)站制作、周寧網(wǎng)頁(yè)制作以及周寧網(wǎng)絡(luò)營(yíng)銷(xiāo)策劃等。多年來(lái)，我們專(zhuān)注于互聯(lián)網(wǎng)行業(yè)，利用自身積累的技術(shù)優(yōu)勢(shì)、行業(yè)經(jīng)驗(yàn)、深度合作伙伴關(guān)系等，向廣大中小型企業(yè)、政府機(jī)構(gòu)等提供互聯(lián)網(wǎng)行業(yè)的解決方案，周寧網(wǎng)站推廣取得了明顯的社會(huì)效益與經(jīng)濟(jì)效益。目前，我們服務(wù)的客戶(hù)以成都為中心已經(jīng)輻射到周寧省份的部分城市，未來(lái)相信會(huì)繼續(xù)擴(kuò)大服務(wù)區(qū)域并繼續(xù)獲得客戶(hù)的支持與信任！今天要介紹的這款軟件TetGen就是一款網(wǎng)格剖分的軟件，算是力學(xué)計(jì)算中的前處理,他能夠?qū)⑤斎氲娜S模型剖分成一個(gè)個(gè)的單元，如下圖：
最左邊的是原三維模型，中間圖為Delaunay算法生成的四面體網(wǎng)格，最右邊的圖為在tetview中查看剖分的結(jié)果。
官網(wǎng)的手冊(cè)里還有一些關(guān)于剖分算法的說(shuō)明，有興趣的可以去看看。
官網(wǎng)：http://tetgen.berlios.de/
Netgen也是一款網(wǎng)格剖分軟件，為奧地利科學(xué)家Joachim Schoeberl負(fù)責(zé)編寫(xiě)的格網(wǎng)（曲面和實(shí)體）剖分程序。是格網(wǎng)劃分技術(shù)中極為先進(jìn)與完善的，在3D格網(wǎng)劃分領(lǐng)域更是具有極大的優(yōu)勢(shì)。
官網(wǎng)：http://www.hpfem.jku.at/netgen/
Stellar的中文意思是恒星，這是一個(gè)博士寫(xiě)的用于優(yōu)化網(wǎng)格的軟件，可以將生成的單元模型進(jìn)行一些smooth、刪除重復(fù)邊的操作。
環(huán)境: ubuntu 12.04 32bit
【定義】三角剖分[1]：假設(shè)V是二維實(shí)數(shù)域上的有限點(diǎn)集，邊e是由點(diǎn)集中的點(diǎn)作為端點(diǎn)構(gòu)成的封閉線(xiàn)段, E為e的集合。那么該點(diǎn)集V的一個(gè)三角剖分T=(V,E)是一個(gè)平面圖G，該平面圖滿(mǎn)足條件：
1.除了端點(diǎn)，平面圖中的邊不包含點(diǎn)集中的任何點(diǎn)。
2.沒(méi)有相交邊。
3.平面圖中所有的面都是三角面，且所有三角面的合集是散點(diǎn)集V的凸包。
在實(shí)際中運(yùn)用的最多的三角剖分是Delaunay三角剖分，它是一種特殊的三角剖分。
先從Delaunay邊說(shuō)起：
【定義】Delaunay邊：假設(shè)E中的一條邊e（兩個(gè)端點(diǎn)為a,b），e若滿(mǎn)足下列條件，則稱(chēng)之為Delaunay邊：存在一個(gè)圓經(jīng)過(guò)a,b兩點(diǎn)，圓內(nèi)(注意是圓內(nèi)，圓上最多三點(diǎn)共圓)不含點(diǎn)集V中任何其他的點(diǎn)，這一特性又稱(chēng)空?qǐng)A特性。
【定義】Delaunay三角剖分：如果點(diǎn)集V的一個(gè)三角剖分T只包含Delaunay邊，那么該三角剖分稱(chēng)為Delaunay三角剖分。
算法描述
Bowyer-Watson算法
的基本步驟是：
1、構(gòu)造一個(gè)超級(jí)三角形，包含所有散點(diǎn)，放入三角形鏈表。
2、將點(diǎn)集中的散點(diǎn)依次插入，在三角形鏈表中找出外接圓包含插入點(diǎn)的三角形（稱(chēng)為該點(diǎn)的影響三角形），刪除影響三角形的公共邊，將插入點(diǎn)同影響三角形的全部頂點(diǎn)連接起來(lái)，完成一個(gè)點(diǎn)在Delaunay三角形鏈表中的插入。
3、根據(jù)優(yōu)化準(zhǔn)則對(duì)局部新形成的三角形優(yōu)化。將形成的三角形放入Delaunay三角形鏈表。
4、循環(huán)執(zhí)行上述第2步，直到所有散點(diǎn)插入完畢。
算法實(shí)現(xiàn)
代碼是當(dāng)時(shí)的隊(duì)友小六子的,注釋比較詳盡。
delaunay.h
#ifndef DELAUNAY_H_INCLUDED #define DELAUNAY_H_INCLUDED #include#include #include #include #include #include #include using namespace std; typedef struct { double x; double y; double z; }Point;//定義點(diǎn)類(lèi) typedef vector PointArray;//定義點(diǎn)類(lèi)的vector容器 typedef struct { int left; int right; int count;//邊的計(jì)數(shù)，如果計(jì)數(shù)為0，則刪除此邊 }Edge;//定義邊類(lèi) typedef vector EdgeArray;//定義邊類(lèi)的vector容器 typedef struct { int v[3];//三角形的三個(gè)頂點(diǎn) Edge s[3];//三角形的三條邊 double xc;//三角形外接圓圓心的x坐標(biāo) double yc;//三角形外接圓圓心的y坐標(biāo) double r;//三角形外接圓的半徑 }Triangle;//定義三角形類(lèi) typedef vector TriangleArray;//定義三角形類(lèi)的vector容器 typedef vector intArray;//定義int類(lèi)的vector容器 class Delaunay//定義Delaunay類(lèi) { public: Delaunay(Point p1,Point p2,Point p3,Point p4);//Delaunay類(lèi)的構(gòu)造函數(shù)，創(chuàng)建外邊框 ~Delaunay();//Delaunay類(lèi)的析構(gòu)函數(shù) bool AddPoint(double xx,double yy,double zz);//向已有剖分圖形中加點(diǎn)的函數(shù) void Delete_Frame();//刪除外邊框 void Boundary_Recover(int fromPoint,int toPoint);//邊界恢復(fù) void output();//輸出ANSYS命令流文件 private: void Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3);//計(jì)算三角形的外接圓圓心坐標(biāo)和半徑 void MakeTriangle(int n1,int n2,int n3);//生成指定頂點(diǎn)的三角形 bool inCircle(double xx,double yy,Triangle currentTris);//判斷點(diǎn)是否在圓內(nèi) void DelTriangle(int n,EdgeArray &BoundEdges);//刪除指定的三角形 PointArray m_Pts;//m_Pts用于存儲(chǔ)所有點(diǎn) EdgeArray m_Edges;//m_Edges用于存儲(chǔ)所有邊 TriangleArray m_Tris;//m_Tris用于存儲(chǔ)所有三角形 }; void GetPoint(double &xx,double &yy,double &zz,string line);//解析從input文件中讀取的每一行數(shù)據(jù) #endif // DELAUNAY_H_INCLUDED
重慶分公司
重慶分公司
#include "delaunay.h" Delaunay::Delaunay(Point p1,Point p2,Point p3,Point p4) { m_Pts.resize(4); m_Pts[0]=p1; m_Pts[1]=p2; m_Pts[2]=p3; m_Pts[3]=p4;//添加四個(gè)外邊框點(diǎn) m_Edges.resize(4); Edge l1={0,1,-1}; Edge l2={1,2,-1}; Edge l3={0,3,-1}; Edge l4={2,3,-1}; m_Edges[0]=l1; m_Edges[1]=l2; m_Edges[2]=l3; m_Edges[3]=l4;//添加四個(gè)外邊框的邊 MakeTriangle(0,1,2); MakeTriangle(0,2,3);//添加初始的兩個(gè)三角形 } Delaunay::~Delaunay()//清空Delaunay類(lèi)的數(shù)據(jù)成員 { m_Pts.resize(0); m_Edges.resize(0); m_Tris.resize(0); } void Delaunay::MakeTriangle(int n1,int n2,int n3) { double x_centre,y_centre,radius; Cal_Centre(x_centre,y_centre,radius,n1,n2,n3);//獲得頂點(diǎn)為n1,n2,n3的三角形的外接圓圓心坐標(biāo)和半徑 Triangle newTriangle={{n1,n2,n3},{{n1,n2,1},{n2,n3,1},{n1,n3,1}},x_centre,y_centre,radius};//生成指定的三角形 m_Tris.push_back(newTriangle);//向m_Tris中添加新構(gòu)造的三角形 int EdgeSzie=(int)m_Edges.size();//獲得目前的邊數(shù) int flag; for (int i=0;i<3;i++) { flag=1; for(int j=0;jBoundEdges[i].left && PtSize-1currentTris.r) return false; else return true; } void Delaunay::DelTriangle(int n,EdgeArray &BoundEdges) { for (int i=0;i<3;i++) { for (int j=0;j<(int)m_Edges.size();j++) { if (m_Edges[j].left==m_Tris[n].s[i].left&&m_Edges[j].right==m_Tris[n].s[i].right) { if (m_Edges[j].count==2)//若要?jiǎng)h除三角形的一邊的計(jì)數(shù)為2，則將其計(jì)數(shù)減1，并將其壓入BoundEdges容器中 { m_Edges[j].count=1; BoundEdges.push_back(m_Edges[j]); } else if (m_Edges[j].count==-1) BoundEdges.push_back(m_Edges[j]);//如果是外邊框，則直接壓入BoundEdges容器中 else if (m_Edges[j].count==1)//如果刪除三角形的一邊的計(jì)數(shù)為1，則刪除該邊，同時(shí)查看BoundEdges中是否有此邊，若有，則刪除 { for (int k=0;k<(int)BoundEdges.size();k++) { if (BoundEdges[k].left==m_Edges[j].left&&BoundEdges[k].right==m_Edges[j].right) { BoundEdges.erase(BoundEdges.begin()+k); break; } } m_Edges.erase(m_Edges.begin()+j); j--; } break; } } } m_Tris.erase(m_Tris.begin()+n);//刪除該三角形 } void Delaunay::output()//向“output.log"文件中寫(xiě)入ANSYS命令流 { ofstream outfile("output.log"); if (!outfile) { cout<<"Unable to output nodes!"; exit(1); } outfile<<"/PREP7"<=(fromPoint-1)&&m_Tris[i].v[2]<=(toPoint-1)) { DelTriangle(i,BoundEdges); BoundEdges.resize(0); i--; } } } #include "delaunay.h" int main() { ifstream infile("input.txt");//打開(kāi)"input.txt"文件 if (!infile)//判斷文件是否正常打開(kāi) { cout<<"Unable to input nodes!"; exit(1); } string line; PointArray p; double xx,yy,zz; int nodeSize; for (int i=0;i<4;i++)//讀入4外邊框點(diǎn) { getline(infile,line); GetPoint(xx,yy,zz,line); Point tmp={xx,yy,zz}; p.push_back(tmp); } Delaunay MyMesh(p[0],p[1],p[2],p[3]);//實(shí)例化Delaunay類(lèi) getline(infile,line);//讀入節(jié)點(diǎn)數(shù)，用于后面循環(huán) char *cstr; cstr=new char[line.size()+1]; strcpy(cstr,line.c_str()); nodeSize=atoi(cstr); for (int i=0;i
測(cè)試一組數(shù)據(jù)后，得到結(jié)果：
編譯tetgen
下載源碼之后cd進(jìn)目錄，然后執(zhí)行
make
編譯完成之后，目錄下就會(huì)生成一個(gè)名為 tetgen 的可執(zhí)行文件。
運(yùn)行tetview
這個(gè)是用于查看網(wǎng)格模型的工具。 因?yàn)檫@個(gè)東西比較老，所以首先要安裝一些比較老的庫(kù)。
g77
下載好之后解壓，cd進(jìn)目錄運(yùn)行：
sudo ./install.sh
stdc++5
sudo apt-get install libstdc++5
將下載好linux版本的tetivew解壓，再將example解壓到相同的目錄，終端cd進(jìn)目錄，執(zhí)行：
./tetview pmdc.1
一切配置正確的話(huà)，tetview就運(yùn)行了。很簡(jiǎn)單的一個(gè)操作界面，按F1沿著plane剖分，效果就像這樣：
網(wǎng)格剖分實(shí)戰(zhàn)
首先打開(kāi)blender，Add->Mesh->Torus,添加一個(gè)圓環(huán)，然后File->Export->Stanford(.ply)，導(dǎo)出ply文件，待會(huì)用于剖分。
將導(dǎo)出的ply模型放到tetgen的目錄，終端執(zhí)行：
./tetgen -p torus.ply
再將生成的文件拷貝到tetiew的目錄下，執(zhí)行
./tetview torus.1.ele
Netgen
這個(gè)東西編譯起來(lái)還是有點(diǎn)頭疼，還在ubuntu的軟件中心有帶，所以直接在軟件中心搜索下載就可以了。
還是選擇用blender導(dǎo)出模型。這里一定要記住，所有用于網(wǎng)格剖分的模型都要是封閉的體模型，不然就會(huì)出現(xiàn)閃退的情況。
這里選擇一個(gè)植物模型，F(xiàn)ile ->Export->stl。記住勾選左邊的ascii。
打開(kāi)netgen，F(xiàn)ile ->Load Geometry,選擇剛才導(dǎo)出的模型。然后點(diǎn)擊工具欄中的GnerateMesh，稍等片刻，得到結(jié)果。
導(dǎo)出單元
首先選擇導(dǎo)出類(lèi)型：
File -> Export File type ->Elmer Format
然后導(dǎo)出：
File-> Export Mesh
Stellar
從官網(wǎng)下載好源碼之后解壓，終端進(jìn)入目錄，運(yùn)行
make
Stellar就編譯好了。
將之前的用tetgen生成的 model.1.node 和 model.1.ele 文件拷貝至Stellar的文件夾，終端執(zhí)行
Stellar model.1
發(fā)現(xiàn)報(bào)錯(cuò)：
Improving mesh.
***** ALERT Input mesh has non-positive worst quality of -0.85263, dying *****
Starting up smoothing, input tet has quality -0.85263
Stellar: ./src/smoothing.c:1640: nonsmooth: Assertion `worstqual > 0.0' failed.
Aborted (core dumped)
發(fā)郵件為問(wèn)作者，說(shuō)是單元模型三角面沒(méi)有遵循右手法則，用meshconvert.py官網(wǎng)給的腳本轉(zhuǎn)化一下就好。
終端執(zhí)行./meshconvert.py model.1.node model.2.node
執(zhí)行完成之后會(huì)生成新的ele,node文件，這時(shí)再在終端運(yùn)行Stellar,
Stellar model.2
原來(lái)的模型有6000多個(gè)頂點(diǎn)，經(jīng)過(guò)大概10分鐘的優(yōu)化，生成了一個(gè)20000點(diǎn)的模型...T T
原因可能是在平滑處理的過(guò)程中插入了很多點(diǎn)，在優(yōu)化結(jié)果中，還會(huì)生成一個(gè)stats文件，里面描述了整個(gè)優(yōu)化過(guò)程。
如果要控制優(yōu)化的過(guò)程的話(huà)，需要自己寫(xiě)配置文件，修改一下官網(wǎng)給的模板就可以了，比如我不想增加單元格的數(shù)量，則關(guān)閉頂點(diǎn)的插入就可以了。
創(chuàng)建一個(gè) conf 文件
#################################### # Stellar mesh improvement options # #################################### # This file contains all the possible options that can currently be set for # mesh improvement. Lines that begin with '#' are comments and are ignored. # Other lines take the form 'option intval floatval stringval', where option # is the name of the option, and intval floatval and stringval are the possible # fields that can be used to set that option. If an option takes an int, only # a value for int needs to be given. If it's a float, a dummy int should be # inserted before the float. If it's a string, a dummy int and a dummy float # should be inserted before the string. This clumsiness is because I don't like # coding string processing in C, and this is the easiest way. Any unset options # will assume their default values. # verbosity: bigger number means more verbose output of improvement. # default = 1 verbosity 0 # use color in verbose improvement output. default = 0 usecolor 1 # just output the mesh unchanged and quit. default = 0 outputandquit 0 ## quality measure options # qualmeasure: selects which quality measure to use as an objective function # for optimizing the tetrahedra. The quality measures are described in # Chapter 2 of Bryan's dissertation. default = 0 # 0 = minimum sine of the dihedral angles (default) # 1 = square root of radius ratio (circumradius divided by inradius) # 2 = V / l_rms^3 (volume divided by cube of root-mean-squared edge length) # 5 = minimum sine with biased obtuse angles qualmeasure 5 # sinewarpfactor: float. for qualmeasure 5 only; sets the factor by which # obtuse angles are scaled relative to acute angles. Default is 0.75 sinewarpfactor 0 0.75 ## termination options # BOTH goal angles must be reached to terminate improvement # goalanglemin: float. terminates improvement early if minimum angle reaches # this value. default = 90.0 (which precludes early termination) goalanglemin 0 90.0 # goalanglemax: float. terminates improvement early if maximum angle reaches # this value. default = 90.0 goalanglemax 0 90.0 ## smoothing options # nonsmooth: enable optimization-based smoothing. default = 1 nonsmooth 1 # facetsmooth: enable smoothing of facet vertices. default = 1 facetsmooth 1 # segmentsmooth: enable smoothing of segment vertices. default = 1 segmentsmooth 1 # usequadrics: enable use of surface quadric error for smoothing fixed boundary # vertices. WARNING: this will allow the domain shape to change slightly. But # even a little play sometimes yields much better meshes. default = 0 usequadrics 0 # quadricoffset: amount to start quadric error at before subtracting. # See alpha in Section 3.2.5 of Bryan's dissertation. default = 0.8 quadricoffset 0 0.8 # quadricscale: amount to scale quadric error before subtracting from offset. # See beta in Section 3.2.5 of Bryan's dissertation. default = 300.0 quadricscale 0 300.0 ## topological transformation options # flip22: enable 2-2 flips (for boundary tets). default = 1 flip22 1 # multifaceremoval: enable multi-face removal. singlefaceremoval might still # be on. default = 1 multifaceremoval 1 # singlefaceremoval: enable single face removal (2-3 and 2-2 flips). Has # no effect when multifaceremoval is enabled. default = 1 singlefaceremoval 1 # edgeremoval: enable edge removal. default = 1 edgeremoval 1 ## edge contraction options # edgecontraction: enable edge contraction. default = 1 edgecontraction 1 ## vertex insertion options # enableinsert: enable ALL vertex insertion (overrides others). default = 1 enableinsert 0 # insertbody: enable just vertex insertion in body (interior). default = 1 insertbody 0 # insertfacet: enable just insertion on facets. default = 1 insertfacet 0 # insertsegment: enable just insertion on segments. default = 1 insertsegment 0 # insertthreshold: on each insertion pass, try vertex insertion in this fraction of the tetrahedra. default = 0.031 (the worst 3.1%) insertthreshold 0 0.031 ## size control options # (See Chapter 6 of Bryan's dissertation.) # sizing: enable control of element sizes. default = 0 sizing 0 # sizingpass: enable edge length correction before quality improvement. # default = 0 sizingpass 0 # targetedgelength: the target edge length for this mesh. If set to 0.0, the # target edge length is initialized automatically to the initial mean edge # length. default = 0.0 targetedgelength 0 0.0 # longerfactor: factor by which an edge can be longer than the target edge # length before being considered "too long". default = 3.0 longerfactor 0 2.0 # shorterfactor: factor by which an edge can be shorter than the target edge # length before being considered "too short" default = 0.33 shorterfactor 0 0.50 ## anisotropy options # (See Chapter 7 of Bryan's dissertation.) # anisotropic: enable anisotropic meshing. default = 0 anisotropic 0 # tensor: which size/anisotropy tensor to use. default = 0 # 0 = identity # 1 = stretch x # 2 = stretch y # 3 = sink # 4 = swirl # 5 = center # 6 = perimeter # 7 = right # 8 = sine tensor 6 ## quality file output options # These files list, for each tetrahedron, the values of the quality measures # minsineout: enable output of .minsine quality file. default = 1 minsineout 1 # minangout: enable output of .minang quality file. default = 0 minangout 0 # maxangout: enable output of .maxang quality file. default = 0 maxangout 0 # vlrmsout: enable output of .vlrms quality file. default = 0 vlrmsout 0 # nrrout: enable output of the .nrr quality file. default = 0 nrrout 0 ## animation options # animate: activate animation file output (a set of output files after each # pass). default = 0 animate 0 # timeseries: when animate = 1, only output .stats. default = 0 timeseries 0 ## output filename option # fileprefix: filename prefix that distinguishes the output files from the # input files. If none specified, an iteration number is appended to the input # filenames. #fileprefix 0 5 ../output/testprefix
再次運(yùn)行，
./Stellar -s conf model.2
運(yùn)行結(jié)果：
頂點(diǎn)從6000多降到了5000多，用tetiew來(lái)查看：
相關(guān)下載
g77
另外有需要云服務(wù)器可以了解下創(chuàng)新互聯(lián)scvps.cn，海內(nèi)外云服務(wù)器15元起步，三天無(wú)理由+7*72小時(shí)售后在線(xiàn)，公司持有idc許可證，提供“云服務(wù)器、裸金屬服務(wù)器、高防服務(wù)器、香港服務(wù)器、美國(guó)服務(wù)器、虛擬主機(jī)、免備案服務(wù)器”等云主機(jī)租用服務(wù)以及企業(yè)上云的綜合解決方案，具有“安全穩(wěn)定、簡(jiǎn)單易用、服務(wù)可用性高、性?xún)r(jià)比高”等特點(diǎn)與優(yōu)勢(shì)，專(zhuān)為企業(yè)上云打造定制，能夠滿(mǎn)足用戶(hù)豐富、多元化的應(yīng)用場(chǎng)景需求。
網(wǎng)站欄目：專(zhuān)注網(wǎng)格剖分-TetGen,NETGEN,Steller-創(chuàng)新互聯(lián)
標(biāo)題URL：http://weahome.cn/article/ceeeod.html
