gctl_toolkits/gmsh2vtk/func_msh2vtk.h

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _FUNC_MSH2VTK_H
#define _FUNC_MSH2VTK_H
#include "iostream"
#include "fstream"
#include "string.h"
#include "sstream"
#include "stdio.h"
#include "stdlib.h"
#include "iomanip"
#include "cmath"
#include "vector"
#include "map"

#include "gctl/core/vector_wrapper.h"
#include "gctl/io.h"
#include "gctl/utility.h"

using namespace std;

typedef map<int,int> typesMap;

struct node
{
	int id;
	double x, y, z;
};
typedef vector<node> nodeArray;

struct element
{
	int id;
	int msh_type;
	int vtk_type;
	int msh_info_num;
	std::vector<int> msh_info;
	std::vector<int> nodes_index;
};
typedef vector<element> elementArray;

struct elementData
{
	std::vector<int> ids;
    string name;
	std::vector<double> val;
};
typedef vector<elementData> elementDataArray;

struct nodeData
{
    std::vector<int> ids;
    string name;
    std::vector<double> val;
};
typedef vector<nodeData> nodeDataArray;

class mshvtk
{
private:
	nodeArray nArray;
	elementArray eArray;
	nodeDataArray ndArray;
	elementDataArray edArray;
	typesMap msh2vtk_type;
	typesMap type_nodenum;
	typesMap::iterator itm;
	int nodenum;
	int elenum;
	bool nopacked;

	std::vector<std::string> outNames;
public:
	mshvtk(){
		msh2vtk_type[1] = 3;
		msh2vtk_type[2] = 5;
		msh2vtk_type[3] = 9;
		msh2vtk_type[4] = 10;
		msh2vtk_type[5] = 12;
		msh2vtk_type[6] = 13; //注意三棱柱顶点序列gmsh和vtk貌似有点区别 上下底面的三角形顶点排序刚好相反

		type_nodenum[1] = 2;
		type_nodenum[2] = 3;
		type_nodenum[3] = 4;
		type_nodenum[4] = 4;
		type_nodenum[5] = 8;
		type_nodenum[6] = 6;

		nopacked = false;
	}
	~mshvtk(){}
	int runtine(std::string parachar);
	int addname(std::string onename);
	int readmsh(std::string inname); //注意这里的文件名都是不带.msh后缀的
	int writevtk(std::string inname); //注意这里的文件名需要添加.vtk后缀
};

int mshvtk::runtine(std::string parachar)
{
	if(readmsh(parachar)) return -1;
	if(writevtk(parachar)) return -1;
	nArray.clear();
	ndArray.clear();
	eArray.clear();
	edArray.clear();
	return 0;
}

int mshvtk::addname(std::string onename)
{
	string s = "\""+onename+"\"";
	outNames.push_back(s);
	return 0;
}

