网格剖分 Delaunay

http://blog.csdn.net/chinamming/article/details/16874371


今天要介绍的这款软件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

  1. #ifndef DELAUNAY_H_INCLUDED   
  2. #define DELAUNAY_H_INCLUDED   
  3. #include <cstdlib>   
  4. #include <iostream>   
  5. #include <cstring>   
  6. #include <string>   
  7. #include <fstream>   
  8. #include <math.h>   
  9. #include <vector>   
  10.   
  11. using namespace std;  
  12.   
  13. typedef struct  
  14. {  
  15.     double x;  
  16.     double y;  
  17.     double z;  
  18. }Point;//定义点类   
  19. typedef vector<Point> PointArray;//定义点类的vector容器  
  20.   
  21. typedef struct  
  22. {  
  23.     int left;  
  24.     int right;  
  25.     int count;//边的计数,如果计数为0,则删除此边  
  26. }Edge;//定义边类   
  27. typedef vector<Edge> EdgeArray;//定义边类的vector容器  
  28.   
  29. typedef struct  
  30. {  
  31.     int v[3];//三角形的三个顶点  
  32.     Edge s[3];//三角形的三条边   
  33.     double xc;//三角形外接圆圆心的x坐标  
  34.     double yc;//三角形外接圆圆心的y坐标  
  35.     double r;//三角形外接圆的半径  
  36. }Triangle;//定义三角形类   
  37. typedef vector<Triangle> TriangleArray;//定义三角形类的vector容器  
  38.   
  39. typedef vector<int> intArray;//定义int类的vector容器  
  40.   
  41. class Delaunay//定义Delaunay类  
  42. {  
  43. public:  
  44.     Delaunay(Point p1,Point p2,Point p3,Point p4);//Delaunay类的构造函数,创建外边框  
  45.     ~Delaunay();//Delaunay类的析构函数  
  46.   
  47.     bool AddPoint(double xx,double yy,double zz);//向已有剖分图形中加点的函数  
  48.     void Delete_Frame();//删除外边框  
  49.     void Boundary_Recover(int fromPoint,int toPoint);//边界恢复  
  50.     void output();//输出ANSYS命令流文件  
  51. private:  
  52.     void Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3);//计算三角形的外接圆圆心坐标和半径  
  53.     void MakeTriangle(int n1,int n2,int n3);//生成指定顶点的三角形  
  54.     bool inCircle(double xx,double yy,Triangle currentTris);//判断点是否在圆内  
  55.     void DelTriangle(int n,EdgeArray &BoundEdges);//删除指定的三角形  
  56.   
  57.     PointArray m_Pts;//m_Pts用于存储所有点  
  58.     EdgeArray m_Edges;//m_Edges用于存储所有边   
  59.     TriangleArray m_Tris;//m_Tris用于存储所有三角形  
  60. };  
  61. void GetPoint(double &xx,double &yy,double &zz,string line);//解析从input文件中读取的每一行数据  
  62. #endif // DELAUNAY_H_INCLUDED  
[cpp] view plain copy
 在CODE上查看代码片派生到我的代码片
  1. #ifndef DELAUNAY_H_INCLUDED  
  2. #define DELAUNAY_H_INCLUDED  
  3. #include <cstdlib>  
  4. #include <iostream>  
  5. #include <cstring>  
  6. #include <string>  
  7. #include <fstream>  
  8. #include <math.h>  
  9. #include <vector>  
  10.   
  11. using namespace std;  
  12.   
  13. typedef struct  
  14. {  
  15.     double x;  
  16.     double y;  
  17.     double z;  
  18. }Point;//定义点类  
  19. typedef vector<Point> PointArray;//定义点类的vector容器  
  20.   
  21. typedef struct  
  22. {  
  23.     int left;  
  24.     int right;  
  25.     int count;//边的计数,如果计数为0,则删除此边  
  26. }Edge;//定义边类  
  27. typedef vector<Edge> EdgeArray;//定义边类的vector容器  
  28.   
  29. typedef struct  
  30. {  
  31.     int v[3];//三角形的三个顶点  
  32.     Edge s[3];//三角形的三条边  
  33.     double xc;//三角形外接圆圆心的x坐标  
  34.     double yc;//三角形外接圆圆心的y坐标  
  35.     double r;//三角形外接圆的半径  
  36. }Triangle;//定义三角形类  
  37. typedef vector<Triangle> TriangleArray;//定义三角形类的vector容器  
  38.   
  39. typedef vector<int> intArray;//定义int类的vector容器  
  40.   
  41. class Delaunay//定义Delaunay类  
  42. {  
  43. public:  
  44.     Delaunay(Point p1,Point p2,Point p3,Point p4);//Delaunay类的构造函数,创建外边框  
  45.     ~Delaunay();//Delaunay类的析构函数  
  46.   
  47.     bool AddPoint(double xx,double yy,double zz);//向已有剖分图形中加点的函数  
  48.     void Delete_Frame();//删除外边框  
  49.     void Boundary_Recover(int fromPoint,int toPoint);//边界恢复  
  50.     void output();//输出ANSYS命令流文件  
  51. private:  
  52.     void Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3);//计算三角形的外接圆圆心坐标和半径  
  53.     void MakeTriangle(int n1,int n2,int n3);//生成指定顶点的三角形  
  54.     bool inCircle(double xx,double yy,Triangle currentTris);//判断点是否在圆内  
  55.     void DelTriangle(int n,EdgeArray &BoundEdges);//删除指定的三角形  
  56.   
  57.     PointArray m_Pts;//m_Pts用于存储所有点  
  58.     EdgeArray m_Edges;//m_Edges用于存储所有边  
  59.     TriangleArray m_Tris;//m_Tris用于存储所有三角形  
  60. };  
  61. void GetPoint(double &xx,double &yy,double &zz,string line);//解析从input文件中读取的每一行数据  
  62. #endif // DELAUNAY_H_INCLUDED  

