gctl_mesh/lib/mesh/mesh.cpp

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  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.
 ******************************************************/

#include "mesh.h"

gctl::base_mesh::base_mesh()
{
    meshtype_ = UNDEFINED;
    meshdim_ = MESH_0D;
    meshname_ = "Untitled";
    meshinfo_ = "Undefined";
    node_num_ = ele_num_ = 0;
	initialized_ = false;
}

gctl::base_mesh::~base_mesh()
{
	clear();
}

void gctl::base_mesh::clear()
{
	meshtype_ = UNDEFINED;
    meshdim_ = MESH_0D;
    meshname_ = "Untitled";
    meshinfo_ = "Undefined";
    node_num_ = ele_num_ = 0;
	initialized_ = false;
    destroy_vector(datalist_);
	return;
}

bool gctl::base_mesh::initiated() const
{
	return initialized_;
}

bool gctl::base_mesh::saved(std::string datname) const
{
	if (datalist_.empty()) return false;
	else
	{
        for (size_t i = 0; i < datalist_.size(); i++)
        {
            if (datalist_[i].name_ == datname) return true;
        }

		return false;
	}
}

gctl::meshdata &gctl::base_mesh::get_data(std::string datname)
{
    return datalist_[data_index(datname)];
}

const gctl::meshdata &gctl::base_mesh::get_data(std::string datname) const
{
    return datalist_[data_index(datname)];
}

gctl::meshdata *gctl::base_mesh::get_data_ptr(std::string datname)
{
	return &datalist_[data_index(datname)];
}

void gctl::base_mesh::remove_data(std::string datname)
{
	int idx = data_index(datname);
    datalist_.erase(datalist_.begin() + idx);
	return;
}

void gctl::base_mesh::show_info(std::ostream &os) const
{
    if (meshtype_ == UNDEFINED) os << "Undefined | ";
	if (meshtype_ == REGULAR_MESH) os << "Regular | ";
	if (meshtype_ == LINEAR_MESH)  os << "Linear | ";
	if (meshtype_ == TRI_TET_MESH) os << "Unstructured | ";
	if (meshtype_ == REGULAR_MESH_SPH) os << "Regular (spherical) | ";
	if (meshtype_ == LINEAR_MESH_SPH)  os << "Linear (spherical) | ";
	if (meshtype_ == TRI_TET_MESH_SPH) os << "Unstructured (spherical) | ";
	if (meshtype_ == REGULAR_GRID) os << "Grid | ";

    if (meshdim_ == MESH_0D) os << "Unknown" << std::endl;
	else if (meshdim_ == MESH_2D) os << "2D" << std::endl;
	else if (meshdim_ == MESH_3D) os << "3D" << std::endl;

	os << "Name: " << meshname_ << std::endl;
	os << "Info: " << meshinfo_ << std::endl;

	show_mesh_dimension(os);

	for (size_t i = 0; i < datalist_.size(); i++)
    {
        datalist_[i].show_info(os);
    }
	return;
}

gctl::mesh_type_e gctl::base_mesh::get_meshtype() const
{
	check_initiated();
	return meshtype_;
}

gctl::mesh_dim_e gctl::base_mesh::get_meshdim() const
{
	check_initiated();
	return meshdim_;
}

int gctl::base_mesh::get_nodenum() const
{
	check_initiated();
	return node_num_;
}

int gctl::base_mesh::get_elenum() const
{
	check_initiated();
	return ele_num_;
}

int gctl::base_mesh::get_datanum() const
{
	check_initiated();
	return datalist_.size();
}

std::string gctl::base_mesh::get_meshname() const
{
	check_initiated();
	return meshname_;
}

void gctl::base_mesh::set_meshname(std::string in_name)
{
	check_initiated();
	meshname_ = in_name;
	return;
}

std::string gctl::base_mesh::get_meshinfo() const
{
	check_initiated();
	return meshinfo_;
}

void gctl::base_mesh::set_meshinfo(std::string in_info)
{
	check_initiated();
	meshinfo_ = in_info;
	return;
}