int mshvtk::readmsh(std::string inname)
{
	int temp_int,ele_num;
	node temp_node;
	element temp_element;
	nodeData temp_node_data;
	elementData temp_element_data;
	string temp_str;
	stringstream temp_ss;

	ifstream infile;
	gctl::open_infile(infile, inname, ".msh");

	while(getline(infile,temp_str))
	{
		if (temp_str == "$Nodes")
		{
			getline(infile,temp_str);
			gctl::str2ss(temp_str, temp_ss);
			temp_ss >> nodenum;
			for (int i = 0; i < nodenum; i++)
			{
				getline(infile,temp_str);
				gctl::str2ss(temp_str, temp_ss);
				temp_ss >> temp_node.id >> temp_node.x >> temp_node.y >> temp_node.z;
				nArray.push_back(temp_node);
			}
		}
		else if (temp_str == "$Elements")
		{
			getline(infile,temp_str);
			gctl::str2ss(temp_str, temp_ss);
			temp_ss >> ele_num;
			for (int i = 0; i < ele_num; i++)
			{
				temp_element.msh_info.clear();
				temp_element.nodes_index.clear();
				getline(infile,temp_str);
				gctl::str2ss(temp_str, temp_ss);
				temp_ss >> temp_element.id >> temp_element.msh_type >> temp_element.msh_info_num;
				for (int j = 0; j < temp_element.msh_info_num; j++)
				{
					temp_ss >> temp_int;
					temp_element.msh_info.push_back(temp_int);
				}
				for (int j = 0; j < type_nodenum[temp_element.msh_type]; j++)
				{
					temp_ss >> temp_int;
					if (nArray[0].id == 1)
					{
						temp_element.nodes_index.push_back(temp_int-1);
						nopacked = true;
					}
					else temp_element.nodes_index.push_back(temp_int);
				}
				temp_element.vtk_type = msh2vtk_type[temp_element.msh_type];
				eArray.push_back(temp_element);
			}
		}
		else if (temp_str == "$ElementData")
		{
			temp_element_data.ids.clear();
			temp_element_data.val.clear();

			getline(infile,temp_str);
			gctl::str2ss(temp_str, temp_ss);
			temp_ss >> temp_int;
			for (int i = 0; i < temp_int; i++)
			{
				getline(infile,temp_str);
				temp_element_data.name = temp_str;
			}
			//跳过6行
			for (int i = 0; i < 6; i++)
				getline(infile,temp_str);
			gctl::str2ss(temp_str, temp_ss);
			temp_ss >> temp_int;

			temp_element_data.ids.resize(temp_int);
			temp_element_data.val.resize(temp_int);
			for (int i = 0; i < temp_int; i++)
			{
				getline(infile,temp_str);
				gctl::str2ss(temp_str, temp_ss);
				temp_ss >> temp_element_data.ids[i] >> temp_element_data.val[i];

				if (nopacked) temp_element_data.ids[i] -= 1;
			}
			edArray.push_back(temp_element_data);
		}
		else if (temp_str == "$NodeData")
		{
			temp_node_data.ids.clear();
			temp_node_data.val.clear();

			getline(infile,temp_str);
			gctl::str2ss(temp_str, temp_ss);
			temp_ss >> temp_int;
			for (int i = 0; i < temp_int; i++)
			{
				getline(infile,temp_str);
				temp_node_data.name = temp_str;
			}
			//跳过6行
			for (int i = 0; i < 6; i++)
				getline(infile,temp_str);
			gctl::str2ss(temp_str, temp_ss);
			temp_ss >> temp_int;

			temp_node_data.ids.resize(temp_int);
			temp_node_data.val.resize(temp_int);
			for (int i = 0; i < temp_int; i++)
			{
				getline(infile,temp_str);
				gctl::str2ss(temp_str, temp_ss);
				temp_ss >> temp_node_data.ids[i] >> temp_node_data.val[i];

				if (nopacked) temp_node_data.ids[i] -= 1;
			}
			ndArray.push_back(temp_node_data);
		}
	}
	infile.close();
	return 0;
}

