gctl_mesh/lib/mesh/meshio.h
2025-04-23 21:55:30 +08:00

734 lines
27 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* 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.
******************************************************/
#ifndef _GCTL_MESH_IO_H
#define _GCTL_MESH_IO_H
#include <unordered_set>
#include <map>
#include "gctl/math/refellipsoid.h"
#include "gctl/io/triangle_io.h"
#include "gctl/io/tetgen_io.h"
#include "gctl/io/gmsh_io.h"
namespace gctl
{
/**
* @brief 无效的索引缺省值。
*/
#define DEFAULT_INVALID_TAG -9999
/**
* @brief 网格单元体名称枚举类型。
*/
enum element_type_enum
{
NotSet,
_2NodeLine,
_3NodeTriangle,
_4NodeQuadrangle,
_4NodeTetrahedron,
_8NodeHexahedron,
_6NodePrism,
_5NodePyramid,
_3NodeSecondOrderLine,
_6NdoeSecondOrderLine,
_9NodeSecondOrderQuadrangle,
_10NodeSecondOrderTetrahedron,
_27NodeSecondOrderHexahedron,
_18NodeSecondOrderPrism,
_14NodeSecondOrderPyramid,
_1NodePoint,
_8NodeSecondOrderQuadrangle,
_20NdoeSecondOrderHexahedron,
_15NodeSecondOrderPrism,
_13NodeSecondOrderPyramid,
_9NodeThirdOrderIncompleteTriangle,
_10NdoeThirdOrderTriangle,
_12NodeFourthOrderIncompleteTriangle,
_15NodeFourthOrderTriangle,
_15NodeFifthOrderCompleteTriangle,
_21NodeFifthOrderCompleteTriangle,
_4NodeThirdOrderEdge,
_5NodeFourthOrderEdge,
_6NodeFifthOrderEdge,
_20NodeThirdOrderTetrahedron,
_35NodeFourthOrderTetrahedron,
_56NodeFifithOrderTetrahedron,
_64NodeThirdOrderHexahedron,
_125NodeFourthOrderHexahedron,
};
/**
* @brief 网格单元体标签类型枚举类型
*
*/
enum element_tag_enum
{
PhysicalTag, // 元素的物理分组标签
GeometryTag, // 元素的几何分组标签
PartitionTag, // 元素的剖分分组标签
NodeTag, // 顶点的标签(仅用于输出顶点标签数据)
};
/**
* @brief 网格单元体结构体
*
*/
struct meshio_element
{
bool enabled; // 单元体是否有效
int id; // 单元体编号
element_type_enum type; // 单元体类型
array<vertex3dc*> vert_ptrs; // 顶点指针数组
array<meshio_element*> neigh_ptrs; // 相邻单元体指针数组
meshio_element();
};
/**
* @brief 网格数据结构体
*
*/
struct meshio_data
{
bool enabled; // 数据体是否有效
mesh_data_type_e d_type; // 数据类型
array<std::string> str_tag; // 字符串类型的标签(默认为一个,即为数据的名称)
array<double> real_tag; // 实数类型的标签默认为一个等于0.0
array<int> int_tag; // 整数类型的标签(最少三个,最后一个为数据的长度)
array<void*> tar_ptrs; // 数据连接的对象指针数组 具体类型为vertex3dc*或meshio_element*
array<double> val; // 数据值(目前仅支持标量数据,后续再添加对矢量数据的支持)
meshio_data();
/**
* @brief 清空数组并重置变量。
*
*/
void clear();
/**
* @brief 检查数据体是否合规
*
*/
bool pass_check();
};
/**
* @brief 网格单元体分组结构体。
*
*/
struct meshio_element_group
{
bool enabled; // 组是否有效
element_type_enum type; // 组内单元体类型
int phys_group; // 物理分组标签
int geom_group; // 几何分组标签
int part_group; // 剖分分组标签
std::string name; // 组名
std::vector<meshio_element*> elem_ptrs; // 组内单元体指针数组
meshio_element_group();
/**
* @brief 将组内所有单元体设置为有效状态。
*
*/
void enable_elements();
/**
* @brief 将组内所有单元体设置为无效状态。
*
*/
void disable_elements();
};
/**
* @brief 网格读写类,这个类实现了多种数据格式的网格文件的读写操作。并具备简单的单元体操作功能。
*
*/
class mesh_io
{
public:
mesh_io();
virtual ~mesh_io();
/**
* @brief 重置(清空)网格数据。
*
*/
void reset();
/**
* @brief 输出网格数据信息至指定流。
*
* @param ss 指定流默认为clog
*/
void info(std::ostream &ss = std::clog);
/**
* @brief 按单元体类型编辑网格单元体组。
*
* @param swt 使能类型Enable或Disable
* @param e_type 单元体类型缺省值值NotSet表示对所有单元体组进行操作
*/
void edit_group(switch_type_e swt, element_type_enum e_type = NotSet);
/**
* @brief 按单元体组名称编辑网格单元体组。
*
* @param swt 使能类型Enable或Disable
* @param grp_name 单元体组名称。
*/
void edit_group(switch_type_e swt, std::string grp_name);
/**
* @brief 按单元体组标签编辑网格单元体组。
*
* @param swt 使能类型Enable或Disable
* @param tag_type 标签类型PhysicalTagGeometryTag或者PartitionTag
* @param tag 标签值。
*/
void edit_group(switch_type_e swt, element_tag_enum tag_type, int tag);
/**
* @brief 按单元体组标签编辑网格单元体组的名称。
*
* @param anchor_type 搜索的标签类型PhysicalTagGeometryTag或者PartitionTag
* @param anchor_group 搜索的标签值。
* @param new_name 单元体组的新名称。
*/
void edit_group(element_tag_enum anchor_type, int anchor_group, std::string new_name);
/**
* @brief 按单元体组标签搜索并编辑网格单元体组的标签。
*
* @param anchor_type 搜索的标签类型PhysicalTagGeometryTag或者PartitionTag
* @param anchor_group 搜索的标签值
* @param tar_type 更改的标签类型PhysicalTagGeometryTag或者PartitionTag
* @param tar_group 更改的标签值
*/
void edit_group(element_tag_enum anchor_type, int anchor_group, element_tag_enum tar_type, int tar_group);
/**
* @brief 返回指定类型与名称的标签值
*
* @param tag_type 查找的标签类型PhysicalTagGeometryTag或者PartitionTag
* @param group_name 查找的元素组名称。
* @return 标签值
*/
int group_tag(element_tag_enum tag_type, std::string group_name);
/**
* @brief 返回所有顶点数组的引用。
*
* @return 顶点数组的引用。
*/
const array<vertex3dc> &get_all_nodes();
/**
* @brief 返回所有单元体数组的引用。
*
* @return 单元体数组的引用。
*/
const array<meshio_element> &get_all_elems();
/**
* @brief 选择所有符合条件的单元体(注意只会统计有效的单元体组)。
*
* @param e_type 单元体类型缺省为NotSet即选择所有单元体类型
*/
void select_elements(element_type_enum e_type = NotSet);
/**
* @brief 选择所有符合条件的单元体(注意只会统计有效的单元体组)。
*
* @param tag_type 标签类型。
* @param tag 标签值。
*/
void select_elements(element_tag_enum tag_type, int tag);
/**
* @brief 选择所有符合条件的单元体(注意只会统计有效的单元体组)。
*
* @param phys_name 单元体组名称
*/
void select_elements(std::string phys_name);
/**
* @brief 选择所有符合条件的单元体(注意只会统计有效的单元体组)。
*
* @param datname 数据名称
* @param dtype 数据类型
*/
void select_elements(std::string dat_name, mesh_data_type_e dtype);
/**
* @brief 返回已选择的单元体所对应的顶点数量使用select_elements函数选择
*
* @return 顶点数量
*/
size_t selected_node_size();
/**
* @brief 返回已选择的单元体的数量使用select_elements函数选择
*
* @return 单元体数量
*/
size_t selected_element_size();
/**
* @brief 返回已选择的单元体所对应的顶点标签(默认的无效标签为-9999
*
* @return 整型数组的引用。
*/
const array<int> &selected_node_tags();
/**
* @brief 返回已选择的单元体所对应的元素标签(默认的无效标签为-9999
*
* @return 整型数组的引用。
*/
const array<int> &selected_element_tags(element_tag_enum tag_type);
/**
* @brief 获取gmsh格式分组表
*
* @param g_groups gmsh格式表
*/
void get_gmsh_physical_groups(std::vector<gmsh_physical_group> &g_groups);
/**
* @brief 将对应标签类型转换为网格数据类型(注意只会转换有效的顶点与单元体的标签)
*
* @param tag_type 标签类型。
*/
void convert_tags_to_data(element_tag_enum tag_type);
/**
* @brief 检查是否存在名为name的数据
*
* @param name 数据名称
* @param type 数据类型
*
* @return 存在则返回数据索引,不存在则返回-1。
*/
int if_saved_data(std::string name, mesh_data_type_e type);
/**
* @brief 获取数据对象的引用(注意此函数会直接返回数据的引用,注意可能会包含无效的顶点和单元体)
*
* @param name 数据名称
* @param type 数据类型
* @return 数据引用
*/
meshio_data &get_data(std::string name, mesh_data_type_e type);
/**
* @brief 获取数据对象的指针
*
* @param name 数据名称
* @param type 数据类型
* @return 数据指针
*/
meshio_data *get_data_ptr(std::string name, mesh_data_type_e type);
/**
* @brief 返回已选择顶点或单元体位置的数据值,会按照已选择顶点或者单元体索引排序。
*
* @param name 数据名称
* @param type 数据类型
*
* @return 输出的数据数组
*/
array<double> get_selected_data(std::string name, mesh_data_type_e type);
/**
* @brief 添加一个顶点数据对象。数据将依次添加到已选择的顶点位置。
*
* @param data 输入的数据数组,长度与已选择的顶点数据相同。
* @param name 新建的数据名称。
* @param dtype 数据类型。
* @param op 数据的添加方式OverWrite或Append
*/
void add_data(std::string name, const array<double> &data,
mesh_data_type_e dtype, output_type_e op = OverWrite);
/**
* @brief 编辑网格数据状态。
*
* @param swt 使能类型Enable或Disable
* @param dataname 数据名称缺省值为null表示对所有数据进行操作
*/
void edit_data(switch_type_e swt, std::string dataname = "null");
/**
* @brief 读入triangle软件输出的网格剖分文件。
*
* @param filename 文件名(.node和.ele文件必须在同一路径下.neigh文件不是必须的文件名不包含后缀名
* @param is_packed 输入文件的索引值是否从0开始。
*/
void read_triangle_ascii(std::string filename, index_packed_e is_packed = Packed);
/**
* @brief 读入tetgen软件输出的网格剖分文件。
*
* @param filename 文件名(.node和.ele文件必须在同一路径下.neigh文件不是必须的文件名不包含后缀名
* @param is_packed 输入文件的索引值是否从0开始。
*/
void read_tetgen_ascii(std::string filename, index_packed_e is_packed = Packed);
/**
* @brief 读入Gmsh软件输出的网格剖分文件只支持v2.2的ASCII文件
*
* @param filename 文件名
* @param is_packed 输入文件的索引值是否从0开始。
*/
void read_gmsh_v2_ascii(std::string filename, index_packed_e is_packed = NotPacked);
/**
* @brief 保存Gmsh软件格式的网格剖分文件只支持v2.2的ASCII文件
*
* @param filename 文件名
* @param is_packed 输入文件的索引值是否从0开始。
*/
void save_gmsh_v2_ascii(std::string filename, index_packed_e is_packed = NotPacked);
/**
* @brief 保存Paraview软件格式的网格剖分文件
*
* @param filename 文件名
*/
void save_vtk_legacy_ascii(std::string filename);
/**
* @brief 导出数据到点云文件
*
* @param filename 输出文件名
* @param dataname 数据名称
* @param out_coor 数据的坐标系Cartesian或Spherical
* @param refr 参考椭球的短半径out_coor为Cartesian时无效
* @param refR 参考椭球的长半径out_coor为Cartesian时无效
*/
void save_data_to_xyz(std::string filename, std::string dataname = "null", coordinate_system_e out_coor = Cartesian, double refr = GCTL_Earth_Radius, double refR = GCTL_Earth_Radius);
/**
* @brief 提取所选择的单元体到外部数组。
*
* @tparam T 单元体类型。
* @param elems 新生成的单元体数组的引用。
* @param nodes 新生成的顶点数组的引用。
*/
template <typename T>
void export_selected_to(array<T> &elems, array<vertex3dc> &nodes);
/**
* @brief 提取所选择的单元体到外部数组。
*
* @tparam T 单元体类型。
* @param elems 新生成的单元体数组的引用。
* @param nodes 新生成的顶点数组的引用。
* @param x_id 二维坐标系统的x索引缺省值为0
* @param y_id 二维坐标系统的y索引缺省值为1
*
* @note 二维坐标系统的索引表示取xyz坐标中的对应项分别赋值给二维坐标的xy值{0,1}即为xoy平面、
* {1,0}即为yox平面、{0,2}即为xoz平面、{1,2}即为yoz平面。
*/
template <typename T>
void export_selected_to(array<T> &elems, array<vertex2dc> &nodes, int x_id = 0, int y_id = 1);
/**
* @brief 导入外部单元体数组。
*
* @tparam T 单元体类型。
* @param elems 单元体数组。
* @param nodes 顶点数组。
*/
template <typename T>
void import_from(const array<T> &elems, const array<vertex3dc> &nodes);
/**
* @brief 导入外部单元体数组。
*
* @tparam T 单元体类型。
* @param elems 单元体数组。
* @param nodes 顶点数组。
* @param x_id 二维x坐标到三维坐标系统的索引缺省值为0
* @param y_id 三维y坐标到三维坐标系统的索引缺省值为1
*
* @note 二维(x,y)坐标映射到三维坐标(x,y,z){0,1}表示(x,y)映射至(x,y,0)、
* {1,0}表示(x,y)映射至(y,x,0)、{0,2}表示(x,y)映射至(x,0,z)、{1,2}表示(x,y)映射至(0,y,z)。
*/
template <typename T>
void import_from(const array<T> &elems, const array<vertex2dc> &nodes, int x_id = 0, int y_id = 1);
protected:
bool initialized_; // 类型是否已经初始化完成
// 有效的顶点、单元体和单元体组的数量
size_t valid_node_size_, valid_elem_size_, valid_group_size_;
array<vertex3dc> nodes_; // 网格顶点 当顶点索引为无效值时将不会被输出
array<meshio_element> elems_; // 网格元素
std::vector<meshio_data> datas_; // 网格数据
std::vector<meshio_element_group> groups_; // 网格单元体组
array<int> nodes_tag_; // 顶点标签
array<int> selected_node_tag_; // 被选择的顶点的标签
array<int> selected_elem_tag_; // 被选择的单元体的标签
std::vector<vertex3dc*> selected_nodes_; // 被选择的顶点
std::vector<meshio_element*> selected_elems_; // 被选择的单元体
std::vector<meshio_element_group*> selected_groups_; // 被选择的单元体组
element_type_enum elem_gmshcode2type_[94]; // gmsh的单元体类型数组 数组索引为gmsh的单元体类型码值
element_type_enum elem_vtkcode2type_[14]; // vtk的单元体类型数组 数组索引为vtk的单元体类型码值
std::map<element_type_enum, int> elem_type2gmshcode_; // 单元体类型到gmsh类型码值的映射
std::map<element_type_enum, int> elem_type2vtkcode_; // 单元体类型到vtk类型码值的映射
std::map<element_type_enum, int> elem_type2size_; // 单元体类型到单元体顶点数量的映射
std::map<element_type_enum, std::string> elem_type2name_; // 单元体类型到单元体名称的映射
std::string elem_name(element_type_enum e_type); // 获取单元体名称字符串
int elem_gmsh_code(element_type_enum e_type); // 获取单元体gmsh类型码值
int elem_vtk_code(element_type_enum e_type); // 获取单元体vtk类型码值
int elem_size(element_type_enum e_type); // 获取单元体顶点数量
element_type_enum elem_gmsh_type(int code); // 获取对应gmsh类型码的单元体类型
element_type_enum elem_vtk_type(int code); // 获取对应vtk类型码的单元体类型
void update_indexing(); // 更新索引(对网格的元素进行操作后需调用)
void sort_groups(); // 对单元体组进行梳理(对网格的元素进行操作后需调用)
element_type_enum match_type(const std::type_info &tinfo); // 根据type_info返回单元体类型
};
template <typename T>
void gctl::mesh_io::export_selected_to(array<T> &elems, array<vertex3dc> &nodes)
{
const std::type_info &tinfo = typeid(T);
element_type_enum oe_type = match_type(tinfo);
size_t n = selected_nodes_.size();
size_t s = selected_elems_.size();
std::map<int, int> node_map;
nodes.resize(n);
for (size_t i = 0; i < n; i++)
{
node_map[selected_nodes_[i]->id] = i; // 原网格顶点索引到新网格顶点索引的映射
nodes[i].id = i;
nodes[i].x = selected_nodes_[i]->x;
nodes[i].y = selected_nodes_[i]->y;
nodes[i].z = selected_nodes_[i]->z;
}
int vnum = elem_size(oe_type);
elems.resize(s);
for (size_t i = 0; i < s; i++)
{
if (selected_elems_[i]->type != oe_type)
throw gctl::invalid_argument("[gctl::mesh_io::export_to] The selected element is not " + elem_name(oe_type) + ".");
elems[i].id = i;
for (size_t v = 0; v < vnum; v++)
elems[i].vert[v] = nodes.get(node_map[selected_elems_[i]->vert_ptrs[v]->id]);
}
node_map.clear();
return;
}
template <typename T>
void gctl::mesh_io::export_selected_to(array<T> &elems, array<vertex2dc> &nodes, int x_id, int y_id)
{
const std::type_info &tinfo = typeid(T);
element_type_enum oe_type = match_type(tinfo);
size_t n = selected_nodes_.size();
size_t s = selected_elems_.size();
std::map<int, int> node_map;
nodes.resize(n);
double xyz_ref[3];
for (size_t i = 0; i < n; i++)
{
node_map[selected_nodes_[i]->id] = i; // 原网格顶点索引到新网格顶点索引的映射
nodes[i].id = i;
xyz_ref[0] = selected_nodes_[i]->x;
xyz_ref[1] = selected_nodes_[i]->y;
xyz_ref[2] = selected_nodes_[i]->z;
nodes[i].x = xyz_ref[x_id];
nodes[i].y = xyz_ref[y_id];
}
int vnum = elem_size(oe_type);
elems.resize(s);
for (size_t i = 0; i < s; i++)
{
if (selected_elems_[i]->type != oe_type)
throw gctl::invalid_argument("[gctl::mesh_io::export_to] The selected element is not " + elem_name(oe_type) + ".");
elems[i].id = i;
for (size_t v = 0; v < vnum; v++)
elems[i].vert[v] = nodes.get(node_map[selected_elems_[i]->vert_ptrs[v]->id]);
}
node_map.clear();
return;
}
template <typename T>
void gctl::mesh_io::import_from(const array<T> &elems, const array<vertex3dc> &nodes)
{
reset(); // 重置网格数据
const std::type_info &tinfo = typeid(T);
element_type_enum oe_type = match_type(tinfo);
int vnum = elem_size(oe_type);
valid_node_size_ = nodes.size();
valid_elem_size_ = elems.size();
nodes_.resize(valid_node_size_);
selected_nodes_.resize(valid_node_size_);
for (size_t i = 0; i < valid_node_size_; i++)
{
nodes_[i].id = i;
nodes_[i].x = nodes[i].x;
nodes_[i].y = nodes[i].y;
nodes_[i].z = nodes[i].z;
selected_nodes_[i] = nodes_.get(i);
}
elems_.resize(valid_elem_size_);
selected_elems_.resize(valid_elem_size_);
for (size_t i = 0; i < valid_elem_size_; i++)
{
elems_[i].id = i;
elems_[i].type = oe_type;
elems_[i].vert_ptrs.resize(vnum);
for (size_t v = 0; v < vnum; v++)
{
elems_[i].vert_ptrs[v] = nodes_.get(elems[i].vert[v]->id);
}
selected_elems_[i] = elems_.get(i);
}
// 所有单元体都属于同一个组
valid_group_size_ = 1;
groups_.resize(1);
groups_[0].enabled = true;
groups_[0].type = oe_type;
groups_[0].phys_group = 0; // 默认组别为0
groups_[0].geom_group = 0; // 默认组别为0
groups_[0].part_group = 3;
for (size_t i = 0; i < elems_.size(); i++)
{
groups_[0].elem_ptrs.push_back(elems_.get(i));
}
groups_[0].enable_elements();
initialized_ = true;
selected_groups_.resize(1);
selected_groups_[0] = &groups_[0];
return;
}
template <typename T>
void gctl::mesh_io::import_from(const array<T> &elems, const array<vertex2dc> &nodes, int x_id, int y_id)
{
reset(); // 重置网格数据
const std::type_info &tinfo = typeid(T);
element_type_enum oe_type = match_type(tinfo);
int vnum = elem_size(oe_type);
valid_node_size_ = nodes.size();
valid_elem_size_ = elems.size();
int xyz_ref[3];
nodes_.resize(valid_node_size_);
selected_nodes_.resize(valid_node_size_);
for (size_t i = 0; i < valid_node_size_; i++)
{
nodes_[i].id = i;
xyz_ref[0] = 0.0;
xyz_ref[1] = 0.0;
xyz_ref[2] = 0.0;
xyz_ref[x_id] = nodes[i].x;
xyz_ref[y_id] = nodes[i].y;
nodes_[i].x = xyz_ref[0];
nodes_[i].y = xyz_ref[1];
nodes_[i].z = xyz_ref[2];
selected_nodes_[i] = nodes_.get(i);
}
elems_.resize(valid_elem_size_);
selected_elems_.resize(valid_elem_size_);
for (size_t i = 0; i < valid_elem_size_; i++)
{
elems_[i].id = i;
elems_[i].type = oe_type;
elems_[i].vert_ptrs.resize(vnum);
for (size_t v = 0; v < vnum; v++)
{
elems_[i].vert_ptrs[v] = nodes_.get(elems[i].vert[v]->id);
}
selected_elems_[i] = elems_.get(i);
}
// 所有单元体都属于同一个组
valid_group_size_ = 1;
groups_.resize(1);
groups_[0].enabled = true;
groups_[0].type = oe_type;
groups_[0].phys_group = 0; // 默认组别为0
groups_[0].geom_group = 0; // 默认组别为0
groups_[0].part_group = 3;
for (size_t i = 0; i < elems_.size(); i++)
{
groups_[0].elem_ptrs.push_back(elems_.get(i));
}
groups_[0].enable_elements();
initialized_ = true;
selected_groups_.resize(1);
selected_groups_[0] = &groups_[0];
return;
}
};
#endif // _GCTL_MESH_IO_H