gctl::meshdata &gctl::base_mesh::add_data(mesh_data_type_e in_loctype, mesh_data_value_e in_valtype, 
    std::string name, double init_val, bool if_output, double nan_val)
{
    check_initiated();
    meshdata new_data(in_loctype, in_valtype, 0, name, if_output, nan_val);

    if (in_loctype == NodeData && in_valtype == Scalar) new_data.datval_.resize(node_num_, init_val);
    else if (in_loctype == ElemData && in_valtype == Scalar) new_data.datval_.resize(ele_num_, init_val);
    else if (in_loctype == NodeData && in_valtype == Vector) new_data.datval_.resize(3*node_num_, init_val);
    else if (in_loctype == ElemData && in_valtype == Vector) new_data.datval_.resize(3*ele_num_, init_val);
    else if (in_loctype == NodeData && in_valtype == Tensor) new_data.datval_.resize(9*node_num_, init_val);
    else if (in_loctype == ElemData && in_valtype == Tensor) new_data.datval_.resize(9*ele_num_, init_val);
    else throw std::runtime_error("[gctl::base_mesh] Fail to initialize the new data.");

    datalist_.push_back(new_data);
    return datalist_.back();
}

gctl::meshdata &gctl::base_mesh::add_data(mesh_data_type_e in_loctype, std::string name, 
    const array<double> &init_arr, bool if_output, double nan_val)
{
    check_initiated();
    meshdata new_data(in_loctype, Scalar, 0, name, if_output, nan_val);

    if (in_loctype == NodeData)
    {
        if (init_arr.size() != node_num_) throw std::runtime_error("[gctl::base_mesh] Invalid input data size.");
        new_data.datval_ = init_arr;
    }
    else if (in_loctype == ElemData)
    {
        if (init_arr.size() != ele_num_) throw std::runtime_error("[gctl::base_mesh] Invalid input data size.");
        new_data.datval_ = init_arr;
    }
    else throw std::runtime_error("[gctl::base_mesh] Invalid input data location.");

    datalist_.push_back(new_data);
    return datalist_.back();
}

gctl::meshdata &gctl::base_mesh::add_data(mesh_data_type_e in_loctype, std::string name, 
    const array<point3dc> &init_arr, bool if_output, double nan_val)
{
    check_initiated();
    meshdata new_data(in_loctype, Vector, 0, name, if_output, nan_val);

    if (in_loctype == NodeData)
    {
        if (init_arr.size() != node_num_) throw std::runtime_error("[gctl::base_mesh] Invalid input data size.");

        new_data.datval_.resize(3*node_num_);
        for (size_t i = 0; i < node_num_; i++)
        {
            new_data.datval_[3*i] = init_arr[i].x;
            new_data.datval_[3*i+1] = init_arr[i].y;
            new_data.datval_[3*i+2] = init_arr[i].z;
        }
    }
    else if (in_loctype == ElemData)
    {
        if (init_arr.size() != ele_num_) throw std::runtime_error("[gctl::base_mesh] Invalid input data size.");
        
        new_data.datval_.resize(3*ele_num_);
        for (size_t i = 0; i < ele_num_; i++)
        {
            new_data.datval_[3*i] = init_arr[i].x;
            new_data.datval_[3*i+1] = init_arr[i].y;
            new_data.datval_[3*i+2] = init_arr[i].z;
        }
    }
    else throw std::runtime_error("[gctl::base_mesh] Invalid input data location.");

    datalist_.push_back(new_data);
    return datalist_.back();
}

