nmWTAI-Platform/3rd/meshgen/example/mesh_example.cpp

#include "stdafx.h"
#include "singlePhaseSolver.h"
#include "windows.h"
#include "iostream"
#include "fstream"
void Ng2VTKWithData(Ng_Mesh *mesh, std::vector<double> vecField, std::string filename)
{
	std::cout << "Convert Ng_Mesh to vtk legacy file version of 4.2" << std::endl;
	int i, nseg, ne, np, matnum;
	int nodes[3];
	double point[2];
	int *EleNodeNum;
	double *NodeCoors;

	np = Ng_GetNP_2D(mesh);
	ne = Ng_GetNE_2D(mesh);
	nseg = Ng_GetNSeg_2D(mesh);

	std::ofstream vtkFile;
	vtkFile.open(filename);

	if (!vtkFile.is_open())
	{
		std::cerr << "Error opening file: " << "mesh.vtk" << std::endl;
		return;
	}

	// Write VTK file header
	vtkFile << "# vtk DataFile Version 4.1" << std::endl;
	vtkFile << "Mesh data" << std::endl;
	vtkFile << "ASCII" << std::endl;
	vtkFile << "DATASET UNSTRUCTURED_GRID" << std::endl;

	// Write points
	vtkFile << "POINTS " << np << " float" << std::endl;
	for (i = 1; i <= np; i++)
	{
		Ng_GetPoint_2D(mesh, i, point);
		vtkFile << point[0] << " " << point[1] << " " << 0 << std::endl;
	}

	// Write cells
	vtkFile << "CELLS " << ne << " " << ne * 4 << std::endl;
	std::vector<int> vecMater;
	for (int i = 1; i <= ne; i++)
	{
		Ng_GetElement_2D(mesh, i, nodes, &matnum);
		vecMater.emplace_back(matnum);
		vtkFile << 3 << " " << nodes[0] - 1 << " " << nodes[1] - 1 << " " << nodes[2] - 1 << std::endl;
	}

	// Write cell types
	vtkFile << "CELL_TYPES " << ne << std::endl;
	for (int i = 1; i <= ne; i++)
	{
		vtkFile << 5 << std::endl; // Assuming getCellType is a method that returns the VTK cell type
	}

	// Write cell data
	vtkFile << "CELL_DATA " << ne << std::endl;
	vtkFile << "SCALARS Material int 1" << std::endl;
	vtkFile << "LOOKUP_TABLE default" << std::endl;
	for (int i = 0; i < vecMater.size(); i++)
	{
		vtkFile << vecMater[i] << std::endl;
	}

	// Write point data
	vtkFile << "POINT_DATA " << np << std::endl;
	vtkFile << "SCALARS Pressure double 1" << std::endl;
	vtkFile << "LOOKUP_TABLE default" << std::endl;
	for (int i = 0; i < vecField.size(); i++)
	{
		vtkFile << vecField[i] << std::endl;
	}
	// write segs
	// vtkFile << "LINES " << nseg << " " << nseg * 3 << std::endl;
	// for (i = 1; i <= nseg; i++)
	//{
	//	Ng_GetSegment_2D(mesh, i, nodes, &matnum);
	//	std::cout << i << " matnum:\t" << matnum << " " << nodes[0] - 1 << " " << nodes[1] - 1 << std::endl;
	//}
	vtkFile.close();
}

