專注網格剖分 - TetGen,NETGEN,Steller

提要

       記得大三那一年有一門課叫做高等有限元，最后的作業就是網格剖分算法的實現，我和同學一起花了些時間做了一個Qt程序，他寫算法，我寫界面，最后成績竟然出奇的拿了90多...

        今天要介紹的這款軟件TetGen就是一款網格剖分的軟件，算是力學計算中的前處理,他能夠將輸入的三維模型剖分成一個個的單元，如下圖：

       最左邊的是原三維模型，中間圖為Delaunay算法生成的四面體網格，最右邊的圖為在tetview中查看剖分的結果。

       官網的手冊里還有一些關于剖分算法的說明，有興趣的可以去看看。

       官網：http://tetgen.berlios.de/

Netgen也是一款網格剖分軟件，為奧地利科學家Joachim Schoeberl負責編寫的格網（曲面和實體）剖分程序。是格網劃分技術中極為先進與完善的，在3D格網劃分領域更是具有極大的優勢。

官網：http://www.hpfem.jku.at/netgen/

Stellar的中文意思是恒星，這是一個博士寫的用于優化網格的軟件，可以將生成的單元模型進行一些smooth、刪除重復邊的操作。

       環境: ubuntu 12.04 32bit 

Delaunay算法簡介及實現

【定義】三角剖分[1]：假設V是二維實數域上的有限點集，邊e是由點集中的點作為端點構成的封閉線段, E為e的集合。那么該點集V的一個三角剖分T=(V,E)是一個平面圖G，該平面圖滿足條件：
1.除了端點，平面圖中的邊不包含點集中的任何點。
2.沒有相交邊。
3.平面圖中所有的面都是三角面，且所有三角面的合集是散點集V的凸包。
在實際中運用的最多的三角剖分是Delaunay三角剖分，它是一種特殊的三角剖分。

先從Delaunay邊說起：
【定義】Delaunay邊：假設E中的一條邊e（兩個端點為a,b），e若滿足下列條件，則稱之為Delaunay邊：存在一個圓經過a,b兩點，圓內(注意是圓內，圓上最多三點共圓)不含點集V中任何其他的點，這一特性又稱空圓特性。

【定義】Delaunay三角剖分：如果點集V的一個三角剖分T只包含Delaunay邊，那么該三角剖分稱為Delaunay三角剖分。

算法描述

Bowyer-Watson算法
的基本步驟是：
1、構造一個超級三角形，包含所有散點，放入三角形鏈表。
2、將點集中的散點依次插入，在三角形鏈表中找出外接圓包含插入點的三角形（稱為該點的影響三角形），刪除影響三角形的公共邊，將插入點同影響三角形的全部頂點連接起來，完成一個點在Delaunay三角形鏈表中的插入。
3、根據優化準則對局部新形成的三角形優化。將形成的三角形放入Delaunay三角形鏈表。
4、循環執行上述第2步，直到所有散點插入完畢。

算法實現

代碼是當時的隊友小六子的,注釋比較詳盡。

delaunay.h