delaunary .cpp

  1. #include "delaunay.h"   
  2.   
  3. Delaunay::Delaunay(Point p1,Point p2,Point p3,Point p4)  
  4. {  
  5.     m_Pts.resize(4);  
  6.     m_Pts[0]=p1;  
  7.     m_Pts[1]=p2;  
  8.     m_Pts[2]=p3;  
  9.     m_Pts[3]=p4;//添加四个外边框点   
  10.     m_Edges.resize(4);  
  11.     Edge l1={0,1,-1};  
  12.     Edge l2={1,2,-1};  
  13.     Edge l3={0,3,-1};  
  14.     Edge l4={2,3,-1};  
  15.     m_Edges[0]=l1;  
  16.     m_Edges[1]=l2;  
  17.     m_Edges[2]=l3;  
  18.     m_Edges[3]=l4;//添加四个外边框的边   
  19.     MakeTriangle(0,1,2);  
  20.     MakeTriangle(0,2,3);//添加初始的两个三角形   
  21. }  
  22.   
  23. Delaunay::~Delaunay()//清空Delaunay类的数据成员  
  24. {  
  25.     m_Pts.resize(0);  
  26.     m_Edges.resize(0);  
  27.     m_Tris.resize(0);  
  28. }  
  29.   
  30. void Delaunay::MakeTriangle(int n1,int n2,int n3)  
  31. {  
  32.     double x_centre,y_centre,radius;  
  33.     Cal_Centre(x_centre,y_centre,radius,n1,n2,n3);//获得顶点为n1,n2,n3的三角形的外接圆圆心坐标和半径  
  34.     Triangle newTriangle={{n1,n2,n3},{{n1,n2,1},{n2,n3,1},{n1,n3,1}},x_centre,y_centre,radius};//生成指定的三角形  
  35.     m_Tris.push_back(newTriangle);//向m_Tris中添加新构造的三角形  
  36.     int EdgeSzie=(int)m_Edges.size();//获得目前的边数  
  37.     int flag;  
  38.     for (int i=0;i<3;i++)  
  39.     {  
  40.         flag=1;  
  41.         for(int j=0;j<EdgeSzie;j++)//通过循环判断新构造的三角形的各边是否已经存在于m_Edges中,如果存在则只增加该边的计数,否则添加新边  
  42.         {  
  43.             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;}  
  44.             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;}  
  45.         }  
  46.         if (flag==1) m_Edges.push_back(newTriangle.s[i]);  
  47.     }  
  48. }  
  49.   
  50. void Delaunay::Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3)  
  51. {  
  52.     double x1,x2,x3,y1,y2,y3;  
  53.     x1=m_Pts[n1].x;  
  54.     y1=m_Pts[n1].y;  
  55.     x2=m_Pts[n2].x;  
  56.     y2=m_Pts[n2].y;  
  57.     x3=m_Pts[n3].x;  
  58.     y3=m_Pts[n3].y;  
  59.     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坐标  
  60.     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坐标  
  61.     radius= sqrt((x1 - x_centre)*(x1 - x_centre) + (y1 - y_centre)*(y1 - y_centre));//计算外接圆的半径  
  62. }  
  63.   
  64. bool Delaunay::AddPoint(double xx,double yy,double zz)  
  65. {  
  66.     EdgeArray BoundEdges;//BoundEdges用于存储在删除三角形后留下的边框,用于构造新的三角形  
  67.     Point newPoint={xx,yy,zz};  
  68.     m_Pts.push_back(newPoint);//向m_Pts中添加新点  
  69.     intArray badTriangle;//badTriangle用于存储不符合空圆规则的三角形的索引号  
  70.     int TriSize=(int)m_Tris.size();//获得目前的三角形数  
  71.     for (int i=0;i<TriSize;i++)//通过循环找到所有不符合空圆规则的三角形,并将其索引号存在badTriangle中  
  72.     {  
  73.         if (inCircle(xx,yy,m_Tris[i])==true) badTriangle.push_back(i);  
  74.     }  
  75.     for (int i=0;i<(int)badTriangle.size();i++)//通过循环删除所有不符合空圆规则的三角形,同时保留边框  
  76.     {  
  77.         DelTriangle(badTriangle[i],BoundEdges);  
  78.         for (int j=i+1;j<(int)badTriangle.size();j++) badTriangle[j]-=1;  
  79.     }  
  80.     int PtSize=(int)m_Pts.size();//获得目前的点数  
  81.     for (int i=0;i<(int)BoundEdges.size();i++)//生成新的三角形  
  82.     {  
  83.         if (PtSize-1<BoundEdges[i].left) MakeTriangle(PtSize-1,BoundEdges[i].left,BoundEdges[i].right);  
  84.         else if (PtSize-1>BoundEdges[i].left && PtSize-1<BoundEdges[i].right) MakeTriangle(BoundEdges[i].left,PtSize-1,BoundEdges[i].right);  
  85.         else MakeTriangle(BoundEdges[i].left,BoundEdges[i].right,PtSize-1);  
  86.     }  
  87.     return true;  
  88. }  
  89.   
  90. bool Delaunay::inCircle(double xx,double yy,Triangle currentTris)//判断点是否在三角形的外接圆内  
  91. {  
  92.     double dis=sqrt((currentTris.xc-xx)*(currentTris.xc-xx)+(currentTris.yc-yy)*(currentTris.yc-yy));  
  93.     if (dis>currentTris.r) return false;  
  94.     else return true;  
  95. }  
  96.   
  97. void Delaunay::DelTriangle(int n,EdgeArray &BoundEdges)  
  98. {  
  99.     for (int i=0;i<3;i++)  
  100.     {  
  101.         for (int j=0;j<(int)m_Edges.size();j++)  
  102.         {  
  103.             if (m_Edges[j].left==m_Tris[n].s[i].left&&m_Edges[j].right==m_Tris[n].s[i].right)  
  104.             {  
  105.                 if (m_Edges[j].count==2)//若要删除三角形的一边的计数为2,则将其计数减1,并将其压入BoundEdges容器中  
  106.                 {  
  107.                     m_Edges[j].count=1;  
  108.                     BoundEdges.push_back(m_Edges[j]);  
  109.                 }  
  110.                 else if (m_Edges[j].count==-1) BoundEdges.push_back(m_Edges[j]);//如果是外边框,则直接压入BoundEdges容器中  
  111.                 else if (m_Edges[j].count==1)//如果删除三角形的一边的计数为1,则删除该边,同时查看BoundEdges中是否有此边,若有,则删除  
  112.                 {  
  113.                     for (int k=0;k<(int)BoundEdges.size();k++)  
  114.                     {  
  115.                         if (BoundEdges[k].left==m_Edges[j].left&&BoundEdges[k].right==m_Edges[j].right)  
  116.                         {  
  117.                             BoundEdges.erase(BoundEdges.begin()+k);  
  118.                             break;  
  119.                         }  
  120.                     }  
  121.                     m_Edges.erase(m_Edges.begin()+j);  
  122.                     j--;  
  123.                 }  
  124.                 break;  
  125.             }  
  126.         }  
  127.     }  
  128.     m_Tris.erase(m_Tris.begin()+n);//删除该三角形  
  129. }  
  130.   
  131. void Delaunay::output()//向“output.log"文件中写入ANSYS命令流  
  132. {  
  133.     ofstream outfile("output.log");  
  134.     if (!outfile)  
  135.     {  
  136.         cout<<"Unable to output nodes!";  
  137.         exit(1);  
  138.     }  
  139.     outfile<<"/PREP7"<<endl;  
  140.     for (int i=0;i<(int)m_Pts.size();i++)  
  141.     {  
  142.         outfile<<"K,"<<i+1<<","<<m_Pts[i].x<<","<<m_Pts[i].y<<","<<m_Pts[i].z<<endl;  
  143.     }  
  144.     for (int i=0;i<(int)m_Edges.size();i++)  
  145.     {  
  146.         outfile<<"L,"<<m_Edges[i].left+1<<","<<m_Edges[i].right+1<<endl;  
  147.     }  
  148.     outfile.close();  
  149. }  
  150. void GetPoint(double &xx,double &yy,double &zz,string line)//从字符串line中解析出点的x,y,z坐标  
  151. {  
  152.     int flag=0;  
  153.     string tmp="";  
  154.     char *cstr;  
  155.     for (int i=(int)line.find(',')+1;i<(int)line.size();i++)  
  156.     {  
  157.         if (line[i]==',')  
  158.         {  
  159.             cstr=new char[tmp.size()+1];  
  160.             strcpy(cstr,tmp.c_str());  
  161.             if (flag==0) {xx=atof(cstr);tmp.resize(0);flag++;}  
  162.             else if (flag==1) {yy=atof(cstr);tmp.resize(0);flag++;}  
  163.             continue;  
  164.         }  
  165.         tmp=tmp+line[i];  
  166.     }  
  167.     if (fabs(xx)<1.0e-6) xx=0.0;  
  168.     if (fabs(yy)<1.0e-6) yy=0.0;  
  169.     if (fabs(zz)<1.0e-6) zz=0.0;  
  170. }  
  171. void Delaunay::Delete_Frame()//删除外边框  
  172. {  
  173.     EdgeArray BoundEdges;  
  174.     for (int i=0;i<4;i++) m_Pts.erase(m_Pts.begin());  
  175.     for (int i=0;i<(int)m_Tris.size();i++)  
  176.     {  
  177.         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)  
  178.         {  
  179.             DelTriangle(i,BoundEdges);  
  180.             BoundEdges.resize(0);  
  181.             i--;  
  182.         }  
  183.         else  
  184.         {  
  185.             for (int j=0;j<3;j++)  
  186.             {  
  187.                 m_Tris[i].v[j]-=4;  
  188.                 m_Tris[i].s[j].left-=4;  
  189.                 m_Tris[i].s[j].right-=4;  
  190.             }  
  191.         }  
  192.     }  
  193.     for (int i=0;i<4;i++) m_Edges.erase(m_Edges.begin());  
  194.     for (int i=0;i<(int)m_Edges.size();i++)  
  195.     {  
  196.         m_Edges[i].left-=4;  
  197.         m_Edges[i].right-=4;  
  198.     }  
  199. }  
  200.   
  201. void Delaunay::Boundary_Recover(int fromPoint,int toPoint)//恢复由指定点组成的边界  
  202. {  
  203.     EdgeArray BoundEdges;  
  204.     for (int i=0;i<(int)m_Tris.size();i++)  
  205.     {  
  206.         if (m_Tris[i].v[0]>=(fromPoint-1)&&m_Tris[i].v[2]<=(toPoint-1))  
  207.         {  
  208.             DelTriangle(i,BoundEdges);  
  209.             BoundEdges.resize(0);  
  210.             i--;  
  211.         }  
  212.     }  
  213. }  