int _tmain(int argc, _TCHAR *argv[])
{
	/************************************************************************/
	/* prepare data for mesh generation                                     */
	/************************************************************************/
	// 多边形外边界，通过指定四条边来定义
	CBoundShape shape;
	shape.boundShape = POLYGON;
	shape.nNumSegs = 4;
	Segment seg1;
	seg1.p1 = Point(-1500, -1500);
	seg1.p2 = Point(1500, -1500);
	Segment seg2;
	seg2.p1 = Point(1500, -1500);
	seg2.p2 = Point(1500, 1500);
	Segment seg3;
	seg3.p1 = Point(1500, 1500);
	seg3.p2 = Point(-1500, 1500);
	Segment seg4;
	seg4.p1 = Point(-1500, 1500);
	seg4.p2 = Point(-1500, -1500);

	shape.vecSegments.push_back(seg1);
	shape.vecSegments.push_back(seg2);
	shape.vecSegments.push_back(seg3);
	shape.vecSegments.push_back(seg4);

	// shape.cCenter = Point(0, 0);
	// shape.dRadius = 1000;

	CBoundWell w1, w2;
	// 第一口井为直井（默认为直井）
	w1.cCenter = Point(0, 0); // 井中心坐标(0,0)
	w1.dRadius = 0.09144;			// 井半径0.09144

	// 第二口井为垂直裂缝井
	w2.wellType = VERTICAL_FRACTURED;
	w2.dHf = 20;							 // 裂缝半长20米
	w2.dTheata = 30;					 // 裂缝角度30°
	w2.cCenter = Point(50, 0); // 井的中心位置（50，0）
	w2.dRadius = 0.1;					 // 井的半径0.1

	// 第三口井为多段压裂水平井
	CBoundWell w3;
	w3.wellType = MULTI_FRACED_HORIZONTAL;
	w3.cCenter = Point(142, -226); // 中心点位置
	w3.dRadius = 0.1;							 // 半径
	w3.dBeta = 20;								 // 水平井角度
	w3.dLength = 200;							 // 水平井长度
	w3.nFracNum = 10;							 // 裂缝数量
	w3.dHf = 20;									 // 裂缝半长

	std::vector<CBoundWell> wells;
	wells.push_back(w1);
	wells.push_back(w2);
	wells.push_back(w3);

	CBoundPolygon limit;
	limit.boundShape = POLYGON;
	limit.nNumSegs = 5;
	limit.vecSegments.push_back(Segment(Point(-516.2335347838829, -58.37577524138442), Point(-592.72401802758, -180.9965861656542)));
	limit.vecSegments.push_back(Segment(Point(-592.72401802758, -180.9965861656542), Point(-485.12090993338666, -280.2195626173259)));
	limit.vecSegments.push_back(Segment(Point(-485.12090993338666, -280.2195626173259), Point(-328.9162066043042, -229.04373034669038)));
	limit.vecSegments.push_back(Segment(Point(-328.9162066043042, -229.04373034669038), Point(-376.36137535544276, -90.33591813134171)));
	limit.vecSegments.push_back(Segment(Point(-376.36137535544276, -90.33591813134171), Point(-516.2335347838829, -58.37577524138442)));
	limit.bAnchor = true;
	std::vector<CBoundPolygon> limits;
	limits.push_back(limit);

	// 断层
	std::vector<Frac> faults;
	faults.push_back(Frac(Point(-267.3, 120), Point(316.7, 150)));
	faults.push_back(Frac(Point(316.7, 150), Point(316.7, -50)));

	// 裂缝
	std::vector<Frac> fracs;
	fracs.push_back(Frac(Point(-326.159, 460.059), Point(165.83631, 459.085)));
	fracs.push_back(Frac(Point(50, -350), Point(358.840, -289.811)));

	/************************************************************************/
	/* load netgen mesh generation library and generate mesh                */
	/************************************************************************/
	HMODULE hMod = LoadLibrary("meshgene.dll");
	if (NULL == hMod)
	{
		std::cout << "meshgene.dll加载失败！\n";
		return 1;
	}

	typedef Ng_Mesh *(*NgMeshGen2DByIn2d)(std::vector<CBoundWell> &wells, const std::vector<CBoundPolygon> limits, const std::vector<Frac> &faults, const std::vector<Frac> &fracs, const CBoundShape &shape);
	typedef Ng_Mesh *(*NgMeshGen2D)(std::string strFileIn2D);
	typedef void *(*Ng2VTK)(Ng_Mesh *mesh);
	// Get the address of the Ng_Mesh function
	NgMeshGen2DByIn2d meshgene1 = (NgMeshGen2DByIn2d)GetProcAddress(hMod, "NgMeshGen2DByIn2d");
	NgMeshGen2D meshgene2 = (NgMeshGen2D)GetProcAddress(hMod, "NgMeshGen2D");
	if (NULL == meshgene1)
	{
		FreeLibrary(hMod);
		std::cout << "meshgene文件加载函数地址获取失败！\n";
		return 1;
	}
	Ng_Mesh *mesh = meshgene1(wells, limits, faults, fracs, shape);

	// Ng_Mesh *mesh=meshgene2("netgen-1.in2d");
	Ng2VTK ng = (Ng2VTK)GetProcAddress(hMod, "Ng2VTK");
	ng(mesh);
	return 0;
}