XIAN-FEM-2026June/3D opticsfem-master/mesh/Mesh_Interface.cpp

#include"Mesh_Base.h"

#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>

#include "../nlohmann/json.hpp"
using json = nlohmann::json;

void Mesh_2D::GetMesh(std::string str,std::string jstr)
{
	std::ifstream fd;
	fd.open(str.c_str());

	if (fd.is_open())
	{
		std::string data;
		//vertex
		fd >> data;
		if (data == "NbrVertex")
			fd >> _mNbrVertex;
		fd >> data;
		if (data == "Vertex")
		{
			_mVertex = Eigen::MatrixXd(_mNbrVertex, 3);
			for (int i = 0; i < _mNbrVertex; i++)
			{
				fd >> _mVertex(i, 0);
				fd >> _mVertex(i, 1);
				_mVertex(i, 2) = 0.;
			}
		}
		//edge
		fd >> data;
		if (data == "NbrEdge")
			fd >> _mNbrEdge;
		fd >> data;
		if (data == "Edge")
		{
			_mEdge = Eigen::MatrixXi(_mNbrEdge, 2);
			for (int i = 0; i < _mNbrEdge; i++)
			{
				fd >> _mEdge(i, 0);
				fd >> _mEdge(i, 1);
			}
		}
		//tri
		fd >> data;
		if (data == "NbrTri")
			fd >> _mNbrTri;
		fd >> data;
		if (data == "Tri")
		{
			_mTri = Eigen::MatrixXi(_mNbrTri, 3);
			for (int i = 0; i < _mNbrTri; i++)
			{
				fd >> _mTri(i, 0);
				fd >> _mTri(i, 1);
				fd >> _mTri(i, 2);
			}
		}
		fd >> data;
		if (data == "EdgeOfTri")
		{
			_mEdgeOfTri = Eigen::MatrixXi(_mNbrTri, 3);
			for (int i = 0; i < _mNbrTri; i++)
			{
				fd >> _mEdgeOfTri(i, 0);
				fd >> _mEdgeOfTri(i, 1);
				fd >> _mEdgeOfTri(i, 2);
			}
		}
		fd >> data;
		if (data == "DomainOfTri")
		{
			_mDomainOfTri = Eigen::VectorXi(_mNbrTri);
			for (int i = 0; i < _mNbrTri; i++)
			{
				fd >> _mDomainOfTri(i);
				_mDomainOfTri(i) = _mDomainOfTri(i) + 1;
			}
		}
		//edges
		fd >> data;
		if (data == "NbrEdges")
			fd >> _mNbrEdges;
		fd >> data;
		if (data == "CoonOfEdges")
		{
			_mCoonOfEdges = Eigen::MatrixXi(_mNbrEdges, 3);
			for (int i = 0; i < _mNbrEdges; i++)
			{
				fd >> _mCoonOfEdges(i, 0);
				fd >> _mCoonOfEdges(i, 1);
				_mCoonOfEdges(i, 0) = _mCoonOfEdges(i, 0);
				_mCoonOfEdges(i, 1) = _mCoonOfEdges(i, 1);
			}
		}
		fd >> data;
		if (data == "DomainOfEdges")
		{
			_mDomainOfEdges = Eigen::VectorXi(_mNbrEdges);
			for (int i = 0; i < _mNbrEdges; i++)
			{
				fd >> _mDomainOfEdges(i);
				_mDomainOfEdges(i) = _mDomainOfEdges(i) + 1;
			}
		}
		fd >> data;
		json js = nlohmann::json::parse(jstr);
		int nbrBoundary = js.at("NbrBoundary");
		if (data == "normal")
		{
			_mNormOfEdges = Eigen::MatrixXd(nbrBoundary, 3);
			for (int i = 0; i < nbrBoundary; i++)
			{
				fd >> _mNormOfEdges(i, 0);
				fd >> _mNormOfEdges(i, 1);
				fd >> _mNormOfEdges(i, 2);
			}
		}
		fd >> data;
		if (data == "NbrCopyOfEdges")
			fd >> _mNbrCopyOfEdges;
		fd >> data;
		if (data == "CopyOfEdges")
		{
			_mCopyOfEdges = Eigen::MatrixXi::Zero(_mNbrCopyOfEdges, 2);
			for (int i = 0; i < _mNbrCopyOfEdges; i++)
			{
				fd >> _mCopyOfEdges(i, 0);
				fd >> _mCopyOfEdges(i, 1);
				_mCopyOfEdges(i, 0) = _mCopyOfEdges(i, 0);
				_mCopyOfEdges(i, 1) = _mCopyOfEdges(i, 1);
			}
		}
		fd >> data;
		if (data == "NbrCopyOfVertexs")
			fd >> _mNbrCopyOfVertexs;
		fd >> data;
		if (data == "CopyOfVertexs")
		{
			_mCopyOfVertex = Eigen::MatrixXi::Zero(_mNbrCopyOfVertexs, 2);
			for (int i = 0; i < _mNbrCopyOfVertexs; i++)
			{
				fd >> _mCopyOfVertex(i, 0);
				fd >> _mCopyOfVertex(i, 1);
				_mCopyOfVertex(i, 0) = _mCopyOfVertex(i, 0);
				_mCopyOfVertex(i, 1) = _mCopyOfVertex(i, 1);
			}
		}
	}
	else
	{
		std::cout << "mesh read fail" << std::endl;
		//<2F><>??<3F><>
	}

}