[cpp] view plain copy
 在CODE上查看代码片派生到我的代码片
  1. #include "delaunay.h"  
  2.   
  3. Delaunay::Delaunay(Point p1,Point p2,Point p3,Point p4)  
  4. {  
  5.     m_Pts.resize(4);  
  6.     m_Pts[0]=p1;  
  7.     m_Pts[1]=p2;  
  8.     m_Pts[2]=p3;  
  9.     m_Pts[3]=p4;//添加四个外边框点  
  10.     m_Edges.resize(4);  
  11.     Edge l1={0,1,-1};  
  12.     Edge l2={1,2,-1};  
  13.     Edge l3={0,3,-1};  
  14.     Edge l4={2,3,-1};  
  15.     m_Edges[0]=l1;  
  16.     m_Edges[1]=l2;  
  17.     m_Edges[2]=l3;  
  18.     m_Edges[3]=l4;//添加四个外边框的边  
  19.     MakeTriangle(0,1,2);  
  20.     MakeTriangle(0,2,3);//添加初始的两个三角形  
  21. }  
  22.   
  23. Delaunay::~Delaunay()//清空Delaunay类的数据成员  
  24. {  
  25.     m_Pts.resize(0);  
  26.     m_Edges.resize(0);  
  27.     m_Tris.resize(0);  
  28. }  
  29.   
  30. void Delaunay::MakeTriangle(int n1,int n2,int n3)  
  31. {  
  32.     double x_centre,y_centre,radius;  
  33.     Cal_Centre(x_centre,y_centre,radius,n1,n2,n3);//获得顶点为n1,n2,n3的三角形的外接圆圆心坐标和半径  
  34.     Triangle newTriangle={{n1,n2,n3},{{n1,n2,1},{n2,n3,1},{n1,n3,1}},x_centre,y_centre,radius};//生成指定的三角形  
  35.     m_Tris.push_back(newTriangle);//向m_Tris中添加新构造的三角形  
  36.     int EdgeSzie=(int)m_Edges.size();//获得目前的边数  
  37.     int flag;  
  38.     for (int i=0;i<3;i++)  
  39.     {  
  40.         flag=1;  
  41.         for(int j=0;j<EdgeSzie;j++)//通过循环判断新构造的三角形的各边是否已经存在于m_Edges中,如果存在则只增加该边的计数,否则添加新边  
  42.         {  
  43.             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;}  
  44.             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;}  
  45.         }  
  46.         if (flag==1) m_Edges.push_back(newTriangle.s[i]);  
  47.     }  
  48. }  
  49.   
  50. void Delaunay::Cal_Centre(double &x_centre,double &y_centre,double &radius,int n1,int n2,int n3)  
  51. {  
  52.     double x1,x2,x3,y1,y2,y3;  
  53.     x1=m_Pts[n1].x;  
  54.     y1=m_Pts[n1].y;  
  55.     x2=m_Pts[n2].x;  
  56.     y2=m_Pts[n2].y;  
  57.     x3=m_Pts[n3].x;  
  58.     y3=m_Pts[n3].y;  
  59.     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坐标  
  60.     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坐标  
  61.     radius= sqrt((x1 - x_centre)*(x1 - x_centre) + (y1 - y_centre)*(y1 - y_centre));//计算外接圆的半径  
  62. }  
  63.   
  64. bool Delaunay::AddPoint(double xx,double yy,double zz)  
  65. {  
  66.     EdgeArray BoundEdges;//BoundEdges用于存储在删除三角形后留下的边框,用于构造新的三角形  
  67.     Point newPoint={xx,yy,zz};  
  68.     m_Pts.push_back(newPoint);//向m_Pts中添加新点  
  69.     intArray badTriangle;//badTriangle用于存储不符合空圆规则的三角形的索引号  
  70.     int TriSize=(int)m_Tris.size();//获得目前的三角形数  
  71.     for (int i=0;i<TriSize;i++)//通过循环找到所有不符合空圆规则的三角形,并将其索引号存在badTriangle中  
  72.     {  
  73.         if (inCircle(xx,yy,m_Tris[i])==true) badTriangle.push_back(i);  
  74.     }  
  75.     for (int i=0;i<(int)badTriangle.size();i++)//通过循环删除所有不符合空圆规则的三角形,同时保留边框  
  76.     {  
  77.         DelTriangle(badTriangle[i],BoundEdges);  
  78.         for (int j=i+1;j<(int)badTriangle.size();j++) badTriangle[j]-=1;  
  79.     }  
  80.     int PtSize=(int)m_Pts.size();//获得目前的点数  
  81.     for (int i=0;i<(int)BoundEdges.size();i++)//生成新的三角形  
  82.     {  
  83.         if (PtSize-1<BoundEdges[i].left) MakeTriangle(PtSize-1,BoundEdges[i].left,BoundEdges[i].right);  
  84.         else if (PtSize-1>BoundEdges[i].left && PtSize-1<BoundEdges[i].right) MakeTriangle(BoundEdges[i].left,PtSize-1,BoundEdges[i].right);  
  85.         else MakeTriangle(BoundEdges[i].left,BoundEdges[i].right,PtSize-1);  
  86.     }  
  87.     return true;  
  88. }  
  89.   
  90. bool Delaunay::inCircle(double xx,double yy,Triangle currentTris)//判断点是否在三角形的外接圆内  
  91. {  
  92.     double dis=sqrt((currentTris.xc-xx)*(currentTris.xc-xx)+(currentTris.yc-yy)*(currentTris.yc-yy));  
  93.     if (dis>currentTris.r) return false;  
  94.     else return true;  
  95. }  
  96.   
  97. void Delaunay::DelTriangle(int n,EdgeArray &BoundEdges)  
  98. {  
  99.     for (int i=0;i<3;i++)  
  100.     {  
  101.         for (int j=0;j<(int)m_Edges.size();j++)  
  102.         {  
  103.             if (m_Edges[j].left==m_Tris[n].s[i].left&&m_Edges[j].right==m_Tris[n].s[i].right)  
  104.             {  
  105.                 if (m_Edges[j].count==2)//若要删除三角形的一边的计数为2,则将其计数减1,并将其压入BoundEdges容器中  
  106.                 {  
  107.                     m_Edges[j].count=1;  
  108.                     BoundEdges.push_back(m_Edges[j]);  
  109.                 }  
  110.                 else if (m_Edges[j].count==-1) BoundEdges.push_back(m_Edges[j]);//如果是外边框,则直接压入BoundEdges容器中  
  111.                 else if (m_Edges[j].count==1)//如果删除三角形的一边的计数为1,则删除该边,同时查看BoundEdges中是否有此边,若有,则删除  
  112.                 {  
  113.                     for (int k=0;k<(int)BoundEdges.size();k++)  
  114.                     {  
  115.                         if (BoundEdges[k].left==m_Edges[j].left&&BoundEdges[k].right==m_Edges[j].right)  
  116.                         {  
  117.                             BoundEdges.erase(BoundEdges.begin()+k);  
  118.                             break;  
  119.                         }  
  120.                     }  
  121.                     m_Edges.erase(m_Edges.begin()+j);  
  122.                     j--;  
  123.                 }  
  124.                 break;  
  125.             }  
  126.         }  
  127.     }  
  128.     m_Tris.erase(m_Tris.begin()+n);//删除该三角形  
  129. }  
  130.   
  131. void Delaunay::output()//向“output.log"文件中写入ANSYS命令流  
  132. {  
  133.     ofstream outfile("output.log");  
  134.     if (!outfile)  
  135.     {  
  136.         cout<<"Unable to output nodes!";  
  137.         exit(1);  
  138.     }  
  139.     outfile<<"/PREP7"<<endl;  
  140.     for (int i=0;i<(int)m_Pts.size();i++)  
  141.     {  
  142.         outfile<<"K,"<<i+1<<","<<m_Pts[i].x<<","<<m_Pts[i].y<<","<<m_Pts[i].z<<endl;  
  143.     }  
  144.     for (int i=0;i<(int)m_Edges.size();i++)  
  145.     {  
  146.         outfile<<"L,"<<m_Edges[i].left+1<<","<<m_Edges[i].right+1<<endl;  
  147.     }  
  148.     outfile.close();  
  149. }  
  150. void GetPoint(double &xx,double &yy,double &zz,string line)//从字符串line中解析出点的x,y,z坐标  
  151. {  
  152.     int flag=0;  
  153.     string tmp="";  
  154.     char *cstr;  
  155.     for (int i=(int)line.find(',')+1;i<(int)line.size();i++)  
  156.     {  
  157.         if (line[i]==',')  
  158.         {  
  159.             cstr=new char[tmp.size()+1];  
  160.             strcpy(cstr,tmp.c_str());  
  161.             if (flag==0) {xx=atof(cstr);tmp.resize(0);flag++;}  
  162.             else if (flag==1) {yy=atof(cstr);tmp.resize(0);flag++;}  
  163.             continue;  
  164.         }  
  165.         tmp=tmp+line[i];  
  166.     }  
  167.     if (fabs(xx)<1.0e-6) xx=0.0;  
  168.     if (fabs(yy)<1.0e-6) yy=0.0;  
  169.     if (fabs(zz)<1.0e-6) zz=0.0;  
  170. }  
  171. void Delaunay::Delete_Frame()//删除外边框  
  172. {  
  173.     EdgeArray BoundEdges;  
  174.     for (int i=0;i<4;i++) m_Pts.erase(m_Pts.begin());  
  175.     for (int i=0;i<(int)m_Tris.size();i++)  
  176.     {  
  177.         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)  
  178.         {  
  179.             DelTriangle(i,BoundEdges);  
  180.             BoundEdges.resize(0);  
  181.             i--;  
  182.         }  
  183.         else  
  184.         {  
  185.             for (int j=0;j<3;j++)  
  186.             {  
  187.                 m_Tris[i].v[j]-=4;  
  188.                 m_Tris[i].s[j].left-=4;  
  189.                 m_Tris[i].s[j].right-=4;  
  190.             }  
  191.         }  
  192.     }  
  193.     for (int i=0;i<4;i++) m_Edges.erase(m_Edges.begin());  
  194.     for (int i=0;i<(int)m_Edges.size();i++)  
  195.     {  
  196.         m_Edges[i].left-=4;  
  197.         m_Edges[i].right-=4;  
  198.     }  
  199. }  
  200.   
  201. void Delaunay::Boundary_Recover(int fromPoint,int toPoint)//恢复由指定点组成的边界  
  202. {  
  203.     EdgeArray BoundEdges;  
  204.     for (int i=0;i<(int)m_Tris.size();i++)  
  205.     {  
  206.         if (m_Tris[i].v[0]>=(fromPoint-1)&&m_Tris[i].v[2]<=(toPoint-1))  
  207.         {  
  208.             DelTriangle(i,BoundEdges);  
  209.             BoundEdges.resize(0);  
  210.             i--;  
  211.         }  
  212.     }  
  213. }  

