文章目录
前言
Gmsh的安装与入门教程可以去看我的前一篇文章有限元剖网格之Gmsh安装与使用入门。本篇文章主要介绍如何借助Gmsh的C++ API剖二维网格(即三角形)。
一、点，三角形的定义
利用Gmsh剖的二维网格，一般是为了后续有限元计算的需要，后面也会给出一篇有限元的文章。此处主要是为了Gmsh前期的输入，由点连线构成多边形，以及最后利用Gmsh得到网格后对三角形信息的存储。
1. 点的定义
class Point
{
public:
Point(const double x = 0, const double y = 0, const double z = 0):_x(x), _y(y), _z(z) { }
~Point() { }
const double getX() const {return _x;}
const double getY() const {return _y;}
const double getZ() const {return _z;}
void setX(const double x) {_x = x;}
void setY(const double y) {_y = y;}
void setZ(const double z) {_z = z;}
bool operator < (const Point& p) const;
protected:
double _x, _y, _z;
};
bool Point::operator < (const Point& p) const
{
double px = p.getX();
double py = p.getY();
if(px == _x) {
return py < _y;
} else {
return px < _x;
}
}
2. 三角形的定义
class Triangle
{
/**
* Triangle vertex index
* _id[0] ~ _id[2] : three vertices index of triangle
* _id[3] is the midpoint between _id[0] and _id[1]
* _id[4] is the midpoint between _id[1] and _id[2]
* _id[5] is the midpoint between _id[2] and _id[0]
*/
public:
Triangle() {}
~Triangle() {}
const vector<unsigned int>& getIds() const {return _ids;}
void clear() {_ids.clear();}
void append(unsigned int id) {_ids.push_back(id);}
void resize(const unsigned int size) {_ids.resize(size);}
private:
vector<unsigned int> _ids;
};
注释也注明了，容器前三个空间就存储了三角形三个顶点的索引，后面存储的是三条边的中点的索引，是为了有限元(Lagrange 二次元)计算方便存储的。
二、区域构造
在此我们可以自己构造一个由多个多边形组成的区域。如下是一个简单的case：
void generateBoudaryAndHoles(vector<Point>& boundary, vector<vector<Point> >& holes)
{
double left, right, bot, top;
double x, y, z;
left = 0;
right = 5;
bot = 0;
top = 4;
z = 0;
Point p;
p = Point(left, bot, z);
boundary.push_back(p);
p = Point(right, bot, z);
boundary.push_back(p);
p = Point(right, top, z);
boundary.push_back(p);
p = Point(left, top, z);
boundary.push_back(p);
p = Point(left, (bot+top)/2.0, z);
boundary.push_back(p);
vector<Point> hole;
left = 1;
right = 2;
bot = 1;
top = 2;
p = Point(left, bot, z);
hole.push_back(p);
p = Point(right, bot, z);
hole.push_back(p);
p = Point(right, top, z);
hole.push_back(p);
p = Point(left, top, z);
hole.push_back(p);
p = Point(left, (bot+top)/2.0, z);
hole.push_back(p);
holes.push_back(hole);
hole.clear();
left = 3;
right = 4;
bot = 1;
top = 2;
p = Point(left, bot, z);
hole.push_back(p);
p = Point(right, bot, z);
hole.push_back(p);
p = Point(right, top, z);
hole.push_back(p);
p = Point(left, top, z);
hole.push_back(p);
holes.push_back(hole);
}
其中boundary存储的是区域的边界点信息，holes存储的是内部所有洞的点的信息。
注意：所有多边形点的信息都是逆时针顺序的
三、调用Gmsh的API剖网格并得到网格信息
以下为main函数中的主要代码 ：调用Gmsh的API。相应位置都有注释，获取所有点的信息，获取单元信息，获取边界上三角形信息，获取洞边上三角形信息。
需要注意：程序中我设置gmsh内点标号即下标是从1开始的，此处我把点和单元信息存储在vector容器中，下标则从0开始，所以在代码的相应位置有-1操作。
1. 剖网格部分代码
int main(int argc, char **argv)
{
gmsh::initialize();
gmsh::option::setNumber("General.Terminal", 1);
gmsh::model::add("demo");
double boundartLc = 1, holeLc = 0.4;
vector<Point> boundary;
vector<vector<Point> > holes;
generateBoudaryAndHoles(boundary, holes);
int dim = 1;
int tag = 1;
double x, y, z;
vector<int> loop, curveLoop, planeSurface;
int boundaryTag;
for(int i = 0; i < boundary.size(); i++)
{
x = boundary[i].getX();
y = boundary[i].getY();
z = boundary[i].getZ();
loop.push_back(gmsh::model::geo::addPoint(x, y, z, boundartLc));
}
for(int i = 0; i < loop.size(); i++)
{
curveLoop.push_back(gmsh::model::geo::addLine(loop[i], i == loop.size()-1 ? loop[0] : loop[i+1]));
}
planeSurface.push_back(gmsh::model::geo::addCurveLoop(curveLoop));
boundaryTag = gmsh::model::addPhysicalGroup(dim, curveLoop);
gmsh::model::setPhysicalName(dim, boundaryTag, "Boundary");
vector<int> bcurveLoop = curveLoop;
vector<vector<int> > recCurveLoop;
vector<int> holesTags;
for(int i = 0; i < holes.size(); i++)
{
loop.clear();
curveLoop.clear();
for(int j = 0; j < holes[i].size(); j++)
{
x = holes[i][j].getX();
y = holes[i][j].getY();
z = holes[i][j].getZ();
loop.push_back(gmsh::model::geo::addPoint(x, y, z, holeLc));
}
for(int j = 0; j < loop.size(); j++)
{
curveLoop.push_back(gmsh::model::geo::addLine(loop[j], j == loop.size()-1 ? loop[0] : loop[j+1]));
}
planeSurface.push_back(-gmsh::model::geo::addCurveLoop(curveLoop));
tag = gmsh::model::addPhysicalGroup(dim, curveLoop);
holesTags.push_back(tag);
string name = "hole" + to_string(i);
gmsh::model::setPhysicalName(dim, tag, name);
recCurveLoop.push_back(curveLoop);
}
gmsh::model::geo::addPlaneSurface(planeSurface);
dim = 2;
gmsh::model::addPhysicalGroup(dim, planeSurface);
gmsh::model::geo::synchronize();
gmsh::model::mesh::generate(2);
gmsh::write("demo.msh");
gmsh::finalize();
return 0;
}
2. 获取网格的所有点的代码
// Nodes
dim = -1;
tag = -1;
bool includeBoundary = false;
vector<double> coord, parametricCoord;
vector<size_t> nodeTags;
gmsh::model::mesh::getNodes(nodeTags, coord, parametricCoord, dim, tag, includeBoundary, true);
cout<<"The number of nodes: "<<coord.size()/3<<endl;
vector<Point> nodes;
nodes.reserve(1000);
Point p;
for(int i = 0; i < coord.size(); i+=3)
{
x = coord[i];
y = coord[i+1];
z = coord[i+2];
p = Point(x, y, z);
nodes.push_back(p);
}
3. 获取网格的所有三角形的代码
//Elements
dim = -1;
tag = -1;
vector<int> elementTypes;
vector<vector<size_t> > elementTags, nodeTagss;
gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTagss, dim, tag);
cout<<"The number of elements: "<<nodeTagss[1].size()/3<<endl;
// boundaryNodesId include holes boundary nodes id
// vector<size_t> boundaryNodesId = elementTags[0];
vector<Triangle> triangles;
triangles.reserve(1000);
unsigned int id1, id2, id3;
for(int i = 0; i < nodeTagss[1].size(); i += 3)
{
Triangle tri;
id1 = nodeTagss[1][i]-1;
id2 = nodeTagss[1][i+1]-1;
id3 = nodeTagss[1][i+2]-1;
tri.append(id1);
tri.append(id2);
tri.append(id3);
triangles.push_back(tri);
}
4. 获取最外围边界边上和内部洞的边上的点的代码
// get boundary nodes id and it's coordinate information
// coord is coordinate information (x1, y1, z1, .....)
dim = 1;
vector<size_t> bouNodesIds;
gmsh::model::mesh::getNodesForPhysicalGroup(dim, boundaryTag, bouNodesIds, coord);
vector<unsigned int> boundaryNodesIds(bouNodesIds.size());
// just because index from start 0
for(int i = 0; i < bouNodesIds.size(); i++)
{
boundaryNodesIds[i] = bouNodesIds[i]-1;
}
// get holes nodes id and it's coordinate information
vector<vector<size_t> > hosNodesIds(holesTags.size());
vector<vector<unsigned int> > holesNodesIds(holesTags.size());
for(int i = 0; i < holesTags.size(); i++)
{
gmsh::model::mesh::getNodesForPhysicalGroup(dim, holesTags[i], hosNodesIds[i], coord);
}
for(int i = 0; i < hosNodesIds.size(); i++)
{
holesNodesIds[i].resize(hosNodesIds[i].size());
for(int j = 0; j < hosNodesIds[i].size(); j++)
{
holesNodesIds[i][j] = hosNodesIds[i][j]-1;
}
}
#if 0
for(int i = 0; i < holesNodesIds.size(); i++)
{
for(int j = 0; j < holesNodesIds[i].size(); j++)
{
cout<<holesNodesIds[i][j]<<" ";
}
cout<<endl;
}
#endif
5. 获取内部洞的边上的三角形的代码
// get holes arround triangle id
vector<vector<int> > holesElementsIds(recCurveLoop.size());
dim = 1;
int dim2 = 2;
vector<size_t> elementTag;
bool strict = true;
double xx, yy, zz;
vector<int> holeElementsIds;
for(int i = 0; i < recCurveLoop.size(); i++)
{
holeElementsIds.clear();
for(int j = 0; j < recCurveLoop[i].size(); j++)
{
//elementTypes.clear();
//elementTags.clear();
//nodeTagss.clear();
tag = recCurveLoop[i][j];
gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTagss, dim, tag);
for(int k = 0; k < nodeTagss[0].size(); k += 2)
{
id1 = nodeTagss[0][k]-1;
id2 = nodeTagss[0][k+1]-1;
xx = (nodes[id1].getX() + nodes[id2].getX())/2.0;
yy = (nodes[id1].getY() + nodes[id2].getY())/2.0;
zz = (nodes[id1].getZ() + nodes[id2].getZ())/2.0;
gmsh::model::mesh::getElementsByCoordinates(xx, yy, zz, elementTag, dim2, strict);
for(int kk = 0; kk < elementTag.size(); kk++)
{
holeElementsIds.push_back(elementTag[kk]-1);
}
}
}
holesElementsIds.push_back(holeElementsIds);
}
#if 0
for(int i = 0; i < holesElementsIds.size(); i++)
{
for(int j = 0; j < holesElementsIds[i].size(); j++)
{
cout<<holesElementsIds[i][j]<<" ";
}
cout<<endl;
}
#endif
注：这里指的三角形是，三角形的有一条边和内部洞的边重合(即三角形有两个顶点落于洞的边上)。
6. 获取最外围的边上的三角形的代码
// get boundary arround triangle id
vector<int> boundaryElementsIds;
dim = 1;
dim2 = 2;
strict = true;
for(int i = 0; i < bcurveLoop.size(); i++)
{
tag = bcurveLoop [i];
gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTagss, dim, tag);
for(int k = 0; k < nodeTagss[0].size(); k += 2)
{
id1 = nodeTagss[0][k]-1;
id2 = nodeTagss[0][k+1]-1;
xx = (nodes[id1].getX() + nodes[id2].getX())/2.0;
yy = (nodes[id1].getY() + nodes[id2].getY())/2.0;
zz = (nodes[id1].getZ() + nodes[id2].getZ())/2.0;
gmsh::model::mesh::getElementsByCoordinates(xx, yy, zz, elementTag, dim2, strict);
for(int kk = 0; kk < elementTag.size(); kk++)
{
boundaryElementsIds.push_back(elementTag[kk]-1);
}
}
}
注：这里指的三角形是，三角形的有一条边和最外围的边重合(即三角形有两个顶点落于最外围的边上)。
四、程序运行结果
1. 程序输出
如上图所示，我们可以很直观的看到点和单元三角形的数量，共有95个点，149个三角形。其中gmsh的info的log信息提到的206个单元是因为他把一维的单元(线段)也加上了。
程序也会在对应目录下产生一个“demo.msh”文件，用图形界面的gmsh软件打开，就能看到带有点标号的图。
2. 带有点标号的结果图
下图中标注的数字为点的编号，从1开始。
因为程序中设置了洞周边三角形尺寸比边界设置的小，所以得到的网格也是在洞的附近网格比较密。
3. 没有点标号的结果图
***注：关于效果图是否有点单元编号以及阴影填充效果，相关设置可以看我Gmsh剖一维网格教程附代码这篇文章，里面有介绍，都是在软件上方的菜单栏中。
五、总结
剖网格也只是前处理，关键也是收集到想要的网格信息，给后面的有限元引擎进行计算。一般是点，单元，以及相应的边界条件信息，在程序中也都有写，且都有注释。可以在程序中加入cout打印输出到终端，再结合图形界面，学习了解各函数API的功能。
以上涉及的完整代码可以去我的GitHub网站下载查看。后续会整理用Gmsh剖三维四面体网格的教程，欢迎大家一起讨论学习交流。
扫一扫
595
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
