使用gmsh作为三维网格生成器 7.26

Gmsh是一个免费的带有内置前后期处理机制的三维有限元网格生成器。其设计的目标是要提供一个快速轻便的具有可控参数功能和先进的可视化能力的网格生成工具。Gmsh主要围绕四个单元：几何，网格，求解和后处理。这些可控参数的输入可以在交互式的图形界面方式，或在ASCII文本文件中使用gmsh自己的脚本语言得以实现。

打算采用OCCT建模->gmsh网格生成与求解->VTK显示的方式来做一个框架

 

OCC->gmsh几何信息输入：

gmsh使用occ作为几何内核 相应的操作类都在namespace gmsh::model::occ中 熟悉相关命令之后可以很方便的直接绘图

也可以直接传TopoDS_Shape值作为几何信息 亦支持STL BREP STEP和IGES等格式的输入

1 STL

STL文件用于表示三角形网格的文件格式 其文件格式非常简单 应用很广泛 是最多快速原型系统所应用的标准文件类型

STL的OCCT生成比较麻烦 原有的类是对三角化面片类Handle_PolyTriangulation进行操作 只能写一个面片的信息

方法为遍历三角化的各个网格 然后输出相应的几何信息

改写后传参为shape 对其中每个面进行遍历三角化后写信息

Standard_Boolean writeASCII(const TopoDS_Shape& shape,FILE* theFile,Handle(Message_ProgressIndicator)& theProgInd, double deflection)
{
	static const Standard_Integer IND_THRESHOLD = 1000; // increment the indicator every 1k triangles
	// note that space after 'solid' is necessary for many systems
	if (fwrite("solid \n", 1, 7, theFile) != 7)
	{
		return Standard_False;
	}
	for (TopExp_Explorer aFaceExplorer(shape, TopAbs_FACE); aFaceExplorer.More(); aFaceExplorer.Next()) {
		TopoDS_Face face = TopoDS::Face(aFaceExplorer.Current());
		TopLoc_Location location;
		BRepMesh_IncrementalMesh(face, deflection);
		Handle_Poly_Triangulation theMesh = BRep_Tool::Triangulation(face, location);
		char aBuffer[512];
		memset(aBuffer, 0, sizeof(aBuffer));

		Message_ProgressSentry aPS(theProgInd, "Triangles", 0,
			theMesh->NbTriangles(), IND_THRESHOLD);

		const TColgp_Array1OfPnt& aNodes = theMesh->Nodes();
		const Poly_Array1OfTriangle& aTriangles = theMesh->Triangles();
		const Standard_Integer NBTriangles = theMesh->NbTriangles();
		Standard_Integer anElem[3] = { 0, 0, 0 };
		for (Standard_Integer aTriIter = 1; aTriIter <= NBTriangles; ++aTriIter)
		{
			const Poly_Triangle& aTriangle = aTriangles(aTriIter);
			aTriangle.Get(anElem[0], anElem[1], anElem[2]);

			const gp_Pnt aP1 = aNodes(anElem[0]);
			const gp_Pnt aP2 = aNodes(anElem[1]);
			const gp_Pnt aP3 = aNodes(anElem[2]);

			const gp_Vec aVec1(aP1, aP2);
			const gp_Vec aVec2(aP1, aP3);
			gp_Vec aVNorm = aVec1.Crossed(aVec2);
			if (aVNorm.SquareMagnitude() > gp::Resolution())
			{
				aVNorm.Normalize();
			}
			else
			{
				aVNorm.SetCoord(0.0, 0.0, 0.0);
			}

			Sprintf(aBuffer,
				" facet normal % 12e % 12e % 12e\n"
				"   outer loop\n"
				"     vertex % 12e % 12e % 12e\n"
				"     vertex % 12e % 12e % 12e\n"
				"     vertex % 12e % 12e % 12e\n"
				"   endloop\n"
				" endfacet\n",
				aVNorm.X(), aVNorm.Y(), aVNorm.Z(),
				aP1.X(), aP1.Y(), aP1.Z(),
				aP2.X(), aP2.Y(), aP2.Z(),
				aP3.X(), aP3.Y(), aP3.Z());

			if (fprintf(theFile, "%s", aBuffer) < 0)
			{
				return Standard_False;
			}

			// update progress only per 1k triangles
			if ((aTriIter % IND_THRESHOLD) == 0)
			{
				aPS.Next();
			}
		}
	}
	if (fwrite("endsolid\n", 1, 9, theFile) != 9)
	{
		return Standard_False;
	}

	return Standard_True;
}

Standard_Boolean writeSTL(const TopoDS_Shape& shape, const TCollection_AsciiString path,double deflection = 0.01) {

	FILE* aFile = OSD_OpenFile(path, "w");
	Handle(Message_ProgressIndicator) theProgInd;

	if (aFile == NULL)
	{
		return Standard_False;
	}

	Standard_Boolean isOK = writeASCII(shape, aFile, theProgInd, deflection);
	fclose(aFile);
	return isOK;
}

 

gmsh读取STL文件后可以直接显示 也可以进行remesh操作 在此处意义不大 既然看了写了 就留作他用

2 TopoDS_Shape输入和BREP STEP和IGES等格式的输入

前者的类为：

importShapesNativePointer(TopoDS_Shape& shape, vector<pair<int, int>& dimTags,bool highestDimOnly)

BREP STEP和IGES文件输入的类为：

importShapes(const string filename, vector<pair<int, int>& dimTags,bool highestDimOnly)

参数分别为Shape/文件名称（通过文件名称进行字符串分割识别出相应的文件类型 采用不同的读文件方法）维数和下标信息 以及是否仅生成最高阶几何信息（默认为否）