main.cpp

  1. #include "delaunay.h"   
  2. int main()  
  3. {  
  4.     ifstream infile("input.txt");//打开"input.txt"文件  
  5.     if (!infile)//判断文件是否正常打开  
  6.     {  
  7.         cout<<"Unable to input nodes!";  
  8.         exit(1);  
  9.     }  
  10.     string line;  
  11.     PointArray p;  
  12.     double xx,yy,zz;  
  13.     int nodeSize;  
  14.     for (int i=0;i<4;i++)//读入4外边框点  
  15.     {  
  16.         getline(infile,line);  
  17.         GetPoint(xx,yy,zz,line);  
  18.         Point tmp={xx,yy,zz};  
  19.         p.push_back(tmp);  
  20.     }  
  21.     Delaunay MyMesh(p[0],p[1],p[2],p[3]);//实例化Delaunay类  
  22.     getline(infile,line);//读入节点数,用于后面循环   
  23.     char *cstr;  
  24.     cstr=new char[line.size()+1];  
  25.     strcpy(cstr,line.c_str());  
  26.     nodeSize=atoi(cstr);  
  27.     for (int i=0;i<nodeSize;i++)//读入每个节点的坐标  
  28.     {  
  29.         getline(infile,line);  
  30.         GetPoint(xx,yy,zz,line);  
  31.         MyMesh.AddPoint(xx,yy,zz);  
  32.     }  
  33.     infile.close();  
  34.     MyMesh.Delete_Frame();//删除外边框   
  35.     MyMesh.Boundary_Recover(203,466);  
  36.     MyMesh.Boundary_Recover(467,487);  
  37.     MyMesh.Boundary_Recover(488,511);  
  38.     MyMesh.Boundary_Recover(512,537);//以上都是恢复指定边界  
  39.     MyMesh.output();//将相应ANSYS命令流输出  
  40.     return 0;  
  41. }  