int mshvtk::writevtk(std::string inname)
{
	int temp_int;
	string temp_dataName;
	int totalCellNum = 0;

	ofstream outfile;
	gctl::open_outfile(outfile,inname,".vtk");

	outfile<<"# vtk DataFile Version 2.0"<<endl<<"Generated by gmsh2vtk"<<endl<<"ASCII"<<endl<<"DATASET UNSTRUCTURED_GRID"<<endl<<"POINTS "<<nArray.size()<<" FLOAT"<<endl;
	for (int i = 0; i < nArray.size(); i++)
	{
		outfile << setprecision(16) << nArray.at(i).x << " " << nArray.at(i).y << " " << nArray.at(i).z << endl;
	}
	//计算一下CELLS的总长 注意类型名算一个 所以要加1
	for (int i = 0; i < eArray.size(); i++)
	{
		totalCellNum += type_nodenum[eArray.at(i).msh_type] + 1;
	}
	outfile<<"CELLS "<<eArray.size()<<" "<<totalCellNum<<endl;
	for (int i = 0; i < eArray.size(); i++)
	{
		outfile << type_nodenum[eArray.at(i).msh_type] << " ";
		for (int j = 0; j < eArray.at(i).nodes_index.size(); j++)
		{
			outfile << eArray.at(i).nodes_index.at(j) << " ";
		}
		outfile << endl;
	}
	outfile << "CELL_TYPES " << eArray.size() << endl;
	for (int i = 0; i < eArray.size(); i++)
	{
		outfile << eArray.at(i).vtk_type << "\n";
	}
	//outfile << endl;

	bool node_notout = true;
	if (!ndArray.empty())
	{
		for (int i = 0; i < ndArray.size(); i++)
		{
			for (int n = 0; n < outNames.size(); n++)
			{
				if (outNames[n] == ndArray[i].name)
				{
					if (node_notout)
					{
						outfile<<"POINT_DATA "<< nArray.size() <<endl;
						node_notout = false;
					}

					gctl::replace_all(temp_dataName, ndArray[i].name," ","_");	

					// If ndArray[i].val.size() is not nArray.size(). Padding with 1e+30
					if (ndArray[i].val.size() == nArray.size())
					{
						//outfile<<"POINT_DATA "<< ndArray[i].val.size() <<endl;
						outfile<<"SCALARS " << temp_dataName << " FLOAT"<<endl<<"LOOKUP_TABLE default"<<endl;
		
						for (int j = 0; j < ndArray[i].val.size(); j++)
						{
							outfile << setprecision(16) << ndArray[i].val[j] << endl;
						}
					}
					else
					{
						//outfile<<"POINT_DATA "<< nArray.size() <<endl;
						outfile<<"SCALARS " << temp_dataName << " FLOAT"<<endl<<"LOOKUP_TABLE default"<<endl;
			
						gctl::_1d_array n_val(nArray.size(), 1e+30);
						for (size_t j = 0; j < ndArray[i].val.size(); j++)
						{
							n_val[ndArray[i].ids[j]] = ndArray[i].val[j];
						}

						for (auto val: n_val)
						{
							outfile << setprecision(16) << val << endl;
						}
					}
				}
			}
		}
	}

	bool elem_notout = true;
	if (!edArray.empty())
	{
		for (int i = 0; i < edArray.size(); i++)
		{
			for (int n = 0; n < outNames.size(); n++)
			{
				if (outNames[n] == edArray[i].name)
				{
					if (elem_notout)
					{
						outfile<<"CELL_DATA "<< eArray.size() <<endl;
						elem_notout = false;
					}

					gctl::replace_all(temp_dataName, edArray[i].name," ","_");

					if (edArray[i].val.size() == eArray.size())
					{
						//outfile<<"CELL_DATA "<< edArray[i].val.size() <<endl;
						outfile<<"SCALARS " << temp_dataName << " FLOAT"<<endl<<"LOOKUP_TABLE default"<<endl;
						for (int j = 0; j < edArray[i].val.size(); j++)
						{
							outfile<<setprecision(16)<<edArray[i].val[j]<<endl;
						}
					}
					else
					{
						//outfile<<"CELL_DATA "<< eArray.size() <<endl;
						outfile<<"SCALARS " << temp_dataName << " FLOAT"<<endl<<"LOOKUP_TABLE default"<<endl;

						gctl::_1d_array e_val(eArray.size(), 1e+30);
						for (size_t j = 0; j < edArray[i].val.size(); j++)
						{
							e_val[edArray[i].ids[j]] = edArray[i].val[j];
						}

						for (auto val: e_val)
						{
							outfile << setprecision(16) << val << endl;
						}
					}
				}
			}
		}
	}

	outfile.close();
	return 0;
}
#endif