提要

    记得大三那一年有一门课叫做高等有限元，最后的作业就是网格剖分算法的实现，我和同学一起花了些时间做了一个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  #include  #include  #include  #include  #include  #include   using namespace std;  typedef struct { 	double x; 	double y; 	double z; }Point;//定义点类 typedef vector PointArray;//定义点类的vector容器  typedef struct { 	int left; 	int right; 	int count;//边的计数，如果计数为0，则删除此边 }Edge;//定义边类 typedef vector EdgeArray;//定义边类的vector容器  typedef struct { 	int v[3];//三角形的三个顶点 	Edge s[3];//三角形的三条边 	double xc;//三角形外接圆圆心的x坐标 	double yc;//三角形外接圆圆心的y坐标 	double r;//三角形外接圆的半径 }Triangle;//定义三角形类 typedef vector TriangleArray;//定义三角形类的vector容器  typedef vector 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

delaunary .cpp

#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;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)//若要删除三角形的一边的计数为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"<=(fromPoint-1)&&m_Tris[i].v[2]<=(toPoint-1))         {             DelTriangle(i,BoundEdges);             BoundEdges.resize(0);             i--;         }     } }

main.cpp

#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  
  
测试一组数据后，得到结果：

  

  
编译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
另外有需要云服务器可以了解下创新互联scvps.cn，海内外云服务器15元起步，三天无理由+7*72小时售后在线，公司持有idc许可证，提供“云服务器、裸金属服务器、高防服务器、香港服务器、美国服务器、虚拟主机、免备案服务器”等云主机租用服务以及企业上云的综合解决方案，具有“安全稳定、简单易用、服务可用性高、性价比高”等特点与优势，专为企业上云打造定制，能够满足用户丰富、多元化的应用场景需求。
                

                名称栏目：专注网格剖分-TetGen,NETGEN,Steller-创新互联                

                新闻来源：http://hbruida.cn/article/ceeeod.html