[cpp] view plain copy
 在CODE上查看代码片派生到我的代码片
  1. #include "delaunay.h"  
  2. int main()  
  3. {  
  4.     ifstream infile("input.txt");//打开"input.txt"文件  
  5.     if (!infile)//判断文件是否正常打开  
  6.     {  
  7.         cout<<"Unable to input nodes!";  
  8.         exit(1);  
  9.     }  
  10.     string line;  
  11.     PointArray p;  
  12.     double xx,yy,zz;  
  13.     int nodeSize;  
  14.     for (int i=0;i<4;i++)//读入4外边框点  
  15.     {  
  16.         getline(infile,line);  
  17.         GetPoint(xx,yy,zz,line);  
  18.         Point tmp={xx,yy,zz};  
  19.         p.push_back(tmp);  
  20.     }  
  21.     Delaunay MyMesh(p[0],p[1],p[2],p[3]);//实例化Delaunay类  
  22.     getline(infile,line);//读入节点数,用于后面循环  
  23.     char *cstr;  
  24.     cstr=new char[line.size()+1];  
  25.     strcpy(cstr,line.c_str());  
  26.     nodeSize=atoi(cstr);  
  27.     for (int i=0;i<nodeSize;i++)//读入每个节点的坐标  
  28.     {  
  29.         getline(infile,line);  
  30.         GetPoint(xx,yy,zz,line);  
  31.         MyMesh.AddPoint(xx,yy,zz);  
  32.     }  
  33.     infile.close();  
  34.     MyMesh.Delete_Frame();//删除外边框  
  35.     MyMesh.Boundary_Recover(203,466);  
  36.     MyMesh.Boundary_Recover(467,487);  
  37.     MyMesh.Boundary_Recover(488,511);  
  38.     MyMesh.Boundary_Recover(512,537);//以上都是恢复指定边界  
  39.     MyMesh.output();//将相应ANSYS命令流输出  
  40.     return 0;  
  41. }  