#ifndef DELAUNAY_H_INCLUDED #define DELAUNAY_H_INCLUDED #include <cstdlib> #include <iostream> #include <cstring> #include <string> #include <fstream> #include <math.h> #include <vector>  using namespace std;  typedef struct { 	double x; 	double y; 	double z; }Point;//定義點類 typedef vector<Point> PointArray;//定義點類的vector容器  typedef struct { 	int left; 	int right; 	int count;//邊的計數，如果計數為0，則刪除此邊 }Edge;//定義邊類 typedef vector<Edge> EdgeArray;//定義邊類的vector容器  typedef struct { 	int v[3];//三角形的三個頂點 	Edge s[3];//三角形的三條邊 	double xc;//三角形外接圓圓心的x坐標 	double yc;//三角形外接圓圓心的y坐標 	double r;//三角形外接圓的半徑 }Triangle;//定義三角形類 typedef vector<Triangle> TriangleArray;//定義三角形類的vector容器  typedef vector<int> intArray;//定義int類的vector容器  class Delaunay//定義Delaunay類 { public: 	Delaunay(Point p1,Point p2,Point p3,Point p4);//Delaunay類的構造函數，創建外邊框 	~Delaunay();//Delaunay類的析構函數  	bool AddPoint(double xx,double yy,double zz);//向已有剖分圖形中加點的函數 	void Delete_Frame();//刪除外邊框 	void Boundary_Recover(int fromPoint,int toPoint);//邊界恢復 	void output();//輸出ANSYS命令流文件 private:     void Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3);//計算三角形的外接圓圓心坐標和半徑     void MakeTriangle(int n1,int n2,int n3);//生成指定頂點的三角形     bool inCircle(double xx,double yy,Triangle currentTris);//判斷點是否在圓內     void DelTriangle(int n,EdgeArray &BoundEdges);//刪除指定的三角形  	PointArray m_Pts;//m_Pts用于存儲所有點 	EdgeArray m_Edges;//m_Edges用于存儲所有邊 	TriangleArray m_Tris;//m_Tris用于存儲所有三角形 }; void GetPoint(double &xx,double &yy,double &zz,string line);//解析從input文件中讀取的每一行數據 #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;//添加四個外邊框點     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;//添加四個外邊框的邊     MakeTriangle(0,1,2);     MakeTriangle(0,2,3);//添加初始的兩個三角形 }  Delaunay::~Delaunay()//清空Delaunay類的數據成員 {     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);//獲得頂點為n1,n2,n3的三角形的外接圓圓心坐標和半徑     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中添加新構造的三角形     int EdgeSzie=(int)m_Edges.size();//獲得目前的邊數     int flag;     for (int i=0;i<3;i++)     {         flag=1;         for(int j=0;j<EdgeSzie;j++)//通過循環判斷新構造的三角形的各邊是否已經存在于m_Edges中，如果存在則只增加該邊的計數，否則添加新邊         {             if (newTriangle.s[i].left==m_Edges[j].left&&newTriangle.s[i].right==m_Edges[j].right&&m_Edges[j].count!=-1) {flag=0;m_Edges[j].count+=1;break;}             else if(newTriangle.s[i].left==m_Edges[j].left&&newTriangle.s[i].right==m_Edges[j].right&&m_Edges[j].count==-1) {flag=0;break;}         }         if (flag==1) m_Edges.push_back(newTriangle.s[i]);     } }  void Delaunay::Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3) {     double x1,x2,x3,y1,y2,y3;     x1=m_Pts[n1].x;     y1=m_Pts[n1].y;     x2=m_Pts[n2].x;     y2=m_Pts[n2].y;     x3=m_Pts[n3].x;     y3=m_Pts[n3].y;     x_centre=((y2-y1)*(y3*y3-y1*y1+x3*x3-x1*x1)-(y3-y1)*(y2*y2-y1*y1+x2*x2-x1*x1))/(2*(x3-x1)*(y2-y1)-2*((x2-x1)*(y3-y1)));//計算外接圓圓心的x坐標     y_centre=((x2-x1)*(x3*x3-x1*x1+y3*y3-y1*y1)-(x3-x1)*(x2*x2-x1*x1+y2*y2-y1*y1))/(2*(y3-y1)*(x2-x1)-2*((y2-y1)*(x3-x1)));//計算外接圓圓心的y坐標     radius= sqrt((x1 - x_centre)*(x1 - x_centre) + (y1 - y_centre)*(y1 - y_centre));//計算外接圓的半徑 }  bool Delaunay::AddPoint(double xx,double yy,double zz) {     EdgeArray BoundEdges;//BoundEdges用于存儲在刪除三角形后留下的邊框，用于構造新的三角形     Point newPoint={xx,yy,zz};     m_Pts.push_back(newPoint);//向m_Pts中添加新點     intArray badTriangle;//badTriangle用于存儲不符合空圓規則的三角形的索引號     int TriSize=(int)m_Tris.size();//獲得目前的三角形數     for (int i=0;i<TriSize;i++)//通過循環找到所有不符合空圓規則的三角形，并將其索引號存在badTriangle中     {         if (inCircle(xx,yy,m_Tris[i])==true) badTriangle.push_back(i);     }     for (int i=0;i<(int)badTriangle.size();i++)//通過循環刪除所有不符合空圓規則的三角形，同時保留邊框     {         DelTriangle(badTriangle[i],BoundEdges);         for (int j=i+1;j<(int)badTriangle.size();j++) badTriangle[j]-=1;     }     int PtSize=(int)m_Pts.size();//獲得目前的點數     for (int i=0;i<(int)BoundEdges.size();i++)//生成新的三角形     {         if (PtSize-1<BoundEdges[i].left) MakeTriangle(PtSize-1,BoundEdges[i].left,BoundEdges[i].right);         else if (PtSize-1>BoundEdges[i].left && PtSize-1<BoundEdges[i].right) MakeTriangle(BoundEdges[i].left,PtSize-1,BoundEdges[i].right);         else MakeTriangle(BoundEdges[i].left,BoundEdges[i].right,PtSize-1);     }     return true; }  bool Delaunay::inCircle(double xx,double yy,Triangle currentTris)//判斷點是否在三角形的外接圓內 {     double dis=sqrt((currentTris.xc-xx)*(currentTris.xc-xx)+(currentTris.yc-yy)*(currentTris.yc-yy));     if (dis>currentTris.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)//若要刪除三角形的一邊的計數為2，則將其計數減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)//如果刪除三角形的一邊的計數為1，則刪除該邊，同時查看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"文件中寫入ANSYS命令流 {     ofstream outfile("output.log");     if (!outfile)     {         cout<<"Unable to output nodes!";         exit(1);     }     outfile<<"/PREP7"<<endl;     for (int i=0;i<(int)m_Pts.size();i++)     {         outfile<<"K,"<<i+1<<","<<m_Pts[i].x<<","<<m_Pts[i].y<<","<<m_Pts[i].z<<endl;     }     for (int i=0;i<(int)m_Edges.size();i++)     {         outfile<<"L,"<<m_Edges[i].left+1<<","<<m_Edges[i].right+1<<endl;     }     outfile.close(); } void GetPoint(double &xx,double &yy,double &zz,string line)//從字符串line中解析出點的x,y,z坐標 {     int flag=0;     string tmp="";     char *cstr;     for (int i=(int)line.find(',')+1;i<(int)line.size();i++)     {         if (line[i]==',')         {             cstr=new char[tmp.size()+1];             strcpy(cstr,tmp.c_str());             if (flag==0) {xx=atof(cstr);tmp.resize(0);flag++;}             else if (flag==1) {yy=atof(cstr);tmp.resize(0);flag++;}             continue;         }         tmp=tmp+line[i];     }     if (fabs(xx)<1.0e-6) xx=0.0;     if (fabs(yy)<1.0e-6) yy=0.0;     if (fabs(zz)<1.0e-6) zz=0.0; } void Delaunay::Delete_Frame()//刪除外邊框 {     EdgeArray BoundEdges;     for (int i=0;i<4;i++) m_Pts.erase(m_Pts.begin());     for (int i=0;i<(int)m_Tris.size();i++)     {         if (m_Tris[i].v[0]==0||m_Tris[i].v[0]==1||m_Tris[i].v[0]==2||m_Tris[i].v[0]==3)         {             DelTriangle(i,BoundEdges);             BoundEdges.resize(0);             i--;         }         else         {             for (int j=0;j<3;j++)             {                 m_Tris[i].v[j]-=4;                 m_Tris[i].s[j].left-=4;                 m_Tris[i].s[j].right-=4;             }         }     }     for (int i=0;i<4;i++) m_Edges.erase(m_Edges.begin());     for (int i=0;i<(int)m_Edges.size();i++)     {         m_Edges[i].left-=4;         m_Edges[i].right-=4;     } }  void Delaunay::Boundary_Recover(int fromPoint,int toPoint)//恢復由指定點組成的邊界 {     EdgeArray BoundEdges;     for (int i=0;i<(int)m_Tris.size();i++)     {         if (m_Tris[i].v[0]>=(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");//打開"input.txt"文件     if (!infile)//判斷文件是否正常打開     {         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外邊框點     {         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]);//實例化Delaunay類     getline(infile,line);//讀入節點數，用于后面循環     char *cstr;     cstr=new char[line.size()+1];     strcpy(cstr,line.c_str());     nodeSize=atoi(cstr);     for (int i=0;i<nodeSize;i++)//讀入每個節點的坐標     {         getline(infile,line);         GetPoint(xx,yy,zz,line);         MyMesh.AddPoint(xx,yy,zz);     }     infile.close();     MyMesh.Delete_Frame();//刪除外邊框     MyMesh.Boundary_Recover(203,466);     MyMesh.Boundary_Recover(467,487);     MyMesh.Boundary_Recover(488,511);     MyMesh.Boundary_Recover(512,537);//以上都是恢復指定邊界     MyMesh.output();//將相應ANSYS命令流輸出     return 0; }

