/******************************************************** * ██████╗ ██████╗████████╗██╗ * ██╔════╝ ██╔════╝╚══██╔══╝██║ * ██║ ███╗██║ ██║ ██║ * ██║ ██║██║ ██║ ██║ * ╚██████╔╝╚██████╗ ██║ ███████╗ * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ * 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 (!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 = std::min(min, datas_[d].val[l]); max = std::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 = std::min(min, datas_[d].val[l]); max = std::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_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::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; } const gctl::array &gctl::mesh_io::get_node_tag() { return nodes_tag_; } int gctl::mesh_io::get_tag(element_tag_enum anchor_type, std::string anchor_name) { for (size_t i = 0; i < groups_.size(); i++) { if (anchor_type == PhysicalTag && groups_[i].name == anchor_name) return groups_[i].phys_group; if (anchor_type == GeometryTag && groups_[i].name == anchor_name) return groups_[i].geom_group; if (anchor_type == PartitionTag && groups_[i].name == anchor_name) return groups_[i].part_group; } return DEFAULT_INVALID_TAG; } const gctl::array &gctl::mesh_io::get_nodes() { return nodes_; } const gctl::array &gctl::mesh_io::get_elems() { return elems_; } size_t gctl::mesh_io::element_size(element_type_enum e_type) { size_t s = 0; for (size_t i = 0; i < groups_.size(); i++) { if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled) { s += groups_[i].elem_ptrs.size(); } } return s; } size_t gctl::mesh_io::element_size(element_tag_enum tag_type, int tag) { size_t s = 0; 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) { s += groups_[i].elem_ptrs.size(); } } return s; } size_t gctl::mesh_io::element_size(std::string phys_name) { size_t s = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].name == phys_name) { s += groups_[i].elem_ptrs.size(); } } return s; } 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; } void gctl::mesh_io::export_elements_to(array &tris, std::string phys_name) { size_t s = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _3NodeTriangle && (groups_[i].name == phys_name || phys_name == "All")) { s += groups_[i].elem_ptrs.size(); } } tris.resize(s); s = 0; for (size_t i = 0; i < elems_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _3NodeTriangle && (groups_[i].name == phys_name || phys_name == "All")) { for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++) { tris[s].id = groups_[i].elem_ptrs[e]->id; tris[s].vert[0] = groups_[i].elem_ptrs[e]->vert_ptrs[0]; tris[s].vert[1] = groups_[i].elem_ptrs[e]->vert_ptrs[1]; tris[s].vert[2] = groups_[i].elem_ptrs[e]->vert_ptrs[2]; s++; } } } return; } void gctl::mesh_io::export_elements_to(array &tris, element_tag_enum tag_type, int tag) { size_t s = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _3NodeTriangle && ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || (tag_type == GeometryTag && groups_[i].geom_group == tag) || (tag_type == PartitionTag && groups_[i].part_group == tag))) { s += groups_[i].elem_ptrs.size(); } } tris.resize(s); s = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _3NodeTriangle && ((tag_type == PhysicalTag && groups_[i].phys_group == tag) || (tag_type == GeometryTag && groups_[i].geom_group == tag) || (tag_type == PartitionTag && groups_[i].part_group == tag))) { for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++) { tris[s].id = groups_[i].elem_ptrs[e]->id; tris[s].vert[0] = groups_[i].elem_ptrs[e]->vert_ptrs[0]; tris[s].vert[1] = groups_[i].elem_ptrs[e]->vert_ptrs[1]; tris[s].vert[2] = groups_[i].elem_ptrs[e]->vert_ptrs[2]; s++; } } } return; } void gctl::mesh_io::export_elements_to(array &tets, std::string phys_name) { size_t s = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && (groups_[i].name == phys_name || phys_name == "All")) { s += groups_[i].elem_ptrs.size(); } } tets.resize(s); s = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && (groups_[i].name == phys_name || phys_name == "All")) { for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++) { tets[s].set(groups_[i].elem_ptrs[e]->vert_ptrs[0], groups_[i].elem_ptrs[e]->vert_ptrs[1], groups_[i].elem_ptrs[e]->vert_ptrs[2], groups_[i].elem_ptrs[e]->vert_ptrs[3], groups_[i].elem_ptrs[e]->id); s++; } } } return; } void gctl::mesh_io::export_elements_to(array &tets, element_tag_enum tag_type, int tag) { size_t s = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && (tag_type == PhysicalTag && groups_[i].phys_group == tag) || (tag_type == GeometryTag && groups_[i].geom_group == tag) || (tag_type == PartitionTag && groups_[i].part_group == tag)) { s += groups_[i].elem_ptrs.size(); } } tets.resize(s); s = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && (tag_type == PhysicalTag && groups_[i].phys_group == tag) || (tag_type == GeometryTag && groups_[i].geom_group == tag) || (tag_type == PartitionTag && groups_[i].part_group == tag)) { for (size_t e = 0; e < groups_[i].elem_ptrs.size(); e++) { tets[s].set(groups_[i].elem_ptrs[e]->vert_ptrs[0], groups_[i].elem_ptrs[e]->vert_ptrs[1], groups_[i].elem_ptrs[e]->vert_ptrs[2], groups_[i].elem_ptrs[e]->vert_ptrs[3], groups_[i].elem_ptrs[e]->id); s++; } } } return; } 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; } 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]; } void gctl::mesh_io::add_node_data(std::string name, const array &data) { size_t s = nodes_.size(); if (data.size()!= s) throw std::runtime_error("[gctl::mesh_io::create_node_data] Incompatible data size."); int d_id = if_saved_data(name, NodeData); if (d_id != -1) { for (size_t i = 0; i < s; i++) { datas_[d_id].val[i] = data[i]; } return; } meshio_data new_data; new_data.enabled = true; new_data.d_type = NodeData; 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); for (size_t i = 0; i < s; i++) { new_data.tar_ptrs[i] = nodes_.get(i); new_data.val[i] = data[i]; } datas_.push_back(new_data); return; } void gctl::mesh_io::add_node_data(std::string name, const array &data, const array &boolen) { size_t s = nodes_.size(); if (data.size()!= s || boolen.size() != s) throw std::runtime_error("[gctl::mesh_io::create_node_data] Incompatible data size."); s = 0; for (size_t i = 0; i < nodes_.size(); i++) { if (boolen[i]) s++; } int d_id = if_saved_data(name, NodeData); if (d_id != -1) { datas_[d_id].val.resize(s); datas_[d_id].tar_ptrs.resize(s, nullptr); datas_[d_id].int_tag[2] = s; s = 0; for (size_t i = 0; i < nodes_.size(); i++) { if (boolen[i]) { datas_[d_id].tar_ptrs[s] = nodes_.get(i); datas_[d_id].val[s] = data[i]; s++; } } return; } meshio_data new_data; new_data.enabled = true; new_data.d_type = NodeData; 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); s = 0; for (size_t i = 0; i < nodes_.size(); i++) { if (boolen[i]) { new_data.tar_ptrs[s] = nodes_.get(i); new_data.val[s] = data[i]; s++; } } datas_.push_back(new_data); return; } void gctl::mesh_io::add_element_data(std::string name, const array &data, element_type_enum e_type) { size_t e = 0; for (size_t i = 0; i < groups_.size(); i++) { if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled) e += groups_[i].elem_ptrs.size(); } if (data.size() != e) throw std::runtime_error("[gctl::mesh_io::create_element_data] Incompatible data size."); int d_id = if_saved_data(name, ElemData); if (d_id != -1) { datas_[d_id].int_tag[2] = e; datas_[d_id].tar_ptrs.resize(e, nullptr); datas_[d_id].val.resize(e, 0.0); e = 0; for (size_t i = 0; i < groups_.size(); i++) { if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled) { for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) { datas_[d_id].tar_ptrs[e] = groups_[i].elem_ptrs[j]; datas_[d_id].val[e] = data[e]; e++; } } } return; } meshio_data new_data; new_data.enabled = true; new_data.d_type = ElemData; 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] = e; new_data.tar_ptrs.resize(e, nullptr); new_data.val.resize(e, 0.0); e = 0; for (size_t i = 0; i < groups_.size(); i++) { if ((e_type == NotSet || groups_[i].type == e_type) && groups_[i].enabled) { for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) { new_data.tar_ptrs[e] = groups_[i].elem_ptrs[j]; new_data.val[e] = data[e]; e++; } } } datas_.push_back(new_data); return; } void gctl::mesh_io::add_element_data(std::string name, const array &data, element_tag_enum tag_type, int tag) { size_t e = 0; 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) { e += groups_[i].elem_ptrs.size(); } } if (data.size() != e) throw std::runtime_error("[gctl::mesh_io::create_element_data] Incompatible data size."); int d_id = if_saved_data(name, ElemData); if (d_id != -1) { array more_elem_ptrs(e, nullptr); array more_val(e, 0.0); e = 0; 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) { for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) { more_elem_ptrs[e] = groups_[i].elem_ptrs[j]; more_val[e] = data[e]; e++; } } } datas_[d_id].int_tag[2] += e; datas_[d_id].tar_ptrs.concat(more_elem_ptrs); datas_[d_id].val.concat(more_val); return; } meshio_data new_data; new_data.enabled = true; new_data.d_type = ElemData; 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] = e; new_data.tar_ptrs.resize(e, nullptr); new_data.val.resize(e, 0.0); e = 0; 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) { for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) { new_data.tar_ptrs[e] = groups_[i].elem_ptrs[j]; new_data.val[e] = data[e]; e++; } } } datas_.push_back(new_data); return; } void gctl::mesh_io::add_element_data(std::string name, std::string phys_name, const array &data) { size_t e = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].name == phys_name) { e += groups_[i].elem_ptrs.size(); } } if (data.size() != e) throw std::runtime_error("[gctl::mesh_io::create_element_data] Incompatible data size."); int d_id = if_saved_data(name, ElemData); if (d_id != -1) { array more_elem_ptrs(e, nullptr); array more_val(e, 0.0); e = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].name == phys_name) { for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) { more_elem_ptrs[e] = groups_[i].elem_ptrs[j]; more_val[e] = data[e]; e++; } } } datas_[d_id].int_tag[2] += e; datas_[d_id].tar_ptrs.concat(more_elem_ptrs); datas_[d_id].val.concat(more_val); return; } meshio_data new_data; new_data.enabled = true; new_data.d_type = ElemData; 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] = e; new_data.tar_ptrs.resize(e, nullptr); new_data.val.resize(e, 0.0); e = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].name == phys_name) { for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) { new_data.tar_ptrs[e] = groups_[i].elem_ptrs[j]; new_data.val[e] = data[e]; e++; } } } datas_.push_back(new_data); return; } void gctl::mesh_io::add_element_data(std::string name, std::string phys_name, double phys_val) { size_t e = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].name == phys_name) { e += groups_[i].elem_ptrs.size(); } } int d_id = if_saved_data(name, ElemData); if (d_id != -1) { array more_elem_ptrs(e, nullptr); array more_val(e, 0.0); e = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].name == phys_name) { for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) { more_elem_ptrs[e] = groups_[i].elem_ptrs[j]; more_val[e] = phys_val; e++; } } } datas_[d_id].int_tag[2] += e; datas_[d_id].tar_ptrs.concat(more_elem_ptrs); datas_[d_id].val.concat(more_val); return; } meshio_data new_data; new_data.enabled = true; new_data.d_type = ElemData; 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] = e; new_data.tar_ptrs.resize(e, nullptr); new_data.val.resize(e, 0.0); e = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].name == phys_name) { for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) { new_data.tar_ptrs[e] = groups_[i].elem_ptrs[j]; new_data.val[e] = phys_val; e++; } } } datas_.push_back(new_data); 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_); 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; } 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); } 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(); for (size_t g = 0; g < groups_.size(); g++) { groups_[g].enabled = true; groups_[g].enable_elements(); } 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_); 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; } 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); } 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); } 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(); for (size_t g = 0; g < groups_.size(); g++) { groups_[g].enabled = true; groups_[g].enable_elements(); } 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); 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; } 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); 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); } } 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(); for (size_t g = 0; g < groups_.size(); g++) { groups_[g].enabled = true; groups_[g].enable_elements(); } 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 = vptr->c2s(); ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*M_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 = pc.c2s(); ps.rad -= ellipse_radius_2d(refR, refr, ps.lat*M_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; }