测试一组数据后,得到结果:



 

编译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 文件

[plain] view plaincopy
  1. ####################################  
  2. # Stellar mesh improvement options #  
  3. ####################################  
  4.   
  5. # This file contains all the possible options that can currently be set for   
  6. # mesh improvement. Lines that begin with '#' are comments and are ignored.  
  7. # Other lines take the form 'option intval floatval stringval', where option  
  8. # is the name of the option, and intval floatval and stringval are the possible  
  9. # fields that can be used to set that option. If an option takes an int, only  
  10. # a value for int needs to be given. If it's a float, a dummy int should be  
  11. # inserted before the float. If it's a string, a dummy int and a dummy float  
  12. # should be inserted before the string. This clumsiness is because I don't like  
  13. # coding string processing in C, and this is the easiest way. Any unset options  
  14. # will assume their default values.  
  15.   
  16. # verbosity: bigger number means more verbose output of improvement.  
  17. # default = 1  
  18. verbosity 0  
  19.   
  20. # use color in verbose improvement output. default = 0  
  21. usecolor 1  
  22.   
  23. # just output the mesh unchanged and quit. default = 0  
  24. outputandquit 0  
  25.   
  26. ## quality measure options  
  27. # qualmeasure: selects which quality measure to use as an objective function  
  28. # for optimizing the tetrahedra. The quality measures are described in  
  29. # Chapter 2 of Bryan's dissertation. default = 0  
  30. # 0 = minimum sine of the dihedral angles (default)  
  31. # 1 = square root of radius ratio (circumradius divided by inradius)  
  32. # 2 = V / l_rms^3 (volume divided by cube of root-mean-squared edge length)  
  33. # 5 = minimum sine with biased obtuse angles  
  34. qualmeasure 5  
  35.   
  36. # sinewarpfactor: float. for qualmeasure 5 only; sets the factor by which  
  37. # obtuse angles are scaled relative to acute angles. Default is 0.75  
  38. sinewarpfactor 0 0.75  
  39.   
  40. ## termination options  
  41. # BOTH goal angles must be reached to terminate improvement  
  42. # goalanglemin: float. terminates improvement early if minimum angle reaches  
  43. # this value. default = 90.0 (which precludes early termination)  
  44. goalanglemin 0 90.0  
  45. # goalanglemax: float. terminates improvement early if maximum angle reaches  
  46. # this value. default = 90.0  
  47. goalanglemax 0 90.0  
  48.   
  49. ## smoothing options  
  50. # nonsmooth: enable optimization-based smoothing. default = 1  
  51. nonsmooth 1  
  52. # facetsmooth: enable smoothing of facet vertices. default = 1  
  53. facetsmooth 1  
  54. # segmentsmooth: enable smoothing of segment vertices. default = 1  
  55. segmentsmooth 1  
  56. # usequadrics: enable use of surface quadric error for smoothing fixed boundary  
  57. # vertices. WARNING: this will allow the domain shape to change slightly.  But  
  58. # even a little play sometimes yields much better meshes. default = 0  
  59. usequadrics 0  
  60. # quadricoffset: amount to start quadric error at before subtracting.  
  61. # See alpha in Section 3.2.5 of Bryan's dissertation. default = 0.8  
  62. quadricoffset 0 0.8  
  63. # quadricscale: amount to scale quadric error before subtracting from offset.  
  64. # See beta in Section 3.2.5 of Bryan's dissertation. default = 300.0  
  65. quadricscale 0 300.0  
  66.   
  67. ## topological transformation options  
  68. # flip22: enable 2-2 flips (for boundary tets). default = 1  
  69. flip22 1  
  70. # multifaceremoval: enable multi-face removal. singlefaceremoval might still  
  71. # be on. default = 1  
  72. multifaceremoval 1  
  73. # singlefaceremoval: enable single face removal (2-3 and 2-2 flips). Has  
  74. # no effect when multifaceremoval is enabled. default = 1  
  75. singlefaceremoval 1  
  76. # edgeremoval: enable edge removal. default = 1  
  77. edgeremoval 1  
  78.   
  79. ## edge contraction options  
  80. # edgecontraction: enable edge contraction. default = 1  
  81. edgecontraction 1  
  82.   
  83. ## vertex insertion options  
  84. # enableinsert: enable ALL vertex insertion (overrides others). default = 1  
  85. enableinsert 0  
  86. # insertbody: enable just vertex insertion in body (interior). default = 1  
  87. insertbody 0  
  88. # insertfacet: enable just insertion on facets. default = 1  
  89. insertfacet 0  
  90. # insertsegment: enable just insertion on segments. default = 1  
  91. insertsegment 0  
  92. # insertthreshold: on each insertion pass, try vertex insertion in this  
  93. fraction of the tetrahedra. default = 0.031 (the worst 3.1%)  
  94. insertthreshold 0 0.031  
  95.   
  96. ## size control options  
  97. # (See Chapter 6 of Bryan's dissertation.)  
  98. # sizing: enable control of element sizes. default = 0  
  99. sizing 0  
  100. # sizingpass: enable edge length correction before quality improvement.  
  101. # default = 0  
  102. sizingpass 0  
  103. # targetedgelength: the target edge length for this mesh. If set to 0.0, the  
  104. # target edge length is initialized automatically to the initial mean edge  
  105. # length. default = 0.0  
  106. targetedgelength 0 0.0  
  107. # longerfactor: factor by which an edge can be longer than the target edge  
  108. # length before being considered "too long". default = 3.0  
  109. longerfactor 0 2.0  
  110. # shorterfactor: factor by which an edge can be shorter than the target edge  
  111. # length before being considered "too short" default = 0.33  
  112. shorterfactor 0 0.50  
  113.   
  114. ## anisotropy options  
  115. # (See Chapter 7 of Bryan's dissertation.)  
  116. # anisotropic: enable anisotropic meshing. default = 0  
  117. anisotropic 0  
  118. # tensor: which size/anisotropy tensor to use. default = 0  
  119. # 0 = identity  
  120. # 1 = stretch x  
  121. # 2 = stretch y  
  122. # 3 = sink  
  123. # 4 = swirl  
  124. # 5 = center  
  125. # 6 = perimeter  
  126. # 7 = right  
  127. # 8 = sine  
  128. tensor 6  
  129.   
  130. ## quality file output options  
  131. # These files list, for each tetrahedron, the values of the quality measures  
  132. # minsineout: enable output of .minsine quality file. default = 1  
  133. minsineout 1  
  134. # minangout: enable output of .minang quality file. default = 0  
  135. minangout 0  
  136. # maxangout: enable output of .maxang quality file. default = 0  
  137. maxangout 0  
  138. # vlrmsout: enable output of .vlrms quality file. default = 0  
  139. vlrmsout 0  
  140. # nrrout: enable output of the .nrr quality file. default = 0  
  141. nrrout 0  
  142.   
  143. ## animation options  
  144. # animate: activate animation file output (a set of output files after each  
  145. # pass). default = 0  
  146. animate 0  
  147. # timeseries: when animate = 1, only output .stats. default = 0  
  148. timeseries 0  
  149.   
  150. ## output filename option  
  151. # fileprefix: filename prefix that distinguishes the output files from the  
  152. # input files. If none specified, an iteration number is appended to the input  
  153. # filenames.  
  154. #fileprefix 0 5 ../output/testprefix  
