From c647ae37e97753a3dfa7ad6fc88bc61ffa1765b3 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Wed, 23 Apr 2025 17:05:28 +0800 Subject: [PATCH] move mesh_io to mesh --- lib/io.h | 1 - lib/io/gmsh_io.h | 2 +- lib/io/mesh_io.cpp | 2024 -------------------------------------------- lib/io/mesh_io.h | 734 ---------------- 4 files changed, 1 insertion(+), 2760 deletions(-) delete mode 100644 lib/io/mesh_io.cpp delete mode 100644 lib/io/mesh_io.h diff --git a/lib/io.h b/lib/io.h index de841e6..8ff08e6 100644 --- a/lib/io.h +++ b/lib/io.h @@ -41,7 +41,6 @@ #include "io/tetgen_io.h" #include "io/triangle_io.h" #include "io/gmsh_io.h" -#include "io/mesh_io.h" #include "io/text_io.h" #include "io/dsv_io.h" diff --git a/lib/io/gmsh_io.h b/lib/io/gmsh_io.h index 0a1c6ec..1ab9fb4 100644 --- a/lib/io/gmsh_io.h +++ b/lib/io/gmsh_io.h @@ -58,7 +58,7 @@ namespace gctl }; /** - * @brief Gmsh文件的IO类，主要为调用下面定义的全局函数 （此结构已被废弃，不再更新。可能在未来移除，建议使用mesh_io类型替代） + * @brief Gmsh文件的IO类，主要为调用下面定义的全局函数 */ class gmshio { diff --git a/lib/io/mesh_io.cpp b/lib/io/mesh_io.cpp deleted file mode 100644 index fcac492..0000000 --- a/lib/io/mesh_io.cpp +++ /dev/null @@ -1,2024 +0,0 @@ -/******************************************************** - * ██████╗ ██████╗████████╗██╗ - * ██╔════╝ ██╔════╝╚══██╔══╝██║ - * ██║ ███╗██║ ██║ ██║ - * ██║ ██║██║ ██║ ██║ - * ╚██████╔╝╚██████╗ ██║ ███████╗ - * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ - * 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 . - * - * 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_io.h" - -gctl::meshio_element::meshio_element() -{ - enabled = false; - id = DEFAULT_INVALID_TAG; - type = NotSet; -} - -gctl::meshio_data::meshio_data() -{ - enabled = false; - d_type = NodeData; -} - -void gctl::meshio_data::clear() -{ - enabled = false; - str_tag.clear(); - real_tag.clear(); - int_tag.clear(); - tar_ptrs.clear(); - val.clear(); - return; -} - -bool gctl::meshio_data::pass_check() -{ - // 检查是否同时连接了顶点和单元体 - if (tar_ptrs.empty()) return false; - if (int_tag[2] != val.size()) return false; - if (tar_ptrs.size() != val.size()) return false; - if (str_tag.empty() || real_tag.empty() || int_tag.size() < 3) return false; - return true; -} - -gctl::meshio_element_group::meshio_element_group() -{ - enabled = false; - type = NotSet; - name = "Untitled"; - phys_group = geom_group = part_group = DEFAULT_INVALID_TAG; -} - -void gctl::meshio_element_group::enable_elements() -{ - for (size_t e = 0; e < elem_ptrs.size(); e++) - { - elem_ptrs[e]->enabled = true; - } - return; -} - -void gctl::meshio_element_group::disable_elements() -{ - for (size_t e = 0; e < elem_ptrs.size(); e++) - { - elem_ptrs[e]->enabled = false; - } - return; -} - -gctl::mesh_io::mesh_io() -{ - valid_node_size_ = 0; - valid_elem_size_ = 0; - valid_group_size_ = 0; - initialized_ = false; - - elem_gmshcode2type_[0] = NotSet; - elem_gmshcode2type_[1] = _2NodeLine; - elem_gmshcode2type_[2] = _3NodeTriangle; - elem_gmshcode2type_[3] = _4NodeQuadrangle; - elem_gmshcode2type_[4] = _4NodeTetrahedron; - elem_gmshcode2type_[5] = _8NodeHexahedron; - elem_gmshcode2type_[6] = _6NodePrism; - elem_gmshcode2type_[7] = _5NodePyramid; - elem_gmshcode2type_[8] = _3NodeSecondOrderLine; - elem_gmshcode2type_[9] = _6NdoeSecondOrderLine; - elem_gmshcode2type_[10] = _9NodeSecondOrderQuadrangle; - elem_gmshcode2type_[11] = _10NodeSecondOrderTetrahedron; - elem_gmshcode2type_[12] = _27NodeSecondOrderHexahedron; - elem_gmshcode2type_[13] = _18NodeSecondOrderPrism; - elem_gmshcode2type_[14] = _14NodeSecondOrderPyramid; - elem_gmshcode2type_[15] = _1NodePoint; - elem_gmshcode2type_[16] = _8NodeSecondOrderQuadrangle; - elem_gmshcode2type_[17] = _20NdoeSecondOrderHexahedron; - elem_gmshcode2type_[18] = _15NodeSecondOrderPrism; - elem_gmshcode2type_[19] = _13NodeSecondOrderPyramid; - elem_gmshcode2type_[20] = _9NodeThirdOrderIncompleteTriangle; - elem_gmshcode2type_[21] = _10NdoeThirdOrderTriangle; - elem_gmshcode2type_[22] = _12NodeFourthOrderIncompleteTriangle; - elem_gmshcode2type_[23] = _15NodeFourthOrderTriangle; - elem_gmshcode2type_[24] = _15NodeFifthOrderCompleteTriangle; - elem_gmshcode2type_[25] = _21NodeFifthOrderCompleteTriangle; - elem_gmshcode2type_[26] = _4NodeThirdOrderEdge; - elem_gmshcode2type_[27] = _5NodeFourthOrderEdge; - elem_gmshcode2type_[28] = _6NodeFifthOrderEdge; - elem_gmshcode2type_[29] = _20NodeThirdOrderTetrahedron; - elem_gmshcode2type_[30] = _35NodeFourthOrderTetrahedron; - elem_gmshcode2type_[31] = _56NodeFifithOrderTetrahedron; - elem_gmshcode2type_[92] = _64NodeThirdOrderHexahedron; - elem_gmshcode2type_[93] = _125NodeFourthOrderHexahedron; - - elem_vtkcode2type_[3] = _2NodeLine; - elem_vtkcode2type_[5] = _3NodeTriangle; - elem_vtkcode2type_[9] = _4NodeQuadrangle; - elem_vtkcode2type_[10] = _4NodeTetrahedron; - elem_vtkcode2type_[12] = _8NodeHexahedron; - elem_vtkcode2type_[13] = _6NodePrism; - - elem_type2gmshcode_[NotSet] = 0; - elem_type2gmshcode_[_2NodeLine] = 1; - elem_type2gmshcode_[_3NodeTriangle] = 2; - elem_type2gmshcode_[_4NodeQuadrangle] = 3; - elem_type2gmshcode_[_4NodeTetrahedron] = 4; - elem_type2gmshcode_[_8NodeHexahedron] = 5; - elem_type2gmshcode_[_6NodePrism] = 6; - elem_type2gmshcode_[_5NodePyramid] = 7; - elem_type2gmshcode_[_3NodeSecondOrderLine] = 8; - elem_type2gmshcode_[_6NdoeSecondOrderLine] = 9; - elem_type2gmshcode_[_9NodeSecondOrderQuadrangle] = 10; - elem_type2gmshcode_[_10NodeSecondOrderTetrahedron] = 11; - elem_type2gmshcode_[_27NodeSecondOrderHexahedron] = 12; - elem_type2gmshcode_[_18NodeSecondOrderPrism] = 13; - elem_type2gmshcode_[_14NodeSecondOrderPyramid] = 14; - elem_type2gmshcode_[_1NodePoint] = 15; - elem_type2gmshcode_[_8NodeSecondOrderQuadrangle] = 16; - elem_type2gmshcode_[_20NdoeSecondOrderHexahedron] = 17; - elem_type2gmshcode_[_15NodeSecondOrderPrism] = 18; - elem_type2gmshcode_[_13NodeSecondOrderPyramid] = 19; - elem_type2gmshcode_[_9NodeThirdOrderIncompleteTriangle] = 20; - elem_type2gmshcode_[_10NdoeThirdOrderTriangle] = 21; - elem_type2gmshcode_[_12NodeFourthOrderIncompleteTriangle] = 22; - elem_type2gmshcode_[_15NodeFourthOrderTriangle] = 23; - elem_type2gmshcode_[_15NodeFifthOrderCompleteTriangle] = 24; - elem_type2gmshcode_[_21NodeFifthOrderCompleteTriangle] = 25; - elem_type2gmshcode_[_4NodeThirdOrderEdge] = 26; - elem_type2gmshcode_[_5NodeFourthOrderEdge] = 27; - elem_type2gmshcode_[_6NodeFifthOrderEdge] = 28; - elem_type2gmshcode_[_20NodeThirdOrderTetrahedron] = 29; - elem_type2gmshcode_[_35NodeFourthOrderTetrahedron] = 30; - elem_type2gmshcode_[_56NodeFifithOrderTetrahedron] = 31; - elem_type2gmshcode_[_64NodeThirdOrderHexahedron] = 92; - elem_type2gmshcode_[_125NodeFourthOrderHexahedron] = 93; - - elem_type2vtkcode_[_2NodeLine] = 3; - elem_type2vtkcode_[_3NodeTriangle] = 5; - elem_type2vtkcode_[_4NodeQuadrangle] = 9; - elem_type2vtkcode_[_4NodeTetrahedron] = 10; - elem_type2vtkcode_[_8NodeHexahedron] = 12; - elem_type2vtkcode_[_6NodePrism] = 13; - - elem_type2size_[_2NodeLine] = 2; - elem_type2size_[_3NodeTriangle] = 3; - elem_type2size_[_4NodeQuadrangle] = 4; - elem_type2size_[_4NodeTetrahedron] = 4; - elem_type2size_[_8NodeHexahedron] = 8; - elem_type2size_[_6NodePrism] = 6; - elem_type2size_[_5NodePyramid] = 5; - elem_type2size_[_3NodeSecondOrderLine] = 3; - elem_type2size_[_6NdoeSecondOrderLine] = 6; - elem_type2size_[_9NodeSecondOrderQuadrangle] = 9; - elem_type2size_[_10NodeSecondOrderTetrahedron] = 10; - elem_type2size_[_27NodeSecondOrderHexahedron] = 27; - elem_type2size_[_18NodeSecondOrderPrism] = 18; - elem_type2size_[_14NodeSecondOrderPyramid] = 14; - elem_type2size_[_1NodePoint] = 1; - elem_type2size_[_8NodeSecondOrderQuadrangle] = 8; - elem_type2size_[_20NdoeSecondOrderHexahedron] = 20; - elem_type2size_[_15NodeSecondOrderPrism] = 15; - elem_type2size_[_13NodeSecondOrderPyramid] = 13; - elem_type2size_[_9NodeThirdOrderIncompleteTriangle] = 9; - elem_type2size_[_10NdoeThirdOrderTriangle] = 10; - elem_type2size_[_12NodeFourthOrderIncompleteTriangle] = 12; - elem_type2size_[_15NodeFourthOrderTriangle] = 15; - elem_type2size_[_15NodeFifthOrderCompleteTriangle] = 15; - elem_type2size_[_21NodeFifthOrderCompleteTriangle] = 21; - elem_type2size_[_4NodeThirdOrderEdge] = 4; - elem_type2size_[_5NodeFourthOrderEdge] = 5; - elem_type2size_[_6NodeFifthOrderEdge] = 6; - elem_type2size_[_20NodeThirdOrderTetrahedron] = 20; - elem_type2size_[_35NodeFourthOrderTetrahedron] = 35; - elem_type2size_[_56NodeFifithOrderTetrahedron] = 56; - elem_type2size_[_64NodeThirdOrderHexahedron] = 64; - elem_type2size_[_125NodeFourthOrderHexahedron] = 125; - - elem_type2name_[_2NodeLine] = "2-node line"; - elem_type2name_[_3NodeTriangle] = "3-node triangle"; - elem_type2name_[_4NodeQuadrangle] = "4-node quadrangle"; - elem_type2name_[_4NodeTetrahedron] = "4-node tetrahedron"; - elem_type2name_[_8NodeHexahedron] = "8-node hexahedron"; - elem_type2name_[_6NodePrism] = "6-node prism"; - elem_type2name_[_5NodePyramid] = "5-node pyramid"; - elem_type2name_[_3NodeSecondOrderLine] = "3-node second order line"; - elem_type2name_[_6NdoeSecondOrderLine] = "6-ndoe second order line"; - elem_type2name_[_9NodeSecondOrderQuadrangle] = "9-node second order quadrangle"; - elem_type2name_[_10NodeSecondOrderTetrahedron] = "10-node second order tetrahedron"; - elem_type2name_[_27NodeSecondOrderHexahedron] = "27-node second order hexahedron"; - elem_type2name_[_18NodeSecondOrderPrism] = "18-node second order prism"; - elem_type2name_[_14NodeSecondOrderPyramid] = "14-node second order pyramid"; - elem_type2name_[_1NodePoint] = "1-node point"; - elem_type2name_[_8NodeSecondOrderQuadrangle] = "8-node second order quadrangle"; - elem_type2name_[_20NdoeSecondOrderHexahedron] = "20-ndoe second order hexahedron"; - elem_type2name_[_15NodeSecondOrderPrism] = "15-node second order prism"; - elem_type2name_[_13NodeSecondOrderPyramid] = "13-node second order pyramid"; - elem_type2name_[_9NodeThirdOrderIncompleteTriangle] = "9-node third order incomplete triangle"; - elem_type2name_[_10NdoeThirdOrderTriangle] = "10-ndoe third order triangle"; - elem_type2name_[_12NodeFourthOrderIncompleteTriangle] = "12-node fourth order incomplete triangle"; - elem_type2name_[_15NodeFourthOrderTriangle] = "15-node fourth order triangle"; - elem_type2name_[_15NodeFifthOrderCompleteTriangle] = "15-node fifth order complete triangle"; - elem_type2name_[_21NodeFifthOrderCompleteTriangle] = "21-node fifth order complete triangle"; - elem_type2name_[_4NodeThirdOrderEdge] = "4-node third order edge"; - elem_type2name_[_5NodeFourthOrderEdge] = "5-node fourth order edge"; - elem_type2name_[_6NodeFifthOrderEdge] = "6-node fifth order edge"; - elem_type2name_[_20NodeThirdOrderTetrahedron] = "20-node third order tetrahedron"; - elem_type2name_[_35NodeFourthOrderTetrahedron] = "35-node fourth order tetrahedron"; - elem_type2name_[_56NodeFifithOrderTetrahedron] = "56-node fifith order tetrahedron"; - elem_type2name_[_64NodeThirdOrderHexahedron] = "64-node third order hexahedron"; - elem_type2name_[_125NodeFourthOrderHexahedron] = "125-node fourth order hexahedron"; -} - -gctl::mesh_io::~mesh_io() -{ - reset(); -} - -void gctl::mesh_io::reset() -{ - if (!selected_nodes_.empty()) selected_nodes_.clear(); - if (!selected_elems_.empty()) selected_elems_.clear(); - if (!selected_groups_.empty()) selected_groups_.clear(); - if (!nodes_.empty()) nodes_.clear(); - if (!elems_.empty()) elems_.clear(); - if (!nodes_tag_.empty()) nodes_tag_.clear(); - - if (!groups_.empty()) - { - for (size_t i = 0; i < groups_.size(); i++) - { - groups_[i].elem_ptrs.clear(); - } - groups_.clear(); - } - - if (!datas_.empty()) - { - for (size_t i = 0; i < datas_.size(); i++) - { - datas_[i].clear(); - } - datas_.clear(); - } - - valid_node_size_ = 0; - valid_elem_size_ = 0; - valid_group_size_ = 0; - initialized_ = false; -} - -void gctl::mesh_io::info(std::ostream &ss) -{ - int inttemp,temp_count; - std::string strtemp; - std::stringstream stemp; - - ss << "number of nodes: " << valid_node_size_ << ", number of elements: " << valid_elem_size_ << std::endl; - ss << "name, element type, physical group, geometrical group, partition group, element#" << std::endl; - for (size_t g = 0; g < groups_.size(); g++) - { - if (!groups_[g].enabled) continue; - - ss << groups_[g].name << ", "; - ss << elem_name(groups_[g].type) << ", "; - ss << groups_[g].phys_group << ", "; - ss << groups_[g].geom_group << ", "; - ss << groups_[g].part_group << ", "; - ss << groups_[g].elem_ptrs.size() << std::endl; - } - - size_t d_size; - double min, max; - for (size_t d = 0; d < datas_.size(); d++) - { - if (!datas_[d].enabled) continue; - - d_size = 0; - min = GCTL_BDL_MAX; - max = GCTL_BDL_MIN; - - if (datas_[d].d_type == NodeData) - { - ss << "nodedata: \""; - for (size_t l = 0; l < datas_[d].val.size(); l++) - { - if (reinterpret_cast(datas_[d].tar_ptrs[l])->id != DEFAULT_INVALID_TAG) - { - min = GCTL_MIN(min, datas_[d].val[l]); - max = GCTL_MAX(max, datas_[d].val[l]); - d_size++; - } - } - } - else - { - ss << "elementdata: \""; - for (size_t l = 0; l < datas_[d].val.size(); l++) - { - if (reinterpret_cast(datas_[d].tar_ptrs[l])->enabled) - { - min = GCTL_MIN(min, datas_[d].val[l]); - max = GCTL_MAX(max, datas_[d].val[l]); - d_size++; - } - } - } - - if (d_size) ss << datas_[d].str_tag[0] << "\", data number: " << d_size << ", mini: " << min << ", maxi = " << max << std::endl; - else ss << datas_[d].str_tag[0] << "\", no valid value found.\n"; - } - - return; -} - -void gctl::mesh_io::edit_group(switch_type_e swt, element_type_enum e_type) -{ - for (size_t g = 0; g < groups_.size(); g++) - { - if ((e_type == NotSet || groups_[g].type == e_type) && swt == Enable) - { - groups_[g].enabled = true; - groups_[g].enable_elements(); - } - else if ((e_type == NotSet || groups_[g].type == e_type) && swt == Disable) - { - groups_[g].enabled = false; - groups_[g].disable_elements(); - } - } - - update_indexing(); - return; -} - -void gctl::mesh_io::edit_group(switch_type_e swt, std::string grp_name) -{ - for (size_t g = 0; g < groups_.size(); g++) - { - if (groups_[g].name == grp_name && swt == Enable) - { - groups_[g].enabled = true; - groups_[g].enable_elements(); - } - else if (groups_[g].name == grp_name && swt == Disable) - { - groups_[g].enabled = false; - groups_[g].disable_elements(); - } - } - - update_indexing(); - return; -} - -void gctl::mesh_io::edit_group(switch_type_e swt, element_tag_enum tag_type, int tag) -{ - for (size_t g = 0; g < groups_.size(); g++) - { - if (((tag_type == PhysicalTag && groups_[g].phys_group == tag) || - (tag_type == GeometryTag && groups_[g].geom_group == tag) || - (tag_type == PartitionTag && groups_[g].part_group == tag)) && - swt == Enable) - { - groups_[g].enabled = true; - groups_[g].enable_elements(); - } - else if (((tag_type == PhysicalTag && groups_[g].phys_group == tag) || - (tag_type == GeometryTag && groups_[g].geom_group == tag) || - (tag_type == PartitionTag && groups_[g].part_group == tag)) && - swt == Disable) - { - groups_[g].enabled = false; - groups_[g].disable_elements(); - } - } - - update_indexing(); - return; -} - -void gctl::mesh_io::edit_group(element_tag_enum anchor_type, int anchor_group, std::string new_name) -{ - for (size_t i = 0; i < groups_.size(); i++) - { - if ((anchor_type == PhysicalTag && groups_[i].phys_group == anchor_group) || - (anchor_type == GeometryTag && groups_[i].geom_group == anchor_group) || - (anchor_type == PartitionTag && groups_[i].part_group == anchor_group)) - { - groups_[i].name = new_name; - } - } - - // No need to update indexing here. - return; -} - -void gctl::mesh_io::edit_group(element_tag_enum anchor_type, int anchor_group, element_tag_enum tar_type, int tar_group) -{ - for (size_t i = 0; i < groups_.size(); i++) - { - if ((anchor_type == PhysicalTag && groups_[i].phys_group == anchor_group) || - (anchor_type == GeometryTag && groups_[i].geom_group == anchor_group) || - (anchor_type == PartitionTag && groups_[i].part_group == anchor_group)) - { - if (tar_type == PhysicalTag) groups_[i].phys_group = tar_group; - if (tar_type == GeometryTag) groups_[i].geom_group = tar_group; - if (tar_type == PartitionTag) groups_[i].part_group = tar_group; - } - } - - // 检查是否与现有分组重叠 有则合并 - sort_groups(); - return; -} - -int gctl::mesh_io::group_tag(element_tag_enum tag_type, std::string group_name) -{ - for (size_t i = 0; i < groups_.size(); i++) - { - if (tag_type == PhysicalTag && groups_[i].name == group_name) return groups_[i].phys_group; - if (tag_type == GeometryTag && groups_[i].name == group_name) return groups_[i].geom_group; - if (tag_type == PartitionTag && groups_[i].name == group_name) return groups_[i].part_group; - } - return DEFAULT_INVALID_TAG; -} - -const gctl::array &gctl::mesh_io::get_all_nodes() -{ - return nodes_; -} - -const gctl::array &gctl::mesh_io::get_all_elems() -{ - return elems_; -} - -void gctl::mesh_io::select_elements(element_type_enum e_type) -{ - if (!selected_nodes_.empty()) selected_nodes_.clear(); - if (!selected_elems_.empty()) selected_elems_.clear(); - if (!selected_groups_.empty()) selected_groups_.clear(); - - int vnum; - for (size_t i = 0; i < groups_.size(); i++) - { - if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled) - { - selected_groups_.push_back(&groups_[i]); - vnum = elem_size(groups_[i].type); - for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++) - { - selected_elems_.push_back(groups_[i].elem_ptrs[e]); - for (size_t v = 0; v < vnum; v++) - { - selected_nodes_.push_back(groups_[i].elem_ptrs[e]->vert_ptrs[v]); - } - } - } - } - - std::vector::iterator pos; - std::sort(selected_nodes_.begin(), selected_nodes_.end()); //排序 - pos = std::unique(selected_nodes_.begin(), selected_nodes_.end()); //获取重复序列开始的位置 - selected_nodes_.erase(pos, selected_nodes_.end()); //删除重复点 - return; -} - -void gctl::mesh_io::select_elements(element_tag_enum tag_type, int tag) -{ - if (!selected_nodes_.empty()) selected_nodes_.clear(); - if (!selected_elems_.empty()) selected_elems_.clear(); - if (!selected_groups_.empty()) selected_groups_.clear(); - - int vnum; - for (size_t i = 0; i < groups_.size(); i++) - { - if ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || - (tag_type == GeometryTag && groups_[i].geom_group == tag) || - (tag_type == PartitionTag && groups_[i].part_group == tag) && - groups_[i].enabled) - { - selected_groups_.push_back(&groups_[i]); - vnum = elem_size(groups_[i].type); - for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++) - { - selected_elems_.push_back(groups_[i].elem_ptrs[e]); - for (size_t v = 0; v < vnum; v++) - { - selected_nodes_.push_back(groups_[i].elem_ptrs[e]->vert_ptrs[v]); - } - } - } - } - - std::vector::iterator pos; - std::sort(selected_nodes_.begin(), selected_nodes_.end()); //排序 - pos = std::unique(selected_nodes_.begin(), selected_nodes_.end()); //获取重复序列开始的位置 - selected_nodes_.erase(pos, selected_nodes_.end()); //删除重复点 - return; -} - -void gctl::mesh_io::select_elements(std::string phys_name) -{ - if (!selected_nodes_.empty()) selected_nodes_.clear(); - if (!selected_elems_.empty()) selected_elems_.clear(); - if (!selected_groups_.empty()) selected_groups_.clear(); - - int vnum; - for (size_t i = 0; i < groups_.size(); i++) - { - if (groups_[i].enabled && groups_[i].name == phys_name) - { - selected_groups_.push_back(&groups_[i]); - vnum = elem_size(groups_[i].type); - for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++) - { - selected_elems_.push_back(groups_[i].elem_ptrs[e]); - for (size_t v = 0; v < vnum; v++) - { - selected_nodes_.push_back(groups_[i].elem_ptrs[e]->vert_ptrs[v]); - } - } - } - } - - std::vector::iterator pos; - std::sort(selected_nodes_.begin(), selected_nodes_.end()); //排序 - pos = std::unique(selected_nodes_.begin(), selected_nodes_.end()); //获取重复序列开始的位置 - selected_nodes_.erase(pos, selected_nodes_.end()); //删除重复点 - return; -} - -void gctl::mesh_io::select_elements(std::string dat_name, mesh_data_type_e dtype) -{ - if (!selected_nodes_.empty()) selected_nodes_.clear(); - if (!selected_elems_.empty()) selected_elems_.clear(); - if (!selected_groups_.empty()) selected_groups_.clear(); - - int d_id = if_saved_data(dat_name, dtype); - if (d_id == -1) throw std::runtime_error("[gctl::mesh_io::select_elements] Data not found."); - - if (dtype == ElemData) - { - int vnum; - meshio_element *e_ptr; - for (size_t i = 0; i < datas_[d_id].tar_ptrs.size(); i++) - { - e_ptr = reinterpret_cast(datas_[d_id].tar_ptrs[i]); - if (e_ptr->enabled) - { - selected_elems_.push_back(e_ptr); - vnum = elem_size(e_ptr->type); - for (size_t v = 0; v < vnum; v++) - { - selected_nodes_.push_back(e_ptr->vert_ptrs[v]); - } - } - } - - std::vector::iterator pos; - std::sort(selected_nodes_.begin(), selected_nodes_.end()); //排序 - pos = std::unique(selected_nodes_.begin(), selected_nodes_.end()); //获取重复序列开始的位置 - selected_nodes_.erase(pos, selected_nodes_.end()); //删除重复点 - } - else // NodeData - { - vertex3dc *v_ptr; - std::unordered_set nodes_set; - for (size_t i = 0; i < datas_[d_id].tar_ptrs.size(); i++) - { - v_ptr = reinterpret_cast(datas_[d_id].tar_ptrs[i]); - if (v_ptr->id != DEFAULT_INVALID_TAG) - { - selected_nodes_.push_back(v_ptr); - nodes_set.insert(v_ptr); - } - } - - int vnum; - bool skip_ele; - for (size_t i = 0; i < elems_.size(); i++) - { - if (elems_[i].enabled) - { - skip_ele = false; - vnum = elem_size(elems_[i].type); - for (size_t v = 0; v < vnum; v++) - { - if (nodes_set.count(elems_[i].vert_ptrs[v]) <= 0) - { - skip_ele = true; - break; - } - } - - if (skip_ele) break; - else selected_elems_.push_back(elems_.get(i)); - } - } - - } - return; -} - -size_t gctl::mesh_io::selected_node_size() -{ - return selected_nodes_.size(); -} - -size_t gctl::mesh_io::selected_element_size() -{ - return selected_elems_.size(); -} - -const gctl::array &gctl::mesh_io::selected_node_tags() -{ - selected_node_tag_.resize(selected_nodes_.size(), DEFAULT_INVALID_TAG); - - if (!nodes_tag_.empty()) - { - for (size_t i = 0; i < selected_node_tag_.size(); i++) - { - selected_node_tag_[i] = nodes_tag_[selected_nodes_[i]->id]; - } - } - return selected_node_tag_; -} - -const gctl::array &gctl::mesh_io::selected_element_tags(element_tag_enum tag_type) -{ - selected_elem_tag_.resize(selected_elems_.size()); - - size_t id = 0; - for (size_t i = 0; i < selected_groups_.size(); i++) - { - for (size_t e = 0; e < selected_groups_[i]->elem_ptrs.size(); e++) - { - if (tag_type == PhysicalTag) selected_elem_tag_[id] = selected_groups_[i]->phys_group; - else if (tag_type == GeometryTag) selected_elem_tag_[id] = selected_groups_[i]->geom_group; - else if (tag_type == PartitionTag) selected_elem_tag_[id] = selected_groups_[i]->part_group; - else selected_elem_tag_[id] = DEFAULT_INVALID_TAG; - id++; - } - } - return selected_elem_tag_; -} - -void gctl::mesh_io::get_gmsh_physical_groups(std::vector &g_groups) -{ - if (!g_groups.empty()) g_groups.clear(); - - gmsh_physical_group tmp_group; - bool not_found; - for (size_t i = 0; i < groups_.size(); i++) - { - if (g_groups.empty() && groups_[i].enabled) - { - tmp_group.name = groups_[i].name; - tmp_group.phys_tag = groups_[i].phys_group; - tmp_group.dim_tag = groups_[i].part_group; - g_groups.push_back(tmp_group); - } - else - { - not_found = true; - for (size_t g = 0; g < g_groups.size(); g++) - { - if (groups_[i].part_group == g_groups[g].dim_tag && - groups_[i].phys_group == g_groups[g].phys_tag) - { - not_found = false; - break; - } - } - - if (not_found && groups_[i].enabled) - { - tmp_group.name = groups_[i].name; - tmp_group.phys_tag = groups_[i].phys_group; - tmp_group.dim_tag = groups_[i].part_group; - g_groups.push_back(tmp_group); - } - } - } - return; -} - -void gctl::mesh_io::convert_tags_to_data(element_tag_enum tag_type) -{ - meshio_data tmp_data; - if (tag_type == NodeTag && (!nodes_tag_.empty())) - { - tmp_data.enabled = true; - tmp_data.d_type = NodeData; - tmp_data.str_tag.resize(1, "Node Tag"); - tmp_data.real_tag.resize(1, 0.0); - tmp_data.int_tag.resize(3, 0); - tmp_data.int_tag[1] = 1; - tmp_data.int_tag[2] = valid_node_size_; - - tmp_data.val.resize(valid_node_size_); - tmp_data.tar_ptrs.resize(valid_node_size_); - - size_t c = 0; - for (size_t i = 0; i < nodes_.size(); i++) - { - if (nodes_[i].id != DEFAULT_INVALID_TAG) - { - tmp_data.val[c] = (double) nodes_tag_[i]; - tmp_data.tar_ptrs[c] = nodes_.get(i); - c++; - } - } - - datas_.push_back(tmp_data); - } - else - { - tmp_data.enabled = true; - tmp_data.d_type = ElemData; - - if (tag_type == PhysicalTag) tmp_data.str_tag.resize(1, "Physical Tag"); - else if (tag_type == GeometryTag) tmp_data.str_tag.resize(1, "Geometry Tag"); - else if (tag_type == PartitionTag) tmp_data.str_tag.resize(1, "Partition Tag"); - - tmp_data.real_tag.resize(1, 0.0); - tmp_data.int_tag.resize(3, 0); - tmp_data.int_tag[1] = 1; - tmp_data.int_tag[2] = valid_elem_size_; - - tmp_data.val.resize(valid_elem_size_); - tmp_data.tar_ptrs.resize(valid_elem_size_); - - size_t c = 0; - for (size_t g = 0; g < groups_.size(); g++) - { - if (groups_[g].enabled) - { - for (size_t e = 0; e < groups_[g].elem_ptrs.size(); e++) - { - if (tag_type == PhysicalTag) tmp_data.val[c] = (double) groups_[g].phys_group; - if (tag_type == GeometryTag) tmp_data.val[c] = (double) groups_[g].geom_group; - if (tag_type == PartitionTag) tmp_data.val[c] = (double) groups_[g].part_group; - tmp_data.tar_ptrs[c] = groups_[g].elem_ptrs[e]; - c++; - } - } - } - - datas_.push_back(tmp_data); - } - return; -} - -int gctl::mesh_io::if_saved_data(std::string name, mesh_data_type_e type) -{ - for (size_t i = 0; i < datas_.size(); i++) - { - if (datas_[i].str_tag.front() == name && - datas_[i].d_type == type) return i; - } - return -1; -} - -gctl::meshio_data &gctl::mesh_io::get_data(std::string name, mesh_data_type_e type) -{ - int id = if_saved_data(name, type); - if (id == -1) throw std::runtime_error("[gctl::mesh_io::get_data] Data not found."); - return datas_[id]; -} - -gctl::meshio_data *gctl::mesh_io::get_data_ptr(std::string name, mesh_data_type_e type) -{ - int id = if_saved_data(name, type); - if (id == -1) throw std::runtime_error("[gctl::mesh_io::get_data] Data not found."); - return &datas_[id]; -} - -gctl::array gctl::mesh_io::get_selected_data(std::string name, mesh_data_type_e type) -{ - int id = if_saved_data(name, type); - if (id == -1) throw std::runtime_error("[gctl::mesh_io::get_data] Data not found."); - - array val; - if (datas_[id].d_type == NodeData) - { - val.resize(selected_nodes_.size(), GCTL_BDL_MAX); - vertex3dc *v_ptr; - for (size_t i = 0; i < selected_nodes_.size(); i++) - { - for (size_t j = 0; j < datas_[id].tar_ptrs.size(); j++) - { - v_ptr = reinterpret_cast(datas_[id].tar_ptrs[j]); - if (v_ptr == selected_nodes_[i]) - { - val[i] = datas_[id].val[j]; - break; - } - } - } - } - else // ElemData - { - val.resize(selected_elems_.size(), GCTL_BDL_MAX); - meshio_element *e_ptr; - for (size_t i = 0; i < selected_elems_.size(); i++) - { - for (size_t j = 0; j < datas_[id].tar_ptrs.size(); j++) - { - e_ptr = reinterpret_cast(datas_[id].tar_ptrs[j]); - if (e_ptr == selected_elems_[i]) - { - val[i] = datas_[id].val[j]; - break; - } - } - } - } - return val; -} - -void gctl::mesh_io::add_data(std::string name, const array &data, - mesh_data_type_e dtype, output_type_e op) -{ - size_t s; - if (dtype == NodeData) s = selected_nodes_.size(); - else s = selected_elems_.size(); - - if (data.size()!= s) throw std::runtime_error("[gctl::mesh_io::add_data] Incompatible data size to the selected nodes or elements."); - - int d_id; - if (dtype == NodeData) d_id = if_saved_data(name, NodeData); - else d_id = if_saved_data(name, ElemData); - - if (d_id != -1) // data exist - { - if (op == Append) - { - array tmp_ptrs(s, nullptr); - if (dtype == NodeData) - { - for (size_t i = 0; i < s; i++) - { - tmp_ptrs[i] = selected_nodes_[i]; - } - } - else - { - for (size_t i = 0; i < s; i++) - { - tmp_ptrs[i] = selected_elems_[i]; - } - } - - datas_[d_id].int_tag[2] += s; - datas_[d_id].tar_ptrs.concat(tmp_ptrs); - datas_[d_id].val.concat(data); - } - else // OverWrite - { - datas_[d_id].int_tag[2] = s; - datas_[d_id].tar_ptrs.resize(s, nullptr); - datas_[d_id].val.resize(s, 0.0); - - if (dtype == NodeData) - { - for (size_t i = 0; i < s; i++) - { - datas_[d_id].tar_ptrs[i] = selected_nodes_[i]; - datas_[d_id].val[i] = data[i]; - } - } - else - { - for (size_t i = 0; i < s; i++) - { - datas_[d_id].tar_ptrs[i] = selected_elems_[i]; - datas_[d_id].val[i] = data[i]; - } - } - } - return; - } - - meshio_data new_data; - new_data.enabled = true; - new_data.str_tag.resize(1, name); - new_data.real_tag.resize(1, 0.0); - new_data.int_tag.resize(3, 0); - new_data.int_tag[1] = 1; - new_data.int_tag[2] = s; - new_data.tar_ptrs.resize(s, nullptr); - new_data.val.resize(s); - - if (dtype == NodeData) - { - new_data.d_type = NodeData; - for (size_t i = 0; i < s; i++) - { - new_data.tar_ptrs[i] = selected_nodes_[i]; - new_data.val[i] = data[i]; - } - } - else - { - new_data.d_type = ElemData; - for (size_t i = 0; i < s; i++) - { - new_data.tar_ptrs[i] = selected_elems_[i]; - new_data.val[i] = data[i]; - } - } - - datas_.push_back(new_data); - return; -} - -void gctl::mesh_io::edit_data(switch_type_e swt, std::string dataname) -{ - if (dataname == "null" && swt == Enable) - { - for (size_t i = 0; i < datas_.size(); i++) - { - datas_[i].enabled = true; - } - } - else if (dataname == "null" && swt == Disable) - { - for (size_t i = 0; i < datas_.size(); i++) - { - datas_[i].enabled = false; - } - } - else - { - for (size_t i = 0; i < datas_.size(); i++) - { - if (datas_[i].str_tag.front() == dataname && swt == Enable) datas_[i].enabled = true; - if (datas_[i].str_tag.front() == dataname && swt == Disable) datas_[i].enabled = false; - } - } - return; -} - -void gctl::mesh_io::read_triangle_ascii(std::string filename, index_packed_e is_packed) -{ - if (initialized_ == true) - { - throw std::runtime_error("[gctl::mesh_io::read_triangle_ascii()] The mesh is already initialized. Run clear() first."); - } - - array node_2d; - array node_tag; - read_Triangle_node(filename, node_2d, is_packed, &node_tag); - if (!node_tag.empty()) nodes_tag_ = node_tag; - - array tri_2d; - matrix tri_attri; - read_Triangle_element(filename, tri_2d, node_2d, is_packed, &tri_attri); - - valid_node_size_ = node_2d.size(); - valid_elem_size_ = tri_2d.size(); - valid_group_size_ = 0; - - 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 = node_2d[i].x; - nodes_[i].y = node_2d[i].y; - nodes_[i].z = 0.0; - selected_nodes_[i] = nodes_.get(i); - } - - elems_.resize(valid_elem_size_); - selected_elems_.resize(valid_elem_size_); - matrix int_tag(valid_elem_size_, 3, DEFAULT_INVALID_TAG); - - for (size_t i = 0; i < valid_elem_size_; i++) - { - elems_[i].id = i; - elems_[i].neigh_ptrs.resize(3, nullptr); - elems_[i].type = _3NodeTriangle; - - int_tag[i][2] = 2; - - elems_[i].vert_ptrs.resize(3); - elems_[i].vert_ptrs[0] = nodes_.get(tri_2d[i].vert[0]->id); - elems_[i].vert_ptrs[1] = nodes_.get(tri_2d[i].vert[1]->id); - elems_[i].vert_ptrs[2] = nodes_.get(tri_2d[i].vert[2]->id); - selected_elems_[i] = elems_.get(i); - } - - if (!tri_attri.empty()) - { - for (size_t i = 0; i < valid_elem_size_; i++) - { - int_tag[i][1] = (int) tri_attri[i][0]; - } - } - - // 整理单元体组 - meshio_element_group tmp_group; - tmp_group.type = elems_[0].type; - tmp_group.phys_group = int_tag[0][0]; - tmp_group.geom_group = int_tag[0][1]; - tmp_group.part_group = int_tag[0][2]; - tmp_group.elem_ptrs.push_back(elems_.get(0)); - groups_.push_back(tmp_group); - - bool not_found; - for (size_t i = 1; i < elems_.size(); i++) - { - not_found = true; - for (size_t g = 0; g < groups_.size(); g++) - { - if (groups_[g].type == elems_[i].type && - groups_[g].phys_group == int_tag[i][0] && - groups_[g].geom_group == int_tag[i][1] && - groups_[g].part_group == int_tag[i][2]) - { - groups_[g].elem_ptrs.push_back(elems_.get(i)); - not_found = false; - break; - } - } - - if (not_found) - { - tmp_group.elem_ptrs.clear(); - tmp_group.type = elems_[i].type; - tmp_group.phys_group = int_tag[i][0]; - tmp_group.geom_group = int_tag[i][1]; - tmp_group.part_group = int_tag[i][2]; - tmp_group.elem_ptrs.push_back(elems_.get(i)); - groups_.push_back(tmp_group); - } - } - - std::string n_file = filename + ".neigh"; - if (!access(n_file.c_str(), R_OK)) - { - read_Triangle_neighbor(filename, tri_2d, is_packed); - - for (size_t i = 0; i < valid_elem_size_; i++) - { - for (size_t j = 0; j < 3; j++) - { - if (tri_2d[i].neigh[j] != nullptr) elems_[i].neigh_ptrs[j] = elems_.get(tri_2d[i].neigh[j]->id); - } - } - } - - // 最后将所有的组别和单元体设置为可用 - valid_group_size_ = groups_.size(); - selected_groups_.resize(valid_group_size_); - for (size_t g = 0; g < groups_.size(); g++) - { - groups_[g].enabled = true; - groups_[g].enable_elements(); - selected_groups_[g] = &groups_[g]; - } - initialized_ = true; - return; -} - -void gctl::mesh_io::read_tetgen_ascii(std::string filename, index_packed_e is_packed) -{ - if (initialized_ == true) - { - throw std::runtime_error("[gctl::mesh_io::read_tetgen_ascii()] The mesh is already initialized. Run clear() first."); - } - - array node_3d; - array node_tag; - read_Tetgen_node(filename, node_3d, is_packed, &node_tag); - if (!node_tag.empty()) nodes_tag_ = node_tag; - - array tet_3d; - array tet_tag; - read_Tetgen_element(filename, tet_3d, node_3d, is_packed, &tet_tag); - - std::string f_file = filename + ".face"; - array tri_3d; - array tri_tag; - if (!access(f_file.c_str(), R_OK)) - { - read_Tetgen_face(filename, tri_3d, node_3d, is_packed, &tri_tag); - } - - valid_node_size_ = node_3d.size(); - valid_elem_size_ = tet_3d.size() + tri_3d.size(); - valid_group_size_ = 0; - - 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 = node_3d[i].x; - nodes_[i].y = node_3d[i].y; - nodes_[i].z = node_3d[i].z; - selected_nodes_[i] = nodes_.get(i); - } - - elems_.resize(valid_elem_size_); - selected_elems_.resize(valid_elem_size_); - matrix int_tag(valid_elem_size_, 3, DEFAULT_INVALID_TAG); - - for (size_t i = 0; i < tet_3d.size(); i++) - { - elems_[i].id = i; - elems_[i].neigh_ptrs.resize(4, nullptr); - elems_[i].type = _4NodeTetrahedron; - - int_tag[i][2] = 3; // 剖分组别为三维 - - elems_[i].vert_ptrs.resize(4); - elems_[i].vert_ptrs[0] = nodes_.get(tet_3d[i].vert[0]->id); - elems_[i].vert_ptrs[1] = nodes_.get(tet_3d[i].vert[1]->id); - elems_[i].vert_ptrs[2] = nodes_.get(tet_3d[i].vert[2]->id); - elems_[i].vert_ptrs[3] = nodes_.get(tet_3d[i].vert[3]->id); - selected_elems_[i] = elems_.get(i); - } - - int offset = tet_3d.size(); - for (size_t i = 0; i < tri_3d.size(); i++) - {; - elems_[i + offset].id = i + offset; - elems_[i + offset].neigh_ptrs.resize(3, nullptr); - elems_[i + offset].type = _3NodeTriangle; - - int_tag[i + offset][2] = 2; - - elems_[i + offset].vert_ptrs.resize(3); - elems_[i + offset].vert_ptrs[0] = nodes_.get(tri_3d[i].vert[0]->id); - elems_[i + offset].vert_ptrs[1] = nodes_.get(tri_3d[i].vert[1]->id); - elems_[i + offset].vert_ptrs[2] = nodes_.get(tri_3d[i].vert[2]->id); - selected_elems_[i + offset] = elems_.get(i + offset); - } - - if (!tet_tag.empty()) - { - for (size_t i = 0; i < tet_3d.size(); i++) - { - int_tag[i][1] = tet_tag[i]; - } - } - - if (!tri_tag.empty()) - { - for (size_t i = 0; i < tri_3d.size(); i++) - { - int_tag[i + offset][1] = tri_tag[i]; - } - } - - // 整理单元体组 - meshio_element_group tmp_group; - tmp_group.type = elems_[0].type; - tmp_group.phys_group = int_tag[0][0]; - tmp_group.geom_group = int_tag[0][1]; - tmp_group.part_group = int_tag[0][2]; - tmp_group.elem_ptrs.push_back(elems_.get(0)); - groups_.push_back(tmp_group); - - bool not_found; - for (size_t i = 1; i < elems_.size(); i++) - { - not_found = true; - for (size_t g = 0; g < groups_.size(); g++) - { - if (groups_[g].type == elems_[i].type && - groups_[g].phys_group == int_tag[i][0] && - groups_[g].geom_group == int_tag[i][1] && - groups_[g].part_group == int_tag[i][2]) - { - groups_[g].elem_ptrs.push_back(elems_.get(i)); - not_found = false; - break; - } - } - - if (not_found) - { - tmp_group.elem_ptrs.clear(); - tmp_group.type = elems_[i].type; - tmp_group.phys_group = int_tag[i][0]; - tmp_group.geom_group = int_tag[i][1]; - tmp_group.part_group = int_tag[i][2]; - tmp_group.elem_ptrs.push_back(elems_.get(i)); - groups_.push_back(tmp_group); - } - } - - std::string n_file = filename + ".neigh"; - if (!access(n_file.c_str(), R_OK)) - { - read_Tetgen_neighbor(filename, tet_3d, is_packed); - - for (size_t i = 0; i < tet_3d.size(); i++) - { - for (size_t j = 0; j < 4; j++) - { - if (tet_3d[i].neigh[j] != nullptr) elems_[i].neigh_ptrs[j] = elems_.get(tet_3d[i].neigh[j]->id); - } - } - } - - // 最后将所有的组别和单元体设置为可用 - valid_group_size_ = groups_.size(); - selected_groups_.resize(valid_group_size_); - for (size_t g = 0; g < groups_.size(); g++) - { - groups_[g].enabled = true; - groups_[g].enable_elements(); - selected_groups_[g] = &groups_[g]; - } - initialized_ = true; - return; -} - -void gctl::mesh_io::read_gmsh_v2_ascii(std::string filename, index_packed_e is_packed) -{ - if (initialized_ == true) - { - throw std::runtime_error("[gctl::mesh_io::read_gmsh_v2_ascii()] The mesh is already initialized. Run clear() first."); - } - - std::ifstream infile; - open_infile(infile, filename, ".msh"); - - std::string tmp_str, tmp_str2; - std::stringstream tmp_ss; - - // 检查mesh文件版本 - while(getline(infile,tmp_str)) - { - if (tmp_str == "$MeshFormat") - { - getline(infile, tmp_str); - if (tmp_str != "2.2 0 8") - { - throw std::runtime_error("[gctl::mesh_io::read_gmsh_v2_ascii()] Unsupported mesh file format."); - } - break; - } - } - - //读入模型空间顶点集和物理分组 - size_t p_size = 0; - size_t n_size = 0; - array phys; - - while(getline(infile,tmp_str)) - { - if (tmp_str == "$PhysicalNames") - { - getline(infile, tmp_str); - gctl::str2ss(tmp_str, tmp_ss); - tmp_ss >> p_size; - - phys.resize(p_size); - for (size_t i = 0; i < p_size; i++) - { - getline(infile, tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> phys[i].dim_tag >> phys[i].phys_tag >> tmp_str2; - // 去掉两端的双引号 - replace_all(phys[i].name, tmp_str2, "\"", ""); - } - } - - if (tmp_str == "$Nodes") - { - getline(infile, tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> n_size; //第一个数为顶点个数 - nodes_.resize(n_size); - selected_nodes_.resize(n_size); - valid_node_size_ = n_size; - - for (int i = 0; i < n_size; i++) - { - getline(infile, tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> nodes_[i].id >> nodes_[i].x >> nodes_[i].y >> nodes_[i].z; - if (is_packed == NotPacked) nodes_[i].id -= 1; - - selected_nodes_[i] = nodes_.get(i); - } - break; - } - } - - // 读入模型空间元素集 - int i_size, type_code, attri_num, vt_idx; - array > file_itag; - while(getline(infile,tmp_str)) - { - if (tmp_str == "$Elements") - { - getline(infile,tmp_str); - gctl::str2ss(tmp_str, tmp_ss); - tmp_ss >> i_size; //第一个数为元素个数 - - valid_elem_size_ = i_size; - elems_.resize(i_size); - selected_elems_.resize(i_size); - file_itag.resize(valid_elem_size_); - - for (size_t i = 0; i < i_size; i++) - { - getline(infile,tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> elems_[i].id >> type_code >> attri_num; - if (is_packed == NotPacked) elems_[i].id -= 1; - - elems_[i].type = elem_gmsh_type(type_code); - elems_[i].vert_ptrs.resize(elem_size(elems_[i].type)); - // 这里暂时不考虑单元体的邻居关系 相关功能不应该由IO类型实现 - //elems_[i].neigh_ptrs.resize(elem_size(elems_[i].type), nullptr); - - // we get at least three tags. see below for tag types - // default tags will be assgined to DEFAULT_INVALID_TAG - file_itag[i].resize(GCTL_MAX(3, attri_num), DEFAULT_INVALID_TAG); - for (size_t a = 0; a < attri_num; a++) - { - tmp_ss >> file_itag[i][a]; - } - - for (size_t v = 0; v < elems_[i].vert_ptrs.size(); v++) - { - tmp_ss >> vt_idx; - if (is_packed == Packed) elems_[i].vert_ptrs[v] = nodes_.get(vt_idx); - else elems_[i].vert_ptrs[v] = nodes_.get(vt_idx - 1); - } - - selected_elems_[i] = elems_.get(i); - } - break; - } - } - - // 读入数据模块 - meshio_data tmp_data; - while(getline(infile,tmp_str)) - { - if (!infile.good()) break; - - if (tmp_str == "$NodeData" || tmp_str == "$ElementData") - { - tmp_data.clear(); - tmp_data.enabled = true; - - if (tmp_str == "$NodeData") tmp_data.d_type = NodeData; - if (tmp_str == "$ElementData") tmp_data.d_type = ElemData; - - getline(infile,tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> i_size; - tmp_data.str_tag.resize(i_size); - for (size_t i = 0; i < i_size; i++) - { - getline(infile, tmp_str); - replace_all(tmp_str2, tmp_str, "\"", ""); - tmp_data.str_tag[i] = tmp_str2; - } - - getline(infile,tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> i_size; - tmp_data.real_tag.resize(i_size); - for (size_t i = 0; i < i_size; i++) - { - getline(infile,tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> tmp_data.real_tag[i]; - } - - getline(infile,tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> i_size; - tmp_data.int_tag.resize(i_size); - for (size_t i = 0; i < i_size; i++) - { - getline(infile,tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> tmp_data.int_tag[i]; - } - - if (tmp_data.d_type == NodeData) - { - tmp_data.tar_ptrs.resize(tmp_data.int_tag.back()); - tmp_data.val.resize(tmp_data.int_tag.back()); - for (size_t i = 0; i < tmp_data.val.size(); i++) - { - getline(infile,tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> vt_idx >> tmp_data.val[i]; - if (is_packed == Packed) tmp_data.tar_ptrs[i] = nodes_.get(vt_idx); - else tmp_data.tar_ptrs[i] = nodes_.get(vt_idx - 1); - } - } - else - { - tmp_data.tar_ptrs.resize(tmp_data.int_tag.back()); - tmp_data.val.resize(tmp_data.int_tag.back()); - for (size_t i = 0; i < tmp_data.val.size(); i++) - { - getline(infile,tmp_str); - str2ss(tmp_str, tmp_ss); - tmp_ss >> vt_idx >> tmp_data.val[i]; - if (is_packed == Packed) tmp_data.tar_ptrs[i] = elems_.get(vt_idx); - else tmp_data.tar_ptrs[i] = elems_.get(vt_idx - 1); - } - } - - datas_.push_back(tmp_data); - continue; - } - } - infile.close(); - - // 整理单元体组 - meshio_element_group tmp_group; - tmp_group.type = elems_[0].type; - tmp_group.phys_group = file_itag[0][0]; - tmp_group.geom_group = file_itag[0][1]; - tmp_group.part_group = file_itag[0][2]; - tmp_group.elem_ptrs.push_back(elems_.get(0)); - groups_.push_back(tmp_group); - - // 元素类型与三个标记都一致的元素将被分到同一组 - bool not_found; - for (size_t i = 1; i < elems_.size(); i++) - { - not_found = true; - for (size_t g = 0; g < groups_.size(); g++) - { - if (groups_[g].type == elems_[i].type && - groups_[g].phys_group == file_itag[i][0] && - groups_[g].geom_group == file_itag[i][1] && - groups_[g].part_group == file_itag[i][2]) - { - groups_[g].elem_ptrs.push_back(elems_.get(i)); - not_found = false; - break; - } - } - - if (not_found) - { - tmp_group.elem_ptrs.clear(); - tmp_group.type = elems_[i].type; - tmp_group.phys_group = file_itag[i][0]; // 物理组 - tmp_group.geom_group = file_itag[i][1]; // 几何组 - tmp_group.part_group = file_itag[i][2]; // 剖分组（一般以元素的维度区分） - tmp_group.elem_ptrs.push_back(elems_.get(i)); - groups_.push_back(tmp_group); - } - } - - if (!phys.empty()) - { - // 遍历所有元素组 按物理组为标准为元素组命名 并以将元素维度赋值为剖分组 - for (size_t g = 0; g < groups_.size(); g++) - { - for (size_t p = 0; p < phys.size(); p++) - { - if (phys[p].phys_tag == groups_[g].phys_group) - { - groups_[g].name = phys[p].name; - groups_[g].part_group = phys[p].dim_tag; - break; - } - } - } - } - - // 最后将所有的组别和单元体设置为可用 - valid_group_size_ = groups_.size(); - selected_groups_.resize(valid_group_size_); - for (size_t g = 0; g < groups_.size(); g++) - { - groups_[g].enabled = true; - groups_[g].enable_elements(); - selected_groups_[g] = &groups_[g]; - } - initialized_ = true; - return; -} - -void gctl::mesh_io::save_gmsh_v2_ascii(std::string filename, index_packed_e is_packed) -{ - if (initialized_ == false) - { - throw std::runtime_error("The object is not initialized. From gctl::mesh_io::save_mesh(...)"); - } - - std::ofstream outfile; - open_outfile(outfile, filename, ".msh"); - - outfile << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n"; - - std::vector gmsh_groups; - get_gmsh_physical_groups(gmsh_groups); - - outfile << "$PhysicalNames\n" << gmsh_groups.size() << "\n"; - for (size_t i = 0; i < gmsh_groups.size(); i++) - { - if (gmsh_groups[i].phys_tag != DEFAULT_INVALID_TAG) outfile << gmsh_groups[i].dim_tag << " " << gmsh_groups[i].phys_tag << " \"" << gmsh_groups[i].name << "\"\n"; - else outfile << gmsh_groups[i].dim_tag << " 0 \"" << gmsh_groups[i].name << "\"\n"; - } - outfile << "$EndPhysicalNames\n"; - - outfile << "$Nodes\n" << valid_node_size_ << std::endl; - for (size_t i = 0; i < nodes_.size(); i++) - { - if (nodes_[i].id == DEFAULT_INVALID_TAG) continue; - - if (is_packed == NotPacked) outfile << nodes_[i].id + 1 << " "; - else outfile << nodes_[i].id << " "; - - outfile << std::setprecision(16) << nodes_[i].x << " " << nodes_[i].y << " " << nodes_[i].z << std::endl; - } - outfile<<"$EndNodes\n"; - - size_t a_size; - outfile << "$Elements\n" << valid_elem_size_ << std::endl; - for (size_t g = 0; g < groups_.size(); g++) - { - if (groups_[g].enabled) - { - for (size_t e = 0; e < groups_[g].elem_ptrs.size(); e++) - { - if (is_packed == NotPacked) outfile << groups_[g].elem_ptrs[e]->id + 1 << " " << elem_gmsh_code(groups_[g].type) << " "; - else outfile << groups_[g].elem_ptrs[e]->id << " " << elem_gmsh_code(groups_[g].type) << " "; - - a_size = 1; - if (groups_[g].geom_group != DEFAULT_INVALID_TAG) a_size++; - if (groups_[g].part_group != DEFAULT_INVALID_TAG) a_size++; - - outfile << a_size; - if (groups_[g].phys_group != DEFAULT_INVALID_TAG) outfile << " " << groups_[g].phys_group; - else outfile << " 0"; - - if (groups_[g].geom_group != DEFAULT_INVALID_TAG) outfile << " " << groups_[g].geom_group; - if (groups_[g].part_group != DEFAULT_INVALID_TAG) outfile << " " << groups_[g].part_group; - - for (size_t v = 0; v < groups_[g].elem_ptrs[e]->vert_ptrs.size(); v++) - { - if (is_packed == NotPacked) outfile << " " << groups_[g].elem_ptrs[e]->vert_ptrs[v]->id + 1; - else outfile << " " << groups_[g].elem_ptrs[e]->vert_ptrs[v]->id; - } - outfile << std::endl; - } - } - } - - outfile << "$EndElements\n"; - - if (!datas_.empty()) - { - size_t d_size; - for (size_t i = 0; i < datas_.size(); i++) - { - if (!datas_[i].enabled) continue; - - if (datas_[i].d_type == NodeData) - { - d_size = 0; - for (size_t n = 0; n < datas_[i].val.size(); n++) - { - if (reinterpret_cast(datas_[i].tar_ptrs[n])->id != DEFAULT_INVALID_TAG) d_size++; - } - - if (d_size) - { - outfile << "$NodeData\n"; - outfile << datas_[i].str_tag.size() << std::endl; - for (size_t a = 0; a < datas_[i].str_tag.size(); a++) - { - outfile << "\"" << datas_[i].str_tag[a] << "\"\n"; - } - - outfile << datas_[i].real_tag.size() << std::endl; - for (size_t a = 0; a < datas_[i].real_tag.size(); a++) - { - outfile << datas_[i].real_tag[a] << "\n"; - } - - outfile << datas_[i].int_tag.size() << std::endl; - for (size_t a = 0; a < datas_[i].int_tag.size() - 1; a++) - { - outfile << datas_[i].int_tag[a] << "\n"; - } - outfile << d_size << "\n"; - - int tmp_id; - for (size_t a = 0; a < datas_[i].val.size(); a++) - { - tmp_id = reinterpret_cast(datas_[i].tar_ptrs[a])->id; - if (tmp_id != DEFAULT_INVALID_TAG && is_packed == Packed) outfile << tmp_id << " " << std::setprecision(12) << datas_[i].val[a] << "\n"; - else if (tmp_id != DEFAULT_INVALID_TAG) outfile << tmp_id + 1 << " " << std::setprecision(12) << datas_[i].val[a] << "\n"; - } - outfile << "$EndNodeData\n"; - } - } - else - { - d_size = 0; - for (size_t n = 0; n < datas_[i].val.size(); n++) - { - if (reinterpret_cast(datas_[i].tar_ptrs[n])->enabled) d_size++; - } - - if (d_size) - { - outfile << "$ElementData\n"; - outfile << datas_[i].str_tag.size() << std::endl; - for (size_t a = 0; a < datas_[i].str_tag.size(); a++) - { - outfile << "\"" << datas_[i].str_tag[a] << "\"\n"; - } - - outfile << datas_[i].real_tag.size() << std::endl; - for (size_t a = 0; a < datas_[i].real_tag.size(); a++) - { - outfile << datas_[i].real_tag[a] << "\n"; - } - - outfile << datas_[i].int_tag.size() << std::endl; - for (size_t a = 0; a < datas_[i].int_tag.size() - 1; a++) - { - outfile << datas_[i].int_tag[a] << "\n"; - } - outfile << d_size << "\n"; - - int tmp_id; - bool d_ok; - for (size_t a = 0; a < datas_[i].val.size(); a++) - { - tmp_id = reinterpret_cast(datas_[i].tar_ptrs[a])->id; - d_ok = reinterpret_cast(datas_[i].tar_ptrs[a])->enabled; - if (d_ok && is_packed == Packed) outfile << tmp_id << " " << std::setprecision(12) << datas_[i].val[a] << "\n"; - else if (d_ok) outfile << tmp_id + 1 << " " << std::setprecision(12) << datas_[i].val[a] << "\n"; - } - outfile << "$EndElementData\n"; - } - } - } - } - - outfile.close(); - return; -} - -void gctl::mesh_io::save_vtk_legacy_ascii(std::string filename) -{ - if (initialized_ == false) - { - throw std::runtime_error("The object is not initialized. From gctl::mesh_io::export_vtk(...)"); - } - - std::ofstream outfile; - open_outfile(outfile, filename, ".vtk"); - - outfile << "# vtk DataFile Version 2.0\nGenerated by gctl::mesh_io\nASCII\nDATASET UNSTRUCTURED_GRID\nPOINTS " << valid_node_size_ <<" FLOAT\n"; - for (int i = 0; i < nodes_.size(); i++) - { - if (nodes_[i].id == DEFAULT_INVALID_TAG) continue; - - outfile << std::setprecision(16) << nodes_[i].x << " " << nodes_[i].y << " " << nodes_[i].z << std::endl; - } - - //计算一下CELLS的总长 注意类型名算一个 所以要加1 - size_t totalCellNum = 0; - for (int i = 0; i < elems_.size(); i++) - { - if (elems_[i].enabled) totalCellNum += elem_size(elems_[i].type) + 1; - } - outfile << "CELLS " << valid_elem_size_ << " " << totalCellNum << std::endl; - - for (int i = 0; i < elems_.size(); i++) - { - if (!elems_[i].enabled) continue; - - outfile << elem_size(elems_[i].type) << " "; - for (int j = 0; j < elems_[i].vert_ptrs.size(); j++) - { - outfile << elems_[i].vert_ptrs[j]->id << " "; - } - outfile << std::endl; - } - - outfile << "CELL_TYPES " << valid_elem_size_ << std::endl; - for (int i = 0; i < elems_.size(); i++) - { - if (elems_[i].enabled) outfile << elem_vtk_code(elems_[i].type) << "\n"; - } - - bool nd_head = true, ed_head = true; - size_t d_size; - std::string temp_dataName; - if (!datas_.empty()) - { - for (size_t i = 0; i < datas_.size(); i++) - { - if (!datas_[i].enabled) continue; - - if (datas_[i].d_type == NodeData) - { - d_size = 0; - for (size_t n = 0; n < datas_[i].val.size(); n++) - { - if (reinterpret_cast(datas_[i].tar_ptrs[n])->id != DEFAULT_INVALID_TAG) d_size++; - } - - if (d_size && d_size == valid_node_size_) - { - if (nd_head) {outfile<<"POINT_DATA "<< d_size << std::endl; nd_head = false; ed_head = true;} - - gctl::replace_all(temp_dataName, datas_[i].str_tag.front()," ","_"); - outfile << "SCALARS " << temp_dataName << " FLOAT\nLOOKUP_TABLE default\n"; - - for (size_t a = 0; a < datas_[i].val.size(); a++) - { - if (reinterpret_cast(datas_[i].tar_ptrs[a])->id != DEFAULT_INVALID_TAG) outfile << std::setprecision(16) << datas_[i].val[a] << "\n"; - } - } - } - else - { - d_size = 0; - for (size_t n = 0; n < datas_[i].val.size(); n++) - { - if (reinterpret_cast(datas_[i].tar_ptrs[n])->enabled) d_size++; - } - - if (d_size && d_size == valid_elem_size_) - { - if (ed_head) {outfile<<"CELL_DATA "<< d_size << std::endl; ed_head = false; nd_head = true;} - - gctl::replace_all(temp_dataName, datas_[i].str_tag.front()," ","_"); - outfile << "SCALARS " << temp_dataName << " FLOAT\nLOOKUP_TABLE default\n"; - - for (size_t a = 0; a < datas_[i].val.size(); a++) - { - if (reinterpret_cast(datas_[i].tar_ptrs[a])->enabled) outfile << std::setprecision(16) << datas_[i].val[a] << "\n"; - } - } - } - } - } - - outfile.close(); - return; -} - -void gctl::mesh_io::save_data_to_xyz(std::string filename, std::string dataname, coordinate_system_e out_coor, double refr, double refR) -{ - int tmp_size; - std::string tmp_name; - point3ds ps; - point3dc pc; - - std::ofstream ofile; - for (size_t i = 0; i < datas_.size(); i++) - { - if (dataname != "null" && (datas_[i].str_tag.front() != dataname || !datas_[i].enabled)) continue; - - gctl::replace_all(tmp_name, datas_[i].str_tag.front()," ","_"); - gctl::replace_all(tmp_name, tmp_name,"/",""); - open_outfile(ofile, filename + "_" + tmp_name, ".txt"); - - if (datas_[i].d_type == NodeData) - { - for (size_t n = 0; n < datas_[i].tar_ptrs.size(); n++) - { - vertex3dc* vptr = reinterpret_cast(datas_[i].tar_ptrs[n]); - if (vptr->id != DEFAULT_INVALID_TAG) - { - if (out_coor == Spherical) - { - ps = c2s(*vptr); - ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*GCTL_Pi/180.0); - - ofile << ps.lon << " " << ps.lat << " " << ps.rad << " " << datas_[i].val[n] << "\n"; - } - else ofile << pc.x << " " << pc.y << " " << pc.z << " " << datas_[i].val[n] << "\n"; - } - } - - } - else - { - for (size_t e = 0; e < datas_[i].tar_ptrs.size(); e++) - { - meshio_element* mptr = reinterpret_cast(datas_[i].tar_ptrs[e]); - if (mptr->enabled) - { - pc = point3dc(0.0, 0.0, 0.0); - tmp_size = mptr->vert_ptrs.size(); - for (size_t v = 0; v < mptr->vert_ptrs.size(); v++) - { - pc.x += mptr->vert_ptrs[v]->x; - pc.y += mptr->vert_ptrs[v]->y; - pc.z += mptr->vert_ptrs[v]->z; - } - pc.x /= tmp_size; - pc.y /= tmp_size; - pc.z /= tmp_size; - - if (out_coor == Spherical) - { - ps = c2s(pc); - ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*GCTL_Pi/180.0); - - ofile << ps.lon << " " << ps.lat << " " << ps.rad << " " << datas_[i].val[e] << "\n"; - } - else ofile << pc.x << " " << pc.y << " " << pc.z << " " << datas_[i].val[e] << "\n"; - } - } - } - ofile.close(); - } - return; -} - -std::string gctl::mesh_io::elem_name(element_type_enum e_type) -{ - return elem_type2name_[e_type]; -} - -int gctl::mesh_io::elem_gmsh_code(element_type_enum e_type) -{ - return elem_type2gmshcode_[e_type]; -} - -int gctl::mesh_io::elem_vtk_code(element_type_enum e_type) -{ - return elem_type2vtkcode_[e_type]; -} - -int gctl::mesh_io::elem_size(element_type_enum e_type) -{ - return elem_type2size_[e_type]; -} - -gctl::element_type_enum gctl::mesh_io::elem_gmsh_type(int code) -{ - return elem_gmshcode2type_[code]; -} - -gctl::element_type_enum gctl::mesh_io::elem_vtk_type(int code) -{ - return elem_vtkcode2type_[code]; -} - -void gctl::mesh_io::update_indexing() -{ - // 重置所有顶点的索引 - for (size_t n = 0; n < nodes_.size(); n++) - { - nodes_[n].id = DEFAULT_INVALID_TAG; - } - - // 确定有效单元个数并初始化对应的顶点索引为1 - int new_id = 0; - for (size_t e = 0; e < elems_.size(); e++) - { - if (elems_[e].enabled) - { - elems_[e].id = new_id; new_id++; - for (size_t n = 0; n < elems_[e].vert_ptrs.size(); n++) - { - elems_[e].vert_ptrs[n]->id = 1; - } - } - } - valid_elem_size_ = new_id; - - // 对所有有效顶点赋新的索引 同时确定有效顶点的个数 - new_id = 0; - for (size_t n = 0; n < nodes_.size(); n++) - { - if (nodes_[n].id == 1) - { - nodes_[n].id = new_id; new_id++; - } - } - valid_node_size_ = new_id; - - // 确定有效单元组的个数 - valid_group_size_ = 0; - for (size_t p = 0; p < groups_.size(); p++) - { - if (groups_[p].enabled) valid_group_size_++; - } - return; -} - -void gctl::mesh_io::sort_groups() -{ - // 拷贝到临时组 - std::vector tmp_groups = groups_; - - // 清空对象 - for (size_t i = 0; i < groups_.size(); i++) - { - groups_[i].elem_ptrs.clear(); - } - groups_.clear(); - - // 整理单元体组 - bool not_found; - meshio_element_group tmp_group; - for (size_t i = 0; i < tmp_groups.size(); i++) - { - tmp_group = tmp_groups[i]; - - not_found = true; - for (size_t g = 0; g < groups_.size(); g++) - { - // 注意名称不是分组的标准 多个组可以共用一个名称 - if (groups_[g].type == tmp_group.type && - groups_[g].phys_group == tmp_group.phys_group && - groups_[g].geom_group == tmp_group.geom_group && - groups_[g].part_group == tmp_group.part_group) - { - for (size_t e = 0; e < tmp_group.elem_ptrs.size(); e++) - { - groups_[g].elem_ptrs.push_back(tmp_group.elem_ptrs[e]); - } - - not_found = false; - } - - } - - if (not_found) - { - groups_.push_back(tmp_group); - } - - tmp_group.elem_ptrs.clear(); - } - - valid_group_size_ = 0; - for (size_t i = 0; i < groups_.size(); i++) - { - if (groups_[i].enabled) valid_group_size_++; - } - - destroy_vector(tmp_groups); - return; -} - -gctl::element_type_enum gctl::mesh_io::match_type(const std::type_info &tinfo) -{ - std::smatch ret; - std::regex pat_triangle("type_triangle"); - std::regex pat_tetrahedron("type_tetrahedron"); - std::regex pat_block("type_block"); - std::regex pat_edge("type_edge"); - std::regex pat_edge2d("type_edge2d"); - std::regex pat_node("type_node"); - std::regex pat_prism("type_prism"); - std::regex pat_rectangle2d("type_rectangle2d"); - std::regex pat_tricone("type_tricone"); - std::regex pat_triangle2d("type_triangle2d"); - - std::string t_str = tinfo.name(); - if (regex_search(t_str, ret, pat_triangle)) return _3NodeTriangle; - else if (regex_search(t_str, ret, pat_tetrahedron)) return _4NodeTetrahedron; - else if (regex_search(t_str, ret, pat_block)) return _8NodeHexahedron; - else if (regex_search(t_str, ret, pat_edge)) return _2NodeLine; - else if (regex_search(t_str, ret, pat_edge2d)) return _2NodeLine; - else if (regex_search(t_str, ret, pat_node)) return _1NodePoint; - else if (regex_search(t_str, ret, pat_prism)) return _6NodePrism; - else if (regex_search(t_str, ret, pat_rectangle2d)) return _4NodeQuadrangle; - else if (regex_search(t_str, ret, pat_tricone)) return _4NodeTetrahedron; - else if (regex_search(t_str, ret, pat_triangle2d)) return _3NodeTriangle; - else return NotSet; -} \ No newline at end of file diff --git a/lib/io/mesh_io.h b/lib/io/mesh_io.h deleted file mode 100644 index 2698039..0000000 --- a/lib/io/mesh_io.h +++ /dev/null @@ -1,734 +0,0 @@ -/******************************************************** - * ██████╗ ██████╗████████╗██╗ - * ██╔════╝ ██╔════╝╚══██╔══╝██║ - * ██║ ███╗██║ ██║ ██║ - * ██║ ██║██║ ██║ ██║ - * ╚██████╔╝╚██████╗ ██║ ███████╗ - * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ - * 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 . - * - * 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 -#include - -#include "../math/refellipsoid.h" -#include "triangle_io.h" -#include "tetgen_io.h" -#include "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 vert_ptrs; // 顶点指针数组 - array neigh_ptrs; // 相邻单元体指针数组 - - meshio_element(); - }; - - /** - * @brief 网格数据结构体 - * - */ - struct meshio_data - { - bool enabled; // 数据体是否有效 - mesh_data_type_e d_type; // 数据类型 - array str_tag; // 字符串类型的标签（默认为一个，即为数据的名称） - array real_tag; // 实数类型的标签（默认为一个，等于0.0） - array int_tag; // 整数类型的标签（最少三个，最后一个为数据的长度） - array tar_ptrs; // 数据连接的对象指针数组 具体类型为vertex3dc*或meshio_element* - array 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 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 标签类型（PhysicalTag，GeometryTag或者PartitionTag）。 - * @param tag 标签值。 - */ - void edit_group(switch_type_e swt, element_tag_enum tag_type, int tag); - - /** - * @brief 按单元体组标签编辑网格单元体组的名称。 - * - * @param anchor_type 搜索的标签类型（PhysicalTag，GeometryTag或者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 搜索的标签类型（PhysicalTag，GeometryTag或者PartitionTag）。 - * @param anchor_group 搜索的标签值 - * @param tar_type 更改的标签类型（PhysicalTag，GeometryTag或者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 查找的标签类型（PhysicalTag，GeometryTag或者PartitionTag）。 - * @param group_name 查找的元素组名称。 - * @return 标签值 - */ - int group_tag(element_tag_enum tag_type, std::string group_name); - - /** - * @brief 返回所有顶点数组的引用。 - * - * @return 顶点数组的引用。 - */ - const array &get_all_nodes(); - - /** - * @brief 返回所有单元体数组的引用。 - * - * @return 单元体数组的引用。 - */ - const array &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 &selected_node_tags(); - - /** - * @brief 返回已选择的单元体所对应的元素标签（默认的无效标签为-9999） - * - * @return 整型数组的引用。 - */ - const array &selected_element_tags(element_tag_enum tag_type); - - /** - * @brief 获取gmsh格式分组表 - * - * @param g_groups gmsh格式表 - */ - void get_gmsh_physical_groups(std::vector &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 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 &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 - void export_selected_to(array &elems, array &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 - void export_selected_to(array &elems, array &nodes, int x_id = 0, int y_id = 1); - - /** - * @brief 导入外部单元体数组。 - * - * @tparam T 单元体类型。 - * @param elems 单元体数组。 - * @param nodes 顶点数组。 - */ - template - void import_from(const array &elems, const array &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 - void import_from(const array &elems, const array &nodes, int x_id = 0, int y_id = 1); - - protected: - bool initialized_; // 类型是否已经初始化完成 - - // 有效的顶点、单元体和单元体组的数量 - size_t valid_node_size_, valid_elem_size_, valid_group_size_; - array nodes_; // 网格顶点 当顶点索引为无效值时将不会被输出 - array elems_; // 网格元素 - std::vector datas_; // 网格数据 - std::vector groups_; // 网格单元体组 - array nodes_tag_; // 顶点标签 - - array selected_node_tag_; // 被选择的顶点的标签 - array selected_elem_tag_; // 被选择的单元体的标签 - std::vector selected_nodes_; // 被选择的顶点 - std::vector selected_elems_; // 被选择的单元体 - std::vector selected_groups_; // 被选择的单元体组 - - element_type_enum elem_gmshcode2type_[94]; // gmsh的单元体类型数组 数组索引为gmsh的单元体类型码值 - element_type_enum elem_vtkcode2type_[14]; // vtk的单元体类型数组 数组索引为vtk的单元体类型码值 - std::map elem_type2gmshcode_; // 单元体类型到gmsh类型码值的映射 - std::map elem_type2vtkcode_; // 单元体类型到vtk类型码值的映射 - std::map elem_type2size_; // 单元体类型到单元体顶点数量的映射 - std::map 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 - void gctl::mesh_io::export_selected_to(array &elems, array &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 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 - void gctl::mesh_io::export_selected_to(array &elems, array &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 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 - void gctl::mesh_io::import_from(const array &elems, const array &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 - void gctl::mesh_io::import_from(const array &elems, const array &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 \ No newline at end of file