測試一組數據后，得到結果：

編譯tetgen

下載源碼之后cd進目錄，然后執行

make

編譯完成之后，目錄下就會生成一個名為 tetgen 的可執行文件。

運行tetview

       這個是用于查看網格模型的工具。 因為這個東西比較老，所以首先要安裝一些比較老的庫。

       g77

下載好之后解壓，cd進目錄運行：

sudo ./install.sh

       stdc++5

 sudo apt-get install libstdc++5 

 將下載好linux版本的tetivew解壓，再將example解壓到相同的目錄，終端cd進目錄，執行：

./tetview pmdc.1

一切配置正確的話，tetview就運行了。很簡單的一個操作界面，按F1沿著plane剖分，效果就像這樣：

網格剖分實戰

首先打開blender，Add->Mesh->Torus,添加一個圓環，然后File->Export->Stanford(.ply)，導出ply文件，待會用于剖分。

將導出的ply模型放到tetgen的目錄，終端執行：

./tetgen -p torus.ply 

再將生成的文件拷貝到tetiew的目錄下，執行

./tetview torus.1.ele

Netgen

這個東西編譯起來還是有點頭疼，還在ubuntu的軟件中心有帶，所以直接在軟件中心搜索下載就可以了。

還是選擇用blender導出模型。這里一定要記住，所有用于網格剖分的模型都要是封閉的體模型，不然就會出現閃退的情況。