[plain] view plain copy
 在CODE上查看代码片派生到我的代码片
  1. ####################################  
  2. # Stellar mesh improvement options #  
  3. ####################################  
  4.   
  5. # This file contains all the possible options that can currently be set for   
  6. # mesh improvement. Lines that begin with '#' are comments and are ignored.  
  7. # Other lines take the form 'option intval floatval stringval', where option  
  8. # is the name of the option, and intval floatval and stringval are the possible  
  9. # fields that can be used to set that option. If an option takes an int, only  
  10. # a value for int needs to be given. If it's a float, a dummy int should be  
  11. # inserted before the float. If it's a string, a dummy int and a dummy float  
  12. # should be inserted before the string. This clumsiness is because I don't like  
  13. # coding string processing in C, and this is the easiest way. Any unset options  
  14. # will assume their default values.  
  15.   
  16. # verbosity: bigger number means more verbose output of improvement.  
  17. # default = 1  
  18. verbosity 0  
  19.   
  20. # use color in verbose improvement output. default = 0  
  21. usecolor 1  
  22.   
  23. # just output the mesh unchanged and quit. default = 0  
  24. outputandquit 0  
  25.   
  26. ## quality measure options  
  27. # qualmeasure: selects which quality measure to use as an objective function  
  28. # for optimizing the tetrahedra. The quality measures are described in  
  29. # Chapter 2 of Bryan's dissertation. default = 0  
  30. # 0 = minimum sine of the dihedral angles (default)  
  31. # 1 = square root of radius ratio (circumradius divided by inradius)  
  32. # 2 = V / l_rms^3 (volume divided by cube of root-mean-squared edge length)  
  33. # 5 = minimum sine with biased obtuse angles  
  34. qualmeasure 5  
  35.   
  36. # sinewarpfactor: float. for qualmeasure 5 only; sets the factor by which  
  37. # obtuse angles are scaled relative to acute angles. Default is 0.75  
  38. sinewarpfactor 0 0.75  
  39.   
  40. ## termination options  
  41. # BOTH goal angles must be reached to terminate improvement  
  42. # goalanglemin: float. terminates improvement early if minimum angle reaches  
  43. # this value. default = 90.0 (which precludes early termination)  
  44. goalanglemin 0 90.0  
  45. # goalanglemax: float. terminates improvement early if maximum angle reaches  
  46. # this value. default = 90.0  
  47. goalanglemax 0 90.0  
  48.   
  49. ## smoothing options  
  50. # nonsmooth: enable optimization-based smoothing. default = 1  
  51. nonsmooth 1  
  52. # facetsmooth: enable smoothing of facet vertices. default = 1  
  53. facetsmooth 1  
  54. # segmentsmooth: enable smoothing of segment vertices. default = 1  
  55. segmentsmooth 1  
  56. # usequadrics: enable use of surface quadric error for smoothing fixed boundary  
  57. # vertices. WARNING: this will allow the domain shape to change slightly.  But  
  58. # even a little play sometimes yields much better meshes. default = 0  
  59. usequadrics 0  
  60. # quadricoffset: amount to start quadric error at before subtracting.  
  61. # See alpha in Section 3.2.5 of Bryan's dissertation. default = 0.8  
  62. quadricoffset 0 0.8  
  63. # quadricscale: amount to scale quadric error before subtracting from offset.  
  64. # See beta in Section 3.2.5 of Bryan's dissertation. default = 300.0  
  65. quadricscale 0 300.0  
  66.   
  67. ## topological transformation options  
  68. # flip22: enable 2-2 flips (for boundary tets). default = 1  
  69. flip22 1  
  70. # multifaceremoval: enable multi-face removal. singlefaceremoval might still  
  71. # be on. default = 1  
  72. multifaceremoval 1  
  73. # singlefaceremoval: enable single face removal (2-3 and 2-2 flips). Has  
  74. # no effect when multifaceremoval is enabled. default = 1  
  75. singlefaceremoval 1  
  76. # edgeremoval: enable edge removal. default = 1  
  77. edgeremoval 1  
  78.   
  79. ## edge contraction options  
  80. # edgecontraction: enable edge contraction. default = 1  
  81. edgecontraction 1  
  82.   
  83. ## vertex insertion options  
  84. # enableinsert: enable ALL vertex insertion (overrides others). default = 1  
  85. enableinsert 0  
  86. # insertbody: enable just vertex insertion in body (interior). default = 1  
  87. insertbody 0  
  88. # insertfacet: enable just insertion on facets. default = 1  
  89. insertfacet 0  
  90. # insertsegment: enable just insertion on segments. default = 1  
  91. insertsegment 0  
  92. # insertthreshold: on each insertion pass, try vertex insertion in this  
  93. fraction of the tetrahedra. default = 0.031 (the worst 3.1%)  
  94. insertthreshold 0 0.031  
  95.   
  96. ## size control options  
  97. # (See Chapter 6 of Bryan's dissertation.)  
  98. # sizing: enable control of element sizes. default = 0  
  99. sizing 0  
  100. # sizingpass: enable edge length correction before quality improvement.  
  101. # default = 0  
  102. sizingpass 0  
  103. # targetedgelength: the target edge length for this mesh. If set to 0.0, the  
  104. # target edge length is initialized automatically to the initial mean edge  
  105. # length. default = 0.0  
  106. targetedgelength 0 0.0  
  107. # longerfactor: factor by which an edge can be longer than the target edge  
  108. # length before being considered "too long". default = 3.0  
  109. longerfactor 0 2.0  
  110. # shorterfactor: factor by which an edge can be shorter than the target edge  
  111. # length before being considered "too short" default = 0.33  
  112. shorterfactor 0 0.50  
  113.   
  114. ## anisotropy options  
  115. # (See Chapter 7 of Bryan's dissertation.)  
  116. # anisotropic: enable anisotropic meshing. default = 0  
  117. anisotropic 0  
  118. # tensor: which size/anisotropy tensor to use. default = 0  
  119. # 0 = identity  
  120. # 1 = stretch x  
  121. # 2 = stretch y  
  122. # 3 = sink  
  123. # 4 = swirl  
  124. # 5 = center  
  125. # 6 = perimeter  
  126. # 7 = right  
  127. # 8 = sine  
  128. tensor 6  
  129.   
  130. ## quality file output options  
  131. # These files list, for each tetrahedron, the values of the quality measures  
  132. # minsineout: enable output of .minsine quality file. default = 1  
  133. minsineout 1  
  134. # minangout: enable output of .minang quality file. default = 0  
  135. minangout 0  
  136. # maxangout: enable output of .maxang quality file. default = 0  
  137. maxangout 0  
  138. # vlrmsout: enable output of .vlrms quality file. default = 0  
  139. vlrmsout 0  
  140. # nrrout: enable output of the .nrr quality file. default = 0  
  141. nrrout 0  
  142.   
  143. ## animation options  
  144. # animate: activate animation file output (a set of output files after each  
  145. # pass). default = 0  
  146. animate 0  
  147. # timeseries: when animate = 1, only output .stats. default = 0  
  148. timeseries 0  
  149.   
  150. ## output filename option  
  151. # fileprefix: filename prefix that distinguishes the output files from the  
  152. # input files. If none specified, an iteration number is appended to the input  
  153. # filenames.  
  154. #fileprefix 0 5 ../output/testprefix  


再次运行,

./Stellar -s conf model.2

运行结果:



顶点从6000多降到了5000多,用tetiew来查看:



展开阅读全文

没有更多推荐了,返回首页