void Mesh_3D::GetMesh(std::string str, std::string jstr)
{
	std::ifstream fd;
	fd.open(str.c_str());

	if (!fd.is_open())
	{
		std::cout << "mesh read fail: " << str << std::endl;
		return;
	}

	std::string data;
	fd >> data;
	if (data == "NbrVertex")
		fd >> _mNbrVertex;
	fd >> data;
	if (data == "Vertex")
	{
		_mVertex = Eigen::MatrixXd(_mNbrVertex, 3);
		for (int i = 0; i < _mNbrVertex; i++)
		{
			fd >> _mVertex(i, 0);
			fd >> _mVertex(i, 1);
			fd >> _mVertex(i, 2);
		}
	}

	fd >> data;
	if (data == "NbrTet")
		fd >> _mNbrTet;
	fd >> data;
	if (data == "Tet")
	{
		_mTet = Eigen::MatrixXi(_mNbrTet, 4);
		for (int i = 0; i < _mNbrTet; i++)
		{
			for (int j = 0; j < 4; j++)
			{
				fd >> _mTet(i, j);
				_mTet(i, j) -= 1;
			}
		}
	}

	fd >> data;
	if (data == "DomainOfTet")
	{
		_mDomainOfTet = Eigen::VectorXi(_mNbrTet);
		for (int i = 0; i < _mNbrTet; i++)
			fd >> _mDomainOfTet(i);
	}

	fd >> data;
	if (data == "NbrEdge")
		fd >> _mNbrEdge;
	fd >> data;
	if (data == "Edge")
	{
		_mEdge = Eigen::MatrixXi(_mNbrEdge, 2);
		for (int i = 0; i < _mNbrEdge; i++)
		{
			fd >> _mEdge(i, 0);
			fd >> _mEdge(i, 1);
			_mEdge(i, 0) -= 1;
			_mEdge(i, 1) -= 1;
		}
	}

	fd >> data;
	if (data == "EdgeOfTet")
	{
		_mEdgeOfTet = Eigen::MatrixXi(_mNbrTet, 6);
		for (int i = 0; i < _mNbrTet; i++)
		{
			for (int j = 0; j < 6; j++)
			{
				fd >> _mEdgeOfTet(i, j);
				_mEdgeOfTet(i, j) -= 1;
			}
		}
	}

	fd >> data;
	if (data == "NbrTri")
		fd >> _mNbrTri;
	fd >> data;
	if (data == "Tri")
	{
		_mTri = Eigen::MatrixXi(_mNbrTri, 3);
		for (int i = 0; i < _mNbrTri; i++)
		{
			for (int j = 0; j < 3; j++)
			{
				fd >> _mTri(i, j);
				_mTri(i, j) -= 1;
			}
		}
	}

	fd >> data;
	if (data == "DomainOfTri")
	{
		_mDomainOfTri = Eigen::VectorXi(_mNbrTri);
		for (int i = 0; i < _mNbrTri; i++)
			fd >> _mDomainOfTri(i);
	}

	fd >> data;
	if (data == "ConnOfTri")
	{
		_mConnOfTri = Eigen::MatrixXi(_mNbrTri, 2);
		for (int i = 0; i < _mNbrTri; i++)
		{
			fd >> _mConnOfTri(i, 0);
			fd >> _mConnOfTri(i, 1);
			_mConnOfTri(i, 0) -= 1;
			_mConnOfTri(i, 1) -= 1;
		}
	}

	fd >> data;
	if (data == "NormOfFace")
	{
		int nbrNorm = 0;
		fd >> nbrNorm;
		int maxDomain = 0;
		std::vector<int> domainIds(nbrNorm);
		std::vector<Eigen::Vector3d> normals(nbrNorm);
		for (int i = 0; i < nbrNorm; i++)
		{
			int domainId = 0;
			fd >> domainId;
			domainIds[i] = domainId;
			normals[i] = Eigen::Vector3d::Zero();
			fd >> normals[i](0);
			fd >> normals[i](1);
			fd >> normals[i](2);
			if (domainId > maxDomain)
				maxDomain = domainId;
		}
		_mNormOfFace = Eigen::MatrixXd::Zero(maxDomain + 1, 3);
		for (int i = 0; i < nbrNorm; i++)
		{
			_mNormOfFace(domainIds[i], 0) = normals[i](0);
			_mNormOfFace(domainIds[i], 1) = normals[i](1);
			_mNormOfFace(domainIds[i], 2) = normals[i](2);
		}
		// getMesh.m defaults (SBCmesh.dat may omit domain 4 used by assembly_inc)
		const int kDefaultNormRows = 15;
		if (_mNormOfFace.rows() < kDefaultNormRows)
		{
			Eigen::MatrixXd expanded = Eigen::MatrixXd::Zero(kDefaultNormRows, 3);
			expanded.topRows(_mNormOfFace.rows()) = _mNormOfFace;
			_mNormOfFace = expanded;
		}
		auto setDefaultNorm = [&](int domainId, double nx, double ny, double nz)
		{
			if (_mNormOfFace.row(domainId).norm() < 1e-30)
			{
				_mNormOfFace(domainId, 0) = nx;
				_mNormOfFace(domainId, 1) = ny;
				_mNormOfFace(domainId, 2) = nz;
			}
		};
		setDefaultNorm(1, -1.0, 0.0, 0.0);
		setDefaultNorm(2, 0.0, -1.0, 0.0);
		setDefaultNorm(3, 0.0, 0.0, -1.0);
		setDefaultNorm(4, 0.0, 0.0, 1.0);
		setDefaultNorm(5, 0.0, 1.0, 0.0);
		setDefaultNorm(14, 1.0, 0.0, 0.0);
	}

	fd.close();
}