其源码均为调用重载的GModelOI_OCC::importShapes方法 观察其源码得知 在读取IGES与STEP的情况下文件读写操作比较复杂 而同时在OCCT中直接写Brep文件也相对比较简单

void generateBrep(TopoDS_Shape& shape)
{
	ofstream dumpFile("test.brep");
	BRepTools::Write(shape, dumpFile);
}

不妨使用Brep来作为OCCT和gmsh的桥梁

Brep参考资料：https://www.cnblogs.com/opencascade/p/3485621.html

 

gmsh->VTK网格信息的显示

原以为要写一个vtk的函数去读.msh文件 后面发现可以直接输出.vtk格式的文件 就比较方便

从文件头可以看到 DATASET为UNSTRUCTURED_GRID  那么在vtk中的数据流就应该是

vtkUnstructuredGrid->vtkAlgorithmOutput->vtkPolyData（Filter）->渲染管线

参考资料：

https://lorensen.github.io/VTKExamples/site/Cxx/IO/ReadLegacyUnstructuredGrid/

https://vtk.org/doc/nightly/html/classvtkUnstructuredGridAlgorithm.html

void viewInVtk(const char* vtkFile) {
	vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New();
	vtkSmartPointer<vtkUnstructuredGridReader> gridReader= vtkSmartPointer<vtkUnstructuredGridReader>::New();
	gridReader->SetFileName(vtkFile);
	gridReader->Update();
	
	/*
	grid->polydata 
	using class vtkPolyDataAlgorithm
	vtkAlgorithmOutput* GetOutputPort()
	vtkPolyData* GetOutput()
	*/
	vtkSmartPointer<vtkSphereSource> sphere =vtkSmartPointer<vtkSphereSource>::New();
	sphere->SetPhiResolution(10);
	sphere->SetThetaResolution(10);
	sphere->SetRadius(.08);
	vtkSmartPointer<vtkGlyph3DMapper> pointMapper =vtkSmartPointer<vtkGlyph3DMapper>::New();
	pointMapper->SetInputConnection(gridReader->GetOutputPort());
	pointMapper->SetSourceConnection(sphere->GetOutputPort());

	vtkSmartPointer<vtkExtractEdges> extractEdges = vtkSmartPointer<vtkExtractEdges>::New();
	extractEdges->SetInputConnection(gridReader->GetOutputPort());
	vtkSmartPointer<vtkTubeFilter> tubesFilter =vtkSmartPointer<vtkTubeFilter>::New();
	tubesFilter->SetInputConnection(extractEdges->GetOutputPort());
	tubesFilter->SetRadius(.03);
	tubesFilter->SetNumberOfSides(10);
	vtkSmartPointer<vtkPolyDataMapper> tubeMapper =vtkSmartPointer<vtkPolyDataMapper>::New();
	tubeMapper->SetInputConnection(tubesFilter->GetOutputPort());
	tubeMapper->Update();

	vtkSmartPointer<vtkDataSetTriangleFilter> triangleFilter = vtkSmartPointer<vtkDataSetTriangleFilter>::New();
	triangleFilter->SetInputConnection(gridReader->GetOutputPort());
	vtkSmartPointer<vtkDataSetMapper> geometryMapper =vtkSmartPointer<vtkDataSetMapper>::New();
	geometryMapper->SetInputConnection(triangleFilter->GetOutputPort());

	vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
	actor->SetMapper(tubeMapper);
	vtkSmartPointer<vtkActor> actor1 = vtkSmartPointer<vtkActor>::New();
	actor1->SetMapper(geometryMapper);

	vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
	renderer->AddActor(actor);
	renderer->AddActor(actor1);

	vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
	renderWindow->SetSize(1200, 900);
	renderWindow->AddRenderer(renderer);

	vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
	renderWindowInteractor->SetRenderWindow(renderWindow);

	renderer->SetBackground(0, 0.2, 0.8);

	renderWindow->Render();
	renderWindowInteractor->Start();
}

int main() {
	generateBrep(MakeBottle(20,30,8));
	generateMesh("test.brep");
	viewInVtk("test.vtk");
	return 0;
}

 

gmsh的相关操作：

1.由于gmsh SDK是使用gcc编译 在使用msvc的时候 要使用C API而不是 C++ API 具体做法是将gmsh.h_cwrap重命名为gmsh.h并替代原gmsh.h文件 可能会带来效率上的轻微损失 （另外有些情况下内存报错不知道是不是和这个有关系）

2.查阅GMSH用户手册 在附录B与附录D中写的比较简略 结合其自带的gui测试一下就明白了

3.其设置通过option类，采用名称（string）+值（double/string）的形式进行表达

例子如下：

void generateMesh(const std::string filename) {
	
	using namespace gmsh;

	double lc = 1;
	initialize();

	option::setNumber("General.Terminal", 1);
	option::setNumber("General.Verbosity", 99);
	option::setNumber("Mesh.Algorithm", 2);
	option::setNumber("Mesh.Algorithm3D", 1);
	option::setNumber("Mesh.CharacteristicLengthMax", lc);
	option::setNumber("Mesh.SurfaceFaces", 1);


	std::vector<pair<int, int>> dimTags;
	model::occ::importShapes("test.brep", dimTags);
	model::occ::synchronize();


	model::mesh::setSize(dimTags, 0.1);
	model::mesh::generate(3);

	write("test.msh");
	write("test.vtk");

	finalize();

}

效果图

gmsh fltk显示

vtk中显示

 

msh文件 在gmsh手册第九章上有详解

从entity行可以看出 从OCC中导入的0 1 2 3维实体数目 （vertex  edge face solid） 这也决定了所能生成网格的最大维数

0-3维网格信息