gctl::meshdata &gctl::base_mesh::add_data(mesh_data_type_e in_loctype, std::string name, 
    const array<tensor> &init_arr, bool if_output, double nan_val)
{
    check_initiated();
    meshdata new_data(in_loctype, Tensor, 0, name, if_output, nan_val);

    if (in_loctype == NodeData)
    {
        if (init_arr.size() != node_num_) throw std::runtime_error("[gctl::base_mesh] Invalid input data size.");

        new_data.datval_.resize(9*node_num_);
        for (size_t i = 0; i < node_num_; i++)
        {
            new_data.datval_[9*i] = init_arr[i].val[0][0];
            new_data.datval_[9*i+1] = init_arr[i].val[0][1];
            new_data.datval_[9*i+2] = init_arr[i].val[0][2];
            new_data.datval_[9*i+3] = init_arr[i].val[1][0];
            new_data.datval_[9*i+4] = init_arr[i].val[1][1];
            new_data.datval_[9*i+5] = init_arr[i].val[1][2];
            new_data.datval_[9*i+6] = init_arr[i].val[2][0];
            new_data.datval_[9*i+7] = init_arr[i].val[2][1];
            new_data.datval_[9*i+8] = init_arr[i].val[2][2];
        }
    }
    else if (in_loctype == ElemData)
    {
        if (init_arr.size() != ele_num_) throw std::runtime_error("[gctl::base_mesh] Invalid input data size.");
        
        new_data.datval_.resize(9*ele_num_);
        for (size_t i = 0; i < ele_num_; i++)
        {
            new_data.datval_[9*i] = init_arr[i].val[0][0];
            new_data.datval_[9*i+1] = init_arr[i].val[0][1];
            new_data.datval_[9*i+2] = init_arr[i].val[0][2];
            new_data.datval_[9*i+3] = init_arr[i].val[1][0];
            new_data.datval_[9*i+4] = init_arr[i].val[1][1];
            new_data.datval_[9*i+5] = init_arr[i].val[1][2];
            new_data.datval_[9*i+6] = init_arr[i].val[2][0];
            new_data.datval_[9*i+7] = init_arr[i].val[2][1];
            new_data.datval_[9*i+8] = init_arr[i].val[2][2];
        }
    }
    else throw std::runtime_error("[gctl::base_mesh] Invalid input data location.");

    datalist_.push_back(new_data);
    return datalist_.back();
}

void gctl::base_mesh::save_gmsh_withdata(std::string filename, output_type_e out_mode, index_packed_e packed)
{
	if (out_mode == OverWrite) save_gmsh(filename, packed);

	std::ofstream outfile;
	gctl::open_outfile(outfile, filename, ".msh", std::ios::out|std::ios::app);

	for (size_t i = 0; i < datalist_.size(); i++)
	{
		if (datalist_[i].loctype_ == NodeData)
		{
			if (datalist_[i].valtype_ == Scalar) gctl::save_gmsh_data(outfile, datalist_[i].name_, datalist_[i].datval_, gctl::NodeData, packed);
			else if (datalist_[i].valtype_ == Vector)
			{
				array<point3dc> vec_data(node_num_);
				for (size_t j = 0; j < node_num_; j++)
				{
					vec_data[j].x = datalist_[i].datval_[3*j];
					vec_data[j].y = datalist_[i].datval_[3*j+1];
					vec_data[j].z = datalist_[i].datval_[3*j+2];
				}

				gctl::save_gmsh_data(outfile, datalist_[i].name_, vec_data, gctl::NodeData, packed);
			}
			//else if (datalist_[i].valtype_ == Tensor)
			//{
			//	array<tensor> ten_data(node_num_);
			//	for (size_t j = 0; j < node_num_; j++)
			//	{
			//		ten_data[j].val[0][0] = datalist_[i].datval_[9*j];
			//		ten_data[j].val[0][1] = datalist_[i].datval_[9*j+1];
			//		ten_data[j].val[0][2] = datalist_[i].datval_[9*j+2];
			//		ten_data[j].val[1][0] = datalist_[i].datval_[9*j+3];
			//		ten_data[j].val[1][1] = datalist_[i].datval_[9*j+4];
			//		ten_data[j].val[1][2] = datalist_[i].datval_[9*j+5];
			//		ten_data[j].val[2][0] = datalist_[i].datval_[9*j+6];
			//		ten_data[j].val[2][1] = datalist_[i].datval_[9*j+7];
			//		ten_data[j].val[2][2] = datalist_[i].datval_[9*j+8];
			//	}
			//
			//	gctl::save_gmsh_data(outfile, datalist_[i].name_, ten_data, gctl::NodeData, packed);
			//}
			else throw std::runtime_error("[gctl::base_mesh] Invalid input data value type.");
		}
		else if (datalist_[i].loctype_ == ElemData)
		{
			if (datalist_[i].valtype_ == Scalar) gctl::save_gmsh_data(outfile, datalist_[i].name_, datalist_[i].datval_, gctl::ElemData, packed);
			else if (datalist_[i].valtype_ == Vector)
			{
				array<point3dc> vec_data(ele_num_);
				for (size_t j = 0; j < ele_num_; j++)
				{
					vec_data[j].x = datalist_[i].datval_[3*j];
					vec_data[j].y = datalist_[i].datval_[3*j+1];
					vec_data[j].z = datalist_[i].datval_[3*j+2];
				}

				gctl::save_gmsh_data(outfile, datalist_[i].name_, vec_data, gctl::ElemData, packed);
			}
			//else if (datalist_[i].valtype_ == Tensor)
			//{
			//	array<tensor> ten_data(ele_num_);
			//	for (size_t j = 0; j < ele_num_; j++)
			//	{
			//		ten_data[j].val[0][0] = datalist_[i].datval_[9*j];
			//		ten_data[j].val[0][1] = datalist_[i].datval_[9*j+1];
			//		ten_data[j].val[0][2] = datalist_[i].datval_[9*j+2];
			//		ten_data[j].val[1][0] = datalist_[i].datval_[9*j+3];
			//		ten_data[j].val[1][1] = datalist_[i].datval_[9*j+4];
			//		ten_data[j].val[1][2] = datalist_[i].datval_[9*j+5];
			//		ten_data[j].val[2][0] = datalist_[i].datval_[9*j+6];
			//		ten_data[j].val[2][1] = datalist_[i].datval_[9*j+7];
			//		ten_data[j].val[2][2] = datalist_[i].datval_[9*j+8];
			//	}
			//
			//	gctl::save_gmsh_data(outfile, datalist_[i].name_, ten_data, gctl::ElemData, packed);
			//}
			else throw std::runtime_error("[gctl::base_mesh] Invalid input data value type.");
		}
		else throw std::runtime_error("[gctl::base_mesh] Invalid input data location.");
	}

	outfile.close();
	return;
}

