From cf8b98af4cfdb4e24fc1ae1d48221a880abfe163 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Wed, 23 Apr 2025 21:55:30 +0800 Subject: [PATCH] tmp --- example/CMakeLists.txt | 21 +- example/meshio_ex.cpp | 109 ++ example/tri2d_meshio_ex.cpp | 55 + lib/mesh.h | 3 + lib/mesh/gctl_mesh_config.h | 2 +- lib/mesh/mesh.cpp | 1 + lib/mesh/mesh.h | 3 +- lib/mesh/meshdata.cpp | 2 + lib/mesh/meshdata.h | 2 +- lib/mesh/meshio.cpp | 2024 ++++++++++++++++++++++++ lib/mesh/meshio.h | 734 +++++++++ lib/mesh/regular_grid.h | 5 +- lib/mesh/tri2d_mesh.h | 10 +- lib/mesh/tri_mesh.cpp | 254 +++ lib/mesh/tri_mesh.h | 85 + tool/gridmanager/gridmanager.cpp | 62 +- tool/gridmanager/gridmanager.h | 2 + tool/gridmanager/readme_gridmanager.md | 4 +- 18 files changed, 3307 insertions(+), 71 deletions(-) create mode 100644 example/meshio_ex.cpp create mode 100644 example/tri2d_meshio_ex.cpp create mode 100644 lib/mesh/meshio.cpp create mode 100644 lib/mesh/meshio.h create mode 100644 lib/mesh/tri_mesh.cpp create mode 100644 lib/mesh/tri_mesh.h diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 63c8227..60b624d 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -12,18 +12,19 @@ macro(add_example name switch) # 链接动态库 target_link_libraries(${name} PRIVATE ${GCTL_LIB}) - target_link_libraries(${name} PRIVATE ${GCTL_GRAPHIC_LIB}) target_link_libraries(${name} PRIVATE gctl_mesh) endif() endmacro() -add_example(mesh_ex1 ON) -add_example(mesh_ex2 ON) -add_example(mesh_ex3 ON) -add_example(mesh_ex4 ON) +add_example(mesh_ex1 OFF) +add_example(mesh_ex2 OFF) +add_example(mesh_ex3 OFF) +add_example(mesh_ex4 OFF) add_example(mesh_ex5 OFF) -add_example(mesh_ex6 ON) -add_example(mesh_ex7 ON) -add_example(mesh_ex8 ON) -add_example(mesh_ex9 ON) -add_example(mesh_ex10 ON) \ No newline at end of file +add_example(mesh_ex6 OFF) +add_example(mesh_ex7 OFF) +add_example(mesh_ex8 OFF) +add_example(mesh_ex9 OFF) +add_example(mesh_ex10 OFF) +add_example(meshio_ex OFF) +add_example(tri2d_meshio_ex ON) \ No newline at end of file diff --git a/example/meshio_ex.cpp b/example/meshio_ex.cpp new file mode 100644 index 0000000..39910c1 --- /dev/null +++ b/example/meshio_ex.cpp @@ -0,0 +1,109 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 "../lib/mesh/mesh_io.h" + +using namespace gctl; + +int main(int argc, char const *argv[]) try +{ + mesh_io mshio; + mshio.read_gmsh_v2_ascii("tmp/test"); + mshio.info(); + + array nodes; + array tris; + + mshio.select_elements("Bz (pT)", NodeData); + mshio.export_selected_to(tris, nodes); + + //meshio_data &d = mshio.get_data("Bz (pT)", NodeData); + //array dat = d.val; + array dat1 = mshio.get_selected_data("Bx (pT)", NodeData); + array dat2 = mshio.get_selected_data("By (pT)", NodeData); + array dat3 = mshio.get_selected_data("Bz (pT)", NodeData); + + mesh_io mshio2; + mshio2.import_from(tris, nodes); + mshio2.add_data("Bx (pT)", dat1, NodeData, OverWrite); + mshio2.add_data("By (pT)", dat2, NodeData, OverWrite); + mshio2.add_data("Bz (pT)", dat3, NodeData, OverWrite); + mshio2.save_gmsh_v2_ascii("tmp"); + mshio2.info(); + +/* + mshio.read_tetgen_ascii("tmp/ex1.1"); + mshio.edit_group(Disable, GeometryTag, 5); + mshio.edit_group(GeometryTag, 1, PhysicalTag, 1); + mshio.edit_group(GeometryTag, 2, PhysicalTag, 2); + mshio.edit_group(GeometryTag, 3, PhysicalTag, 3); + mshio.edit_group(GeometryTag, 4, PhysicalTag, 4); + mshio.edit_group(GeometryTag, 1, "Boundary"); + mshio.edit_group(GeometryTag, 2, "Body1"); + mshio.edit_group(GeometryTag, 3, "Body2"); + mshio.edit_group(GeometryTag, 4, "Body3"); + mshio.save_gmsh_v2_ascii("tmp/ex1.1"); +*/ +/* + mshio.read_gmsh_v2_ascii("tmp/ex1.1"); + mshio.convert_tags_to_data(GeometryTag); + + array body_val(mshio.element_size("Body2"), 2.0); + mshio.add_element_data("BodyValue", "Body2", body_val); + + array body_val2(mshio.element_size("Body3"), 1.0); + mshio.add_element_data("BodyValue", "Body3", body_val2); + + mshio.save_gmsh_v2_ascii("tmp/ex1.2"); + //mshio.save_vtk_legacy_ascii("tmp/ex1.1"); + mshio.info(); + + const array &nodes = mshio.get_nodes(); + + array body2_tets; + mshio.export_elements_to(body2_tets, "All"); + + gmshio gio; + gio.init_file("tmp.msh", Output); + gio.set_packed(NotPacked, Output); + gio.save_mesh(body2_tets, nodes); +*/ +/* + mshio.read_gmsh_v2_ascii("tmp/wjb.1"); + mshio.edit_group(Disable); + mshio.edit_group(Enable, GeometryTag, 3); + mshio.edit_group(Enable, GeometryTag, 8); + mshio.edit_group(Enable, GeometryTag, 9); + + mshio.save_gmsh_v2_ascii("tmp/wjb.2"); +*/ + return 0; +} +catch(std::exception &e) +{ + GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0); +} \ No newline at end of file diff --git a/example/tri2d_meshio_ex.cpp b/example/tri2d_meshio_ex.cpp new file mode 100644 index 0000000..db0ad93 --- /dev/null +++ b/example/tri2d_meshio_ex.cpp @@ -0,0 +1,55 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 "../lib/mesh/tri2d_mesh.h" + + using namespace gctl; + + int main(int argc, char const *argv[]) try + { + mesh_io mshio; + mshio.read_gmsh_v2_ascii("tmp/tri2d"); + mshio.info(); + + array nodes; + array tris; + + mshio.select_elements(_3NodeTriangle); + mshio.export_selected_to(tris, nodes); + _1d_array data = mshio.get_selected_data("Model slowness (s/km)", ElemData); + + triangle2d_mesh tri2dmesh; + tri2dmesh.init("Test", "This is a test.", nodes, tris); + tri2dmesh.add_data(ElemData, "Test Data", data); + tri2dmesh.show_info(); + tri2dmesh.save_gmsh_withdata("tri2d_out", OverWrite); + return 0; + } + catch(std::exception &e) + { + GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0); + } \ No newline at end of file diff --git a/lib/mesh.h b/lib/mesh.h index 69e7545..affae57 100644 --- a/lib/mesh.h +++ b/lib/mesh.h @@ -29,6 +29,8 @@ #define _GCTL_MESH_H #include "mesh/gctl_mesh_config.h" +#include "mesh/meshdata.h" +#include "mesh/meshio.h" #include "mesh/mesh.h" #include "mesh/regular_grid.h" #include "mesh/regular_mesh_2d.h" @@ -36,6 +38,7 @@ #include "mesh/linear_mesh_2d.h" #include "mesh/linear_mesh_3d.h" #include "mesh/tri2d_mesh.h" +#include "mesh/tri_mesh.h" #include "mesh/tet_mesh.h" #include "mesh/regular_mesh_sph_3d.h" diff --git a/lib/mesh/gctl_mesh_config.h b/lib/mesh/gctl_mesh_config.h index c2b6d06..8b15abd 100644 --- a/lib/mesh/gctl_mesh_config.h +++ b/lib/mesh/gctl_mesh_config.h @@ -1,3 +1,3 @@ -#define GCTL_MESH_INSTALL_PREFIX "/usr/local" +#define GCTL_MESH_INSTALL_PREFIX "/opt/stow/gctl_mesh" #define GCTL_MESH_EXPRTK #define GCTL_MESH_WAVELIB diff --git a/lib/mesh/mesh.cpp b/lib/mesh/mesh.cpp index 915a18e..afeb691 100644 --- a/lib/mesh/mesh.cpp +++ b/lib/mesh/mesh.cpp @@ -60,6 +60,7 @@ void gctl::base_mesh::clear() node_num_ = ele_num_ = 0; initialized_ = false; destroy_vector(datalist_); + datalist_.reserve(100); return; } diff --git a/lib/mesh/mesh.h b/lib/mesh/mesh.h index e6ff3b6..543e97b 100644 --- a/lib/mesh/mesh.h +++ b/lib/mesh/mesh.h @@ -29,8 +29,7 @@ #define _GCTL_BASE_MESH_H #include "meshdata.h" -#include "gctl/io.h" -#include "gctl/math.h" +#include "meshio.h" namespace gctl { diff --git a/lib/mesh/meshdata.cpp b/lib/mesh/meshdata.cpp index a0681a9..43afdd9 100644 --- a/lib/mesh/meshdata.cpp +++ b/lib/mesh/meshdata.cpp @@ -63,6 +63,8 @@ void gctl::meshdata::create(mesh_data_type_e in_loctype, mesh_data_value_e in_valtype, size_t size, std::string name, bool if_output, double nan_val) { + if (name == "") throw std::runtime_error("[gctl::meshdata::create] Invalid data name."); + name_ = name; loctype_ = in_loctype; valtype_ = in_valtype; diff --git a/lib/mesh/meshdata.h b/lib/mesh/meshdata.h index 78d6490..ba43288 100644 --- a/lib/mesh/meshdata.h +++ b/lib/mesh/meshdata.h @@ -101,7 +101,7 @@ namespace gctl bool if_output = true, double nan_val = GCTL_BDL_MAX); /** - * @brief 返回适量数据。 + * @brief 返回矢量数据。 */ array export_vector() const; diff --git a/lib/mesh/meshio.cpp b/lib/mesh/meshio.cpp new file mode 100644 index 0000000..3d9941d --- /dev/null +++ b/lib/mesh/meshio.cpp @@ -0,0 +1,2024 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 "meshio.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/mesh/meshio.h b/lib/mesh/meshio.h new file mode 100644 index 0000000..7477b26 --- /dev/null +++ b/lib/mesh/meshio.h @@ -0,0 +1,734 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 "gctl/math/refellipsoid.h" +#include "gctl/io/triangle_io.h" +#include "gctl/io/tetgen_io.h" +#include "gctl/io/gmsh_io.h" + +namespace gctl +{ + /** + * @brief 无效的索引缺省值。 + */ + #define DEFAULT_INVALID_TAG -9999 + + /** + * @brief 网格单元体名称枚举类型。 + */ + enum element_type_enum + { + NotSet, + _2NodeLine, + _3NodeTriangle, + _4NodeQuadrangle, + _4NodeTetrahedron, + _8NodeHexahedron, + _6NodePrism, + _5NodePyramid, + _3NodeSecondOrderLine, + _6NdoeSecondOrderLine, + _9NodeSecondOrderQuadrangle, + _10NodeSecondOrderTetrahedron, + _27NodeSecondOrderHexahedron, + _18NodeSecondOrderPrism, + _14NodeSecondOrderPyramid, + _1NodePoint, + _8NodeSecondOrderQuadrangle, + _20NdoeSecondOrderHexahedron, + _15NodeSecondOrderPrism, + _13NodeSecondOrderPyramid, + _9NodeThirdOrderIncompleteTriangle, + _10NdoeThirdOrderTriangle, + _12NodeFourthOrderIncompleteTriangle, + _15NodeFourthOrderTriangle, + _15NodeFifthOrderCompleteTriangle, + _21NodeFifthOrderCompleteTriangle, + _4NodeThirdOrderEdge, + _5NodeFourthOrderEdge, + _6NodeFifthOrderEdge, + _20NodeThirdOrderTetrahedron, + _35NodeFourthOrderTetrahedron, + _56NodeFifithOrderTetrahedron, + _64NodeThirdOrderHexahedron, + _125NodeFourthOrderHexahedron, + }; + + /** + * @brief 网格单元体标签类型枚举类型 + * + */ + enum element_tag_enum + { + PhysicalTag, // 元素的物理分组标签 + GeometryTag, // 元素的几何分组标签 + PartitionTag, // 元素的剖分分组标签 + NodeTag, // 顶点的标签（仅用于输出顶点标签数据） + }; + + /** + * @brief 网格单元体结构体 + * + */ + struct meshio_element + { + bool enabled; // 单元体是否有效 + int id; // 单元体编号 + element_type_enum type; // 单元体类型 + array 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 diff --git a/lib/mesh/regular_grid.h b/lib/mesh/regular_grid.h index 3875011..884b17f 100644 --- a/lib/mesh/regular_grid.h +++ b/lib/mesh/regular_grid.h @@ -28,7 +28,10 @@ #ifndef _GCTL_REGULAR_GRID_H #define _GCTL_REGULAR_GRID_H -#include "gctl/graphic.h" +#include "gctl/math/interpolate.h" +#include "gctl/graphic/gmt.h" +#include "gctl/io/surfer_io.h" +#include "gctl/io/netcdf_io.h" #include "mesh.h" #ifdef GCTL_MESH_EXPRTK diff --git a/lib/mesh/tri2d_mesh.h b/lib/mesh/tri2d_mesh.h index 551e7ff..2014bd6 100644 --- a/lib/mesh/tri2d_mesh.h +++ b/lib/mesh/tri2d_mesh.h @@ -25,8 +25,8 @@ * Also add information on how to contact you by electronic and paper mail. ******************************************************/ -#ifndef _GCTL_TRI_MESH_H -#define _GCTL_TRI_MESH_H +#ifndef _GCTL_TRI2D_MESH_H +#define _GCTL_TRI2D_MESH_H #include "mesh.h" @@ -69,7 +69,7 @@ namespace gctl void load_gmsh_groups(); /** - * @brief 将模型物理分组转换为模型数据组 + * @brief 将模型物理分组转换为模型数据组 * * @param datname 赋值给模型数据组的名称（如果数组不存在则新建） * @param phynames 物理组名称 @@ -82,7 +82,9 @@ namespace gctl array elems_; array groups_; _2i_vector elems_tag_; + + mesh_io fio_; }; }; -#endif //_GCTL_TRI_MESH_H \ No newline at end of file +#endif //_GCTL_TRI2D_MESH_H \ No newline at end of file diff --git a/lib/mesh/tri_mesh.cpp b/lib/mesh/tri_mesh.cpp new file mode 100644 index 0000000..29453af --- /dev/null +++ b/lib/mesh/tri_mesh.cpp @@ -0,0 +1,254 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 "tri_mesh.h" + +void gctl::triangle_mesh::init(std::string in_name, std::string in_info, const array &in_nodes_, + const array &in_triangles) +{ + check_initiated(true); // 检查是否已经初始化 + base_mesh::init(TRI_TET_MESH, MESH_3D, in_name, in_info); + + node_num_ = in_nodes_.size(); + ele_num_ = in_triangles.size(); + + nodes_.resize(node_num_); + for (int i = 0; i < node_num_; i++) + { + nodes_[i] = in_nodes_[i]; + } + + elems_.resize(ele_num_); + for (int i = 0; i < ele_num_; i++) + { + elems_[i].id = i; + for (int j = 0; j < 3; j++) + { + elems_[i].vert[j] = nodes_.get(in_triangles[i].vert[j]->id); + } + } + + initialized_ = true; + return; +} + +void gctl::triangle_mesh::show_mesh_dimension(std::ostream &os) const +{ + os << "Node num: " << node_num_ << std::endl; + os << "Elem num: " << ele_num_ << std::endl; + return; +} + +void gctl::triangle_mesh::load_binary(std::string filename) +{ + check_initiated(true); // 检查是否已经初始化 + + std::ifstream infile; + gctl::open_infile(infile, filename, ".2m", std::ios::in|std::ios::binary); + + // 读入网格头信息 + load_headinfo(infile, TRI_TET_MESH, MESH_3D); + + // 读入网格信息 + infile.read((char*)&node_num_, sizeof(int)); + infile.read((char*)&ele_num_, sizeof(int)); + + nodes_.resize(node_num_); + elems_.resize(ele_num_); + + for (int i = 0; i < node_num_; i++) + { + infile.read((char*)nodes_.get(i), sizeof(gctl::vertex3dc)); + } + + int in_index; + for (int i = 0; i < ele_num_; i++) + { + elems_[i].id = i; + for (int j = 0; j < 3; j++) + { + infile.read((char*)&in_index, sizeof(int)); + elems_[i].vert[j] = nodes_.get(in_index); + } + } + + // 读入模型数据单元 + load_datablock(infile); + + infile.close(); + initialized_ = true; + return; +} + +void gctl::triangle_mesh::save_binary(std::string filename) +{ + check_initiated(); // 检查是否已经初始化 + + std::ofstream outfile; + gctl::open_outfile(outfile, filename, ".2m", std::ios::out|std::ios::binary); + + // 首先输出网格的头信息 + save_headinfo(outfile); + + // 输出网格信息 + outfile.write((char*)&node_num_, sizeof(int)); + outfile.write((char*)&ele_num_, sizeof(int)); + + for (int i = 0; i < node_num_; i++) + { + outfile.write((char*)nodes_.get(i), sizeof(gctl::vertex3dc)); + } + + int in_index; + for (int i = 0; i < ele_num_; i++) + { + for (int j = 0; j < 3; j++) + { + in_index = elems_[i].vert[j]->id; + outfile.write((char*)&in_index, sizeof(int)); + } + } + + // 输出的模型数据单元 + save_datablock(outfile); + + outfile.close(); + return; +} + +gctl::triangle_mesh::triangle_mesh() : base_mesh::base_mesh(){} + +gctl::triangle_mesh::triangle_mesh(std::string in_name, std::string in_info, const array &in_nodes_, + const array &in_triangles) +{ + init(in_name, in_info, in_nodes_, in_triangles); +} + +gctl::triangle_mesh::~triangle_mesh(){} + +const gctl::array &gctl::triangle_mesh::get_nodes() const +{ + check_initiated(); // 检查是否已经初始化 + return nodes_; +} + +const gctl::array &gctl::triangle_mesh::get_elements() const +{ + check_initiated(); // 检查是否已经初始化 + return elems_; +} + +void gctl::triangle_mesh::load_gmsh(std::string filename, index_packed_e packed) +{ + meshio_.init_file(filename, Input); + meshio_.set_packed(packed, Input); + meshio_.read_mesh(elems_, nodes_, &elems_tag_); + + // 设置名称与信息等 + node_num_ = nodes_.size(); + ele_num_ = elems_.size(); + base_mesh::init(TRI_TET_MESH, MESH_3D, filename, "Imported from a .msh file."); + initialized_ = true; + return; +} + +void gctl::triangle_mesh::load_gmsh_groups() +{ + check_initiated(); // 检查是否已经初始化 + meshio_.read_physical_groups(groups_); + return; +} + +void gctl::triangle_mesh::save_gmsh(std::string filename, index_packed_e packed) +{ + meshio_.init_file(filename, Output); + meshio_.set_packed(packed, Output); + meshio_.save_mesh(elems_, nodes_); + return; +} + +void gctl::triangle_mesh::groups2data(std::string datname, _1s_vector phynames, _1d_vector phyvals) +{ + check_initiated(); // 检查是否已经初始化 + + if (phynames.size() != phyvals.size() || phynames.empty()) + { + throw std::invalid_argument("The size of phynames and phyvals must be the same and non-empty."); + } + + if (groups_.empty() || elems_tag_.empty()) + { + throw std::invalid_argument("No physical groups found."); + } + + if (elems_tag_.empty()) throw std::invalid_argument("No physical tags found."); + for (size_t i = 0; i < ele_num_; i++) + { + if (elems_tag_[i].empty()) throw std::invalid_argument("No physical tags found."); + } + + array tag_idx(phynames.size()); + for (size_t i = 0; i < phynames.size(); i++) + { + tag_idx[i] = meshio_.physical_name2tag(groups_, phynames[i]); + } + + if (saved(datname)) + { + meshdata &dat = get_data(datname); + if (dat.loctype_ != ElemData) throw std::invalid_argument("The data type must be ElemData."); + + for (size_t i = 0; i < ele_num_; i++) + { + for (size_t e = 0; e < phynames.size(); e++) + { + if (elems_tag_[i][0] == tag_idx[e]) + { + dat.datval_[i] = phyvals[e]; + break; + } + } + } + } + else + { + meshdata &newdat = add_data(ElemData, Scalar, datname, 0.0); + + for (size_t i = 0; i < ele_num_; i++) + { + for (size_t e = 0; e < phynames.size(); e++) + { + if (elems_tag_[i][0] == tag_idx[e]) + { + newdat.datval_[i] = phyvals[e]; + break; + } + } + } + } + return; +} \ No newline at end of file diff --git a/lib/mesh/tri_mesh.h b/lib/mesh/tri_mesh.h new file mode 100644 index 0000000..494ffd7 --- /dev/null +++ b/lib/mesh/tri_mesh.h @@ -0,0 +1,85 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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_TRI_MESH_H +#define _GCTL_TRI_MESH_H + +#include "mesh.h" + +namespace gctl +{ + class triangle_mesh : public base_mesh + { + public: + + /** + * 以下为mesh类型的虚函数实现 + */ + + void load_gmsh(std::string filename, index_packed_e packed = Packed); + void save_gmsh(std::string filename, index_packed_e packed = Packed); + + void load_binary(std::string filename); + void save_binary(std::string filename); + + void show_mesh_dimension(std::ostream &os) const; + + /** + * 以下为triangle_mesh的专有函数 + */ + + triangle_mesh(); + triangle_mesh(std::string in_name, std::string in_info, const array &in_nodes, + const array &in_triangles); + virtual ~triangle_mesh(); + + void init(std::string in_name, std::string in_info, const array &in_nodes, + const array &in_triangles); + + const array &get_nodes() const; + const array &get_elements() const; + + void load_gmsh_groups(); + + /** + * @brief 将模型物理分组转换为模型数据组 + * + * @param datname 赋值给模型数据组的名称（如果数组不存在则新建） + * @param phynames 物理组名称 + * @param phyvals 物理组值 + */ + void groups2data(std::string datname, _1s_vector phynames, _1d_vector phyvals); + + protected: + array nodes_; + array elems_; + array groups_; + _2i_vector elems_tag_; + }; +}; + +#endif //_GCTL_TRI2D_MESH_H \ No newline at end of file diff --git a/tool/gridmanager/gridmanager.cpp b/tool/gridmanager/gridmanager.cpp index 326c415..942d0fe 100644 --- a/tool/gridmanager/gridmanager.cpp +++ b/tool/gridmanager/gridmanager.cpp @@ -443,7 +443,7 @@ void save_text(const std::vector &cmd_units) void data_cloud(const std::vector &cmd_units) { - // data-cloud \ node|cell \ \ \ \ [\] + // data-cloud \ node|cell \ \ \ \ if (cmd_units.size() < 7) throw std::runtime_error("data-cloud: insufficient parameters."); gctl::array copy_str(7, "null"); @@ -456,34 +456,16 @@ void data_cloud(const std::vector &cmd_units) if (copy_str[1] == "node") d_type = NodeData; else if (copy_str[1] == "cell") d_type = ElemData; else throw std::runtime_error("open-surfer: invalid grid data type."); - - text_descriptor desc; - desc.file_name_ = copy_str[2]; - desc.file_ext_ = ".txt"; - desc.col_str_ = "0,1,2"; - desc.delimiter_ = ' '; - desc.att_sym_ = '#'; - desc.head_num_ = 0; - if (copy_str[6] != "null") - { - parse_string_to_value(copy_str[6], '/', true, desc.col_str_, - desc.delimiter_, desc.att_sym_, desc.head_num_); - } - std::vector posix_vec, posiy_vec, data_vec; - read_text2vectors(desc, posix_vec, posiy_vec, data_vec); + geodsv_io text_in; + text_in.load_text(copy_str[2]); - array posi_arr(posix_vec.size()); - array posi_val; + array posi_arr; + _1d_array data_vec; + text_in.get_column_point2dc(posi_arr, 0, 1); + text_in.get_column(data_vec, 2); - posi_val.input(data_vec); - for (size_t i = 0; i < posi_arr.size(); i++) - { - posi_arr[i].x = posix_vec[i]; - posi_arr[i].y = posiy_vec[i]; - } - - rg.load_data_cloud(posi_arr, posi_val, atof(copy_str[3].c_str()), atof(copy_str[4].c_str()), atof(copy_str[5].c_str()), copy_str[0], d_type); + rg.load_data_cloud(posi_arr, data_vec, atof(copy_str[3].c_str()), atof(copy_str[4].c_str()), atof(copy_str[5].c_str()), copy_str[0], d_type); return; } @@ -694,7 +676,7 @@ void get_stats(const std::vector &cmd_units) void get_profile(const std::vector &cmd_units) { - // profile \ \ {\,\ \,\ \}|{\ [\]} + // profile \ \ {\,\ \,\ \}|{\} if (cmd_units.size() < 4) throw std::runtime_error("profile: insufficient parameters."); gctl::array copy_str(5, "null"); @@ -714,29 +696,9 @@ void get_profile(const std::vector &cmd_units) } else // File exist { - text_descriptor desc; - desc.file_name_ = copy_str[2]; - desc.file_ext_ = ".txt"; - desc.att_sym_ = '#'; - desc.col_str_ = "0,1"; - desc.delimiter_ = ' '; - desc.head_num_ = 0; - - if (copy_str[3] != "null") - { - gctl::parse_string_to_value(copy_str[3], '/', true, desc.col_str_, desc.delimiter_, - desc.att_sym_, desc.head_num_); - } - - std::vector xs, ys; - gctl::read_text2vectors(desc, xs, ys); - - xys.resize(xs.size()); - for (size_t i = 0; i < xys.size(); i++) - { - xys[i].x = xs[i]; - xys[i].y = ys[i]; - } + geodsv_io fio; + fio.load_text(copy_str[2]); + fio.get_column_point2dc(xys, 0, 1); rg.extract_points(copy_str[0], xys, p_data); } diff --git a/tool/gridmanager/gridmanager.h b/tool/gridmanager/gridmanager.h index 186c592..87d4e3f 100644 --- a/tool/gridmanager/gridmanager.h +++ b/tool/gridmanager/gridmanager.h @@ -28,6 +28,8 @@ #ifndef GCTL_MESH_GRIDMANAGER_H #define GCTL_MESH_GRIDMANAGER_H +#include "gctl/io/term_io.h" +#include "gctl/io/dsv_io.h" #include "../../lib/mesh.h" using namespace gctl; diff --git a/tool/gridmanager/readme_gridmanager.md b/tool/gridmanager/readme_gridmanager.md index 6315131..a5c657e 100644 --- a/tool/gridmanager/readme_gridmanager.md +++ b/tool/gridmanager/readme_gridmanager.md @@ -22,8 +22,8 @@ Mathematic expresstions could be used to specify the grid ranges. 4. `save gmsh ` Save the grid using the Gmsh legacy format. This command will write all data that are allowed for outputing. 5. `save text ` Save the grid data to text file. -#### cloud \ node|cell \ \ \ \ [\] -Import a randomly distributed data could (namely a xyz data). Note that a grid must exist before the use this command. 'node' or 'cell' indicates whether the input could data should be node or cell registed. \ is the input data file. \,\ and \ defines an oval of which the command is used to interpolate the could data to grid points. Additional options for reading the \ could be given by the [\] argument which has a format of \/\/\/\. For example 0,1,2/,/#/0 means that the intput data could is stored int the first three columns separated by commas. The file has no head records and any line starts with '#' is an annotation. +#### cloud \ node|cell \ \ \ \ +Import a randomly distributed data could (namely a xyz data). Note that a grid must exist before the use this command. 'node' or 'cell' indicates whether the input could data should be node or cell registed. \ is the input data file. \,\ and \ defines an oval of which the command is used to interpolate the could data to grid points. #### gradient \ \ dx|dy \ Calculate directional gradient(s) of the selected data. The new gradient data will have the same