這里選擇一個植物模型，File ->Export->stl。記住勾選左邊的ascii。

打開netgen，File ->Load Geometry,選擇剛才導出的模型。然后點擊工具欄中的GnerateMesh，稍等片刻，得到結果。

導出單元

首先選擇導出類型：
File -> Export File type ->Elmer Format

然后導出：

File-> Export Mesh

Stellar

從官網下載好源碼之后解壓，終端進入目錄，運行

make

Stellar就編譯好了。

將之前的用tetgen生成的 model.1.node 和 model.1.ele 文件拷貝至Stellar的文件夾，終端執行

Stellar model.1

發現報錯：

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)

發郵件為問作者，說是單元模型三角面沒有遵循右手法則，用meshconvert.py官網給的腳本轉化一下就好。

終端執行./meshconvert.py model.1.node model.2.node

執行完成之后會生成新的ele,node文件，這時再在終端運行Stellar,

Stellar model.2

原來的模型有6000多個頂點，經過大概10分鐘的優化，生成了一個20000點的模型...T T

原因可能是在平滑處理的過程中插入了很多點，在優化結果中，還會生成一個stats文件，里面描述了整個優化過程。

如果要控制優化的過程的話，需要自己寫配置文件，修改一下官網給的模板就可以了，比如我不想增加單元格的數量，則關閉頂點的插入就可以了。

創建一個 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

再次運行，

./Stellar -s conf model.2

運行結果：

頂點從6000多降到了5000多，用tetiew來查看：

相關下載

g77