void gctl::base_mesh::save_gmsh_withdata(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
{
	if (out_mode == OverWrite) save_gmsh(filename, packed);

	std::ofstream outfile;
	gctl::open_outfile(outfile, filename, ".msh", std::ios::out|std::ios::app);

	const meshdata &data = get_data(datname);
	if (data.loctype_ == NodeData)
	{
		if (data.valtype_ == Scalar) gctl::save_gmsh_data(outfile, data.name_, data.datval_, gctl::NodeData, packed);
		else if (data.valtype_ == Vector)
		{
			array<point3dc> vec_data(node_num_);
			for (size_t j = 0; j < node_num_; j++)
			{
				vec_data[j].x = data.datval_[3*j];
				vec_data[j].y = data.datval_[3*j+1];
				vec_data[j].z = data.datval_[3*j+2];
			}

			gctl::save_gmsh_data(outfile, data.name_, vec_data, gctl::NodeData, packed);
		}
		//else if (data.valtype_ == Tensor)
		//{
		//	array<tensor> ten_data(node_num_);
		//	for (size_t j = 0; j < node_num_; j++)
		//	{
		//		ten_data[j].val[0][0] = data.datval_[9*j];
		//		ten_data[j].val[0][1] = data.datval_[9*j+1];
		//		ten_data[j].val[0][2] = data.datval_[9*j+2];
		//		ten_data[j].val[1][0] = data.datval_[9*j+3];
		//		ten_data[j].val[1][1] = data.datval_[9*j+4];
		//		ten_data[j].val[1][2] = data.datval_[9*j+5];
		//		ten_data[j].val[2][0] = data.datval_[9*j+6];
		//		ten_data[j].val[2][1] = data.datval_[9*j+7];
		//		ten_data[j].val[2][2] = data.datval_[9*j+8];
		//	}
		//
		//	gctl::save_gmsh_data(outfile, data.name_, ten_data, gctl::NodeData, packed);
		//}
		else throw std::runtime_error("[gctl::base_mesh] Invalid input data value type.");
	}
	else if (data.loctype_ == ElemData)
	{
		if (data.valtype_ == Scalar) gctl::save_gmsh_data(outfile, data.name_, data.datval_, gctl::ElemData, packed);
		else if (data.valtype_ == Vector)
		{
			array<point3dc> vec_data(ele_num_);
			for (size_t j = 0; j < ele_num_; j++)
			{
				vec_data[j].x = data.datval_[3*j];
				vec_data[j].y = data.datval_[3*j+1];
				vec_data[j].z = data.datval_[3*j+2];
			}

			gctl::save_gmsh_data(outfile, data.name_, vec_data, gctl::ElemData, packed);
		}
		//else if (data.valtype_ == Tensor)
		//{
		//	array<tensor> ten_data(ele_num_);
		//	for (size_t j = 0; j < ele_num_; j++)
		//	{
		//		ten_data[j].val[0][0] = data.datval_[9*j];
		//		ten_data[j].val[0][1] = data.datval_[9*j+1];
		//		ten_data[j].val[0][2] = data.datval_[9*j+2];
		//		ten_data[j].val[1][0] = data.datval_[9*j+3];
		//		ten_data[j].val[1][1] = data.datval_[9*j+4];
		//		ten_data[j].val[1][2] = data.datval_[9*j+5];
		//		ten_data[j].val[2][0] = data.datval_[9*j+6];
		//		ten_data[j].val[2][1] = data.datval_[9*j+7];
		//		ten_data[j].val[2][2] = data.datval_[9*j+8];
		//	}
		//
		//	gctl::save_gmsh_data(outfile, data.name_, ten_data, gctl::ElemData, packed);
		//}
		else throw std::runtime_error("[gctl::base_mesh] Invalid input data value type.");
	}
	else throw std::runtime_error("[gctl::base_mesh] Invalid input data location.");

	outfile.close();
	return;
}

/**
 * 以下是类的私有函数，可以简单一些
 */

void gctl::base_mesh::check_initiated(bool inverse) const
{
	if (inverse == true)
	{
		if (initialized_) throw std::runtime_error("[gctl::base_mesh] Mesh already initialized.");
		else return;
	}
	
	if (!initialized_) throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
    else return;
}

void gctl::base_mesh::init(mesh_type_e in_type, mesh_dim_e in_dim, std::string in_name, std::string in_info)
{
	if (in_name == "") throw std::runtime_error("[gctl::base_mesh] The input name is empty.");
	if (in_info == "") throw std::runtime_error("[gctl::base_mesh] The input info. is empty.");

	meshname_ = in_name;
	meshinfo_ = in_info;
	meshtype_ = in_type;
	meshdim_  = in_dim;
	return;
}

int gctl::base_mesh::data_index(std::string datname) const
{
    if (datalist_.empty()) throw std::runtime_error("[gctl::base_mesh] No data saved.");

    for (size_t i = 0; i < datalist_.size(); i++)
    {
        if (datalist_[i].name_ == datname) return i;
    }

	throw gctl::runtime_error("[gctl::base_mesh] No data found by the name: " + datname);
}

void gctl::base_mesh::load_headinfo(std::ifstream &infile, mesh_type_e expected_type, mesh_dim_e expected_dim)
{
	// 读入网格头信息
	infile.read((char*)&meshtype_, sizeof(int));
	infile.read((char*)&meshdim_, sizeof(int));
	if (meshdim_ != expected_dim || meshtype_ != expected_type)
	{
		infile.close();
		throw std::runtime_error("[gctl::base_mesh] Invalid input mesh type.");
	}

	int info_size;
	infile.read((char*)&info_size, sizeof(int));
	meshname_.resize(info_size);
	infile.read((char*)meshname_.c_str(), info_size);

	infile.read((char*)&info_size, sizeof(int));
	meshinfo_.resize(info_size);
	infile.read((char*)meshinfo_.c_str(), info_size);
	return;
}

void gctl::base_mesh::load_datablock(std::ifstream &infile)
{
	int in_num;
	meshdata new_data;

	infile.read((char*)&in_num, sizeof(int));
	for (int i = 0; i < in_num; i++)
	{
		new_data.clear();
		new_data.load_binary(infile);
		datalist_.push_back(new_data);
	}
	return;
}

void gctl::base_mesh::save_headinfo(std::ofstream &outfile)
{
	// 首先输出网格的类型和维度
	outfile.write((char*)&meshtype_, sizeof(int));
	outfile.write((char*)&meshdim_, sizeof(int));
	// 输出网格名称与信息
	int info_size = meshname_.size();
	outfile.write((char*)&info_size, sizeof(int));
	outfile.write((char*)meshname_.c_str(), info_size);
	info_size = meshinfo_.size();
	outfile.write((char*)&info_size, sizeof(int));
	outfile.write((char*)meshinfo_.c_str(), info_size);
	return;
}

void gctl::base_mesh::save_datablock(std::ofstream &outfile)
{	
	int num = 0;
	for (size_t i = 0; i < datalist_.size(); i++)
	{
		if (datalist_[i].output_ok_) num++;
	}

	outfile.write((char*) &num, sizeof(int));
	for (size_t i = 0; i < datalist_.size(); i++)
	{
		if (datalist_[i].output_ok_)
		{
			datalist_[i].save_binary(outfile); // 输出数据块
		}
	}
	return;
}