tmp

CMakeLists.txt

@@ -1,20 +1,18 @@
 cmake_minimum_required(VERSION 3.15.2)
 # 设置项目名称与语言
-project(GCTL_MESH VERSION 1.0)
+project(GCTL_MESH VERSION 2.0)
 # 添加配置配件编写的函数
 include(CMakePackageConfigHelpers)
 
+# ExprTK库在macOS 15.4中编译会有错误 添加以下命令对应检查项
+add_compile_options(-Wno-missing-template-arg-list-after-template-kw)
+
 message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME})
 message(STATUS "Install prefix: " ${CMAKE_INSTALL_PREFIX})
 message(STATUS "Processor: " ${CMAKE_HOST_SYSTEM_PROCESSOR})
 
 find_package(GCTL REQUIRED)
-find_package(GCTL_GRAPHIC REQUIRED)
 message(STATUS "GCTL Version: " ${GCTL_VERSION})
-message(STATUS "GCTL_GRAPHIC Version: " ${GCTL_GRAPHIC_VERSION})
-#if(${GCTL_VERSION} LESS 1.0)
-#    message(FATAL_ERROR "GCTL's version must be v1.0 or bigger.")
-#endif()
 
 option(GCTL_MESH_EXPRTK "Use the exprtk library." ON)
 option(GCTL_MESH_WAVELIB "Use the WaveLib library" ON)

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)
+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)

example/mesh_ex1.cpp

@@ -45,7 +45,10 @@ int main(int argc, char *argv[])
 	std::cout << "data name: " << data3.name_ << std::endl;
 	data3.datval_.show();
 
+	rgd.diff("data-4", "data-1", "data-3");
+
 	rgd.remove_data("data-1");
 	rgd.show_info();
+	rgd.save_gmsh_withdata("ex1_out", gctl::OverWrite, gctl::NotPacked);
 	return 0;
 }

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 <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "../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<vertex3dc> nodes;
+    array<triangle> tris;
+
+    mshio.select_elements("Bz (pT)", NodeData);
+    mshio.export_selected_to(tris, nodes);
+
+    //meshio_data &d = mshio.get_data("Bz (pT)", NodeData);
+    //array<double> dat = d.val;
+    array<double> dat1 = mshio.get_selected_data("Bx (pT)", NodeData);
+    array<double> dat2 = mshio.get_selected_data("By (pT)", NodeData);
+    array<double> 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<double> body_val(mshio.element_size("Body2"), 2.0);
+    mshio.add_element_data("BodyValue", "Body2", body_val);
+
+    array<double> 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<vertex3dc> &nodes = mshio.get_nodes();
+
+    array<tetrahedron> 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);
+}

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 <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+ #include "../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<vertex2dc> nodes;
+     array<triangle2d> 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);
+ }

lib/CMakeLists.txt

@@ -27,9 +27,7 @@ set_target_properties(gctl_mesh_static PROPERTIES CXX_STANDARD 17 CXX_STANDARD_R
 
 #连接动态库
 target_link_libraries(gctl_mesh PUBLIC ${GCTL_LIB})
-target_link_libraries(gctl_mesh PUBLIC ${GCTL_GRAPHIC_LIB})
 target_link_libraries(gctl_mesh_static ${GCTL_LIB})
-target_link_libraries(gctl_mesh_static ${GCTL_GRAPHIC_LIB})
 
 if(GCTL_MESH_WAVELIB)
 	target_link_libraries(gctl_mesh PUBLIC ${WaveLib_LIB})

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"
 

lib/mesh/linear_mesh_2d.cpp

@@ -276,9 +276,8 @@ const gctl::array<double> *gctl::linear_mesh_2d::get_ysizes() const
 
 void gctl::linear_mesh_2d::save_gmsh(std::string filename, index_packed_e packed)
 {
-	std::ofstream outfile;
-	gctl::open_outfile(outfile, filename, ".msh");
-	gctl::save2gmsh(outfile, elements, nodes, packed);
-	outfile.close();
+	meshio_.init_file(filename, Output);
+	meshio_.set_packed(packed, Output);
+	meshio_.save_mesh(elements, nodes);
 	return;
 }

lib/mesh/linear_mesh_2d.h

@@ -77,6 +77,6 @@ namespace gctl
 		array<vertex2dc> nodes;
 		array<rectangle2d> elements;
 	};
-}
+};
 
 #endif //_GCTL_LINEAR_MESH_2D_H

lib/mesh/linear_mesh_3d.cpp

@@ -366,9 +366,8 @@ const gctl::array<double> *gctl::linear_mesh_3d::get_zsizes() const
 
 void gctl::linear_mesh_3d::save_gmsh(std::string filename, index_packed_e packed)
 {
-	std::ofstream outfile;
-	gctl::open_outfile(outfile, filename, ".msh");
-	gctl::save2gmsh(outfile, elements, nodes, packed);
-	outfile.close();
+	meshio_.init_file(filename, Output);
+	meshio_.set_packed(packed, Output);
+	meshio_.save_mesh(elements, nodes);
 	return;
 }

lib/mesh/linear_mesh_3d.h

@@ -84,6 +84,6 @@ namespace gctl
 		array<vertex3dc> nodes;
 		array<block> elements;
 	};
-}
+};
 
 #endif //_GCTL_LINEAR_MESH_3D_H

lib/mesh/mesh.cpp

@@ -35,6 +35,15 @@ gctl::base_mesh::base_mesh()
     meshinfo_ = "Undefined";
     node_num_ = ele_num_ = 0;
 	initialized_ = false;
+	// 注意这里我们一定要预先为datalist_分配空间，因为
+	// datalist_为vector类型，所以元素增加时可能造成
+	// 原有的迭代器失效造成访问失败。一个简单的例子是做
+	// 网格加法时，我们通过get_data获取两个网格数据并将
+	// 计算结果通过add_data保存，这时候就可能出现vector
+	// 扩容造成get_data获取的数据无法访问，访问失败。
+	// 默认可以保存100个网格数据，应该是够用了。同时
+	// 在add_data中添加相应的错误提示。
+	datalist_.reserve(100);
 }
 
 gctl::base_mesh::~base_mesh()
@@ -51,6 +60,7 @@ void gctl::base_mesh::clear()
     node_num_ = ele_num_ = 0;
 	initialized_ = false;
     destroy_vector(datalist_);
+	datalist_.reserve(100);
 	return;
 }
 
@@ -182,6 +192,8 @@ gctl::meshdata &gctl::base_mesh::add_data(mesh_data_type_e in_loctype, mesh_data
     std::string name, double init_val, bool if_output, double nan_val)
 {
     check_initiated();
+	if (datalist_.size() == 100) throw std::runtime_error("[gctl::base_mesh] Maximal data number reached.");
+
     meshdata new_data(in_loctype, in_valtype, 0, name, if_output, nan_val);
 
     if (in_loctype == NodeData && in_valtype == Scalar) new_data.datval_.resize(node_num_, init_val);
@@ -200,6 +212,8 @@ gctl::meshdata &gctl::base_mesh::add_data(mesh_data_type_e in_loctype, std::stri
     const array<double> &init_arr, bool if_output, double nan_val)
 {
     check_initiated();
+	if (datalist_.size() == 100) throw std::runtime_error("[gctl::base_mesh] Maximal data number reached.");
+
     meshdata new_data(in_loctype, Scalar, 0, name, if_output, nan_val);
 
     if (in_loctype == NodeData)
@@ -222,6 +236,8 @@ gctl::meshdata &gctl::base_mesh::add_data(mesh_data_type_e in_loctype, std::stri
     const array<point3dc> &init_arr, bool if_output, double nan_val)
 {
     check_initiated();
+	if (datalist_.size() == 100) throw std::runtime_error("[gctl::base_mesh] Maximal data number reached.");
+
     meshdata new_data(in_loctype, Vector, 0, name, if_output, nan_val);
 
     if (in_loctype == NodeData)
@@ -258,6 +274,8 @@ gctl::meshdata &gctl::base_mesh::add_data(mesh_data_type_e in_loctype, std::stri
     const array<tensor> &init_arr, bool if_output, double nan_val)
 {
     check_initiated();
+	if (datalist_.size() == 100) throw std::runtime_error("[gctl::base_mesh] Maximal data number reached.");
+
     meshdata new_data(in_loctype, Tensor, 0, name, if_output, nan_val);
 
     if (in_loctype == NodeData)

lib/mesh/mesh.h

@@ -29,8 +29,7 @@
 #define _GCTL_BASE_MESH_H
 
 #include "meshdata.h"
-#include "gctl/io.h"
-#include "gctl/algorithm.h"
+#include "meshio.h"
 
 namespace gctl
 {
@@ -377,6 +376,6 @@ namespace gctl
 		 */
 		void save_datablock(std::ofstream &outfile);
 	};
-}
+};
 
 #endif //_GCTL_BASE_MESH_H

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;

lib/mesh/meshdata.h

@@ -29,9 +29,8 @@
 #define _GCTL_MESHDATA_H
 
 #include "gctl_mesh_config.h"
-#include "gctl/core.h"
-#include "gctl/geometry/point3c.h"
-#include "gctl/geometry/tensor.h"
+#include "gctl/poly/point3c.h"
+#include "gctl/poly/tensor.h"
 
 namespace gctl
 {
@@ -102,7 +101,7 @@ namespace gctl
 			bool if_output = true, double nan_val = GCTL_BDL_MAX);
 
 		/**
-		 * @brief 返回适量数据。
+		 * @brief 返回矢量数据。
 		 */
 		array<point3dc> export_vector() const;
 
@@ -139,6 +138,6 @@ namespace gctl
 		 */
 		void save_binary(std::ofstream &outfile);
     };
-}
+};
 
 #endif // _GCTL_MESHDATA_H

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

lib/mesh/regular_grid.cpp

@@ -652,10 +652,9 @@ void gctl::regular_grid::save_surfer_grid(std::string filename, std::string datn
 
 void gctl::regular_grid::save_gmsh(std::string filename, index_packed_e packed)
 {
-	std::ofstream outfile;
-	gctl::open_outfile(outfile, filename, ".msh");
-	gctl::save2gmsh(outfile, elements, nodes, packed);
-	outfile.close();
+	meshio_.init_file(filename, Output);
+	meshio_.set_packed(packed, Output);
+	meshio_.save_mesh(elements, nodes);
 	return;
 }
 
@@ -951,7 +950,6 @@ void gctl::regular_grid::diff(std::string newname, std::string datname, std::str
 	}
 
 	meshdata &new_data = add_data(data_type1, value_type1, newname, 0.0, true, GCTL_BDL_MAX);
-
 	for (size_t i = 0; i < new_data.datval_.size(); i++)
 	{
 		new_data.datval_[i] = data_1.datval_[i] - data_2.datval_[i];

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
@@ -328,6 +331,6 @@ namespace gctl
 		array<vertex2dc> nodes; ///< 规则网格的节点数组
 		array<rectangle2d> elements; ///< 规则网格的单元数组
 	};
-}
+};
 
 #endif //_GCTL_REGULAR_GRID_H

lib/mesh/regular_mesh_2d.cpp

@@ -213,9 +213,8 @@ double gctl::regular_mesh_2d::get_ysize() const
 
 void gctl::regular_mesh_2d::save_gmsh(std::string filename, index_packed_e packed)
 {
-	std::ofstream outfile;
-	gctl::open_outfile(outfile, filename, ".msh");
-	gctl::save2gmsh(outfile, elements, nodes, packed);
-	outfile.close();
+	meshio_.init_file(filename, Output);
+	meshio_.set_packed(packed, Output);
+	meshio_.save_mesh(elements, nodes);
 	return;
 }

lib/mesh/regular_mesh_2d.h

@@ -74,6 +74,6 @@ namespace gctl
 		array<vertex2dc> nodes;
 		array<rectangle2d> elements;
 	};
-}
+};
 
 #endif //_GCTL_REGULAR_MESH_2D_H

lib/mesh/regular_mesh_3d.cpp

@@ -260,9 +260,8 @@ double gctl::regular_mesh_3d::get_zsize() const
 
 void gctl::regular_mesh_3d::save_gmsh(std::string filename, index_packed_e packed)
 {
-	std::ofstream outfile;
-	gctl::open_outfile(outfile, filename, ".msh");
-	gctl::save2gmsh(outfile, elements, nodes, packed);
-	outfile.close();
+	meshio_.init_file(filename, Output);
+	meshio_.set_packed(packed, Output);
+	meshio_.save_mesh(elements, nodes);
 	return;
 }

lib/mesh/regular_mesh_3d.h

@@ -77,6 +77,6 @@ namespace gctl
 		array<vertex3dc> nodes;
 		array<block> elements;
 	};
-}
+};
 
 #endif //_GCTL_REGULAR_MESH_3D_H

lib/mesh/regular_mesh_sph_3d.cpp

@@ -271,16 +271,16 @@ void gctl::regular_mesh_sph_3d::save_gmsh(std::string filename, index_packed_e p
 	gctl::array<vertex3dc> tmp_nodes(nodes.size());
 	for (int i = 0; i < nodes.size(); i++)
 	{
-		tmp_c = nodes[i].s2c();
+		tmp_c = s2c(nodes[i]);
 		tmp_nodes[i].id = i;
 		tmp_nodes[i].x = tmp_c.x;
 		tmp_nodes[i].y = tmp_c.y;
 		tmp_nodes[i].z = tmp_c.z;
 	}
-	std::ofstream outfile;
-	gctl::open_outfile(outfile, filename, ".msh");
-	gctl::save2gmsh(outfile, elements, tmp_nodes, packed);
-	outfile.close();
+
+	meshio_.init_file(filename, Output);
+	meshio_.set_packed(packed, Output);
+	meshio_.save_mesh(elements, tmp_nodes);
 	tmp_nodes.clear();
 	return;
 }

lib/mesh/regular_mesh_sph_3d.h

@@ -76,6 +76,6 @@ namespace gctl
         array<vertex3ds> nodes;
         array<tesseroid> elements;
     };
-}
+};
 
 #endif // _GCTL_REGULAR_MESH_SPH_3D_H

lib/mesh/tet_mesh.cpp

@@ -197,9 +197,8 @@ void gctl::tetrahedron_mesh::load_gmsh(std::string filename, index_packed_e pack
 
 void gctl::tetrahedron_mesh::save_gmsh(std::string filename, index_packed_e packed)
 {
-	std::ofstream outfile;
-	gctl::open_outfile(outfile, filename, ".msh");
-	gctl::save2gmsh(outfile, elements, nodes, packed);
-	outfile.close();
+	meshio_.init_file(filename, Output);
+	meshio_.set_packed(packed, Output);
+	meshio_.save_mesh(elements, nodes);
 	return;
 }

lib/mesh/tet_mesh.h

@@ -69,6 +69,6 @@ namespace gctl
 		array<vertex3dc> nodes;
 		array<tetrahedron> elements;
 	};
-}
+};
 
 #endif //_GCTL_TET_MESH_H

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"
 
@@ -82,7 +82,9 @@ namespace gctl
 		array<triangle2d> elems_;
 		array<gmsh_physical_group> groups_;
 		_2i_vector elems_tag_;
-	};
-}
 
-#endif //_GCTL_TRI_MESH_H
+		mesh_io fio_;
+	};
+};
+
+#endif //_GCTL_TRI2D_MESH_H

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 <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "tri_mesh.h"
+
+void gctl::triangle_mesh::init(std::string in_name, std::string in_info, const array<vertex3dc> &in_nodes_, 
+    const array<triangle> &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<vertex3dc> &in_nodes_, 
+    const array<triangle> &in_triangles)
+{
+    init(in_name, in_info, in_nodes_, in_triangles);
+}
+
+gctl::triangle_mesh::~triangle_mesh(){}
+
+const gctl::array<gctl::vertex3dc> &gctl::triangle_mesh::get_nodes() const
+{
+    check_initiated(); // 检查是否已经初始化
+    return nodes_;
+}
+
+const gctl::array<gctl::triangle> &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<int> 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;
+}

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 <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_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<vertex3dc> &in_nodes, 
+            const array<triangle> &in_triangles);
+        virtual ~triangle_mesh();
+
+        void init(std::string in_name, std::string in_info, const array<vertex3dc> &in_nodes, 
+            const array<triangle> &in_triangles);
+
+        const array<vertex3dc> &get_nodes() const;
+        const array<triangle> &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<vertex3dc> nodes_;
+        array<triangle> elems_;
+        array<gmsh_physical_group> groups_;
+        _2i_vector elems_tag_;
+    };
+};
+
+#endif //_GCTL_TRI2D_MESH_H

tool/gridmanager/gridmanager.cpp

@@ -443,7 +443,7 @@ void save_text(const std::vector<std::string> &cmd_units)
 
 void data_cloud(const std::vector<std::string> &cmd_units)
 {
-    // data-cloud \<new-data-name\> node|cell \<data-could-file\> \<search-x\> \<search-y\> \<search-deg\> [\<text_descriptor\>]
+    // data-cloud \<new-data-name\> node|cell \<data-could-file\> \<search-x\> \<search-y\> \<search-deg\>
     if (cmd_units.size() < 7) throw std::runtime_error("data-cloud: insufficient parameters.");
 
     gctl::array<std::string> copy_str(7, "null");
@@ -457,33 +457,15 @@ void data_cloud(const std::vector<std::string> &cmd_units)
     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_);
-    }
+    geodsv_io text_in;
+    text_in.load_text(copy_str[2]);
 
-    std::vector<double> posix_vec, posiy_vec, data_vec;
-    read_text2vectors(desc, posix_vec, posiy_vec, data_vec);
+    array<point2dc> posi_arr;
+    _1d_array data_vec;
+    text_in.get_column_point2dc(posi_arr, 0, 1);
+    text_in.get_column(data_vec, 2);
 
-    array<point2dc> posi_arr(posix_vec.size());
-    array<double> posi_val;
-
-    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;
 }
 
@@ -530,31 +512,19 @@ void wavelet(const std::vector<std::string> &cmd_units)
 
 void sum_data(const std::vector<std::string> &cmd_units)
 {
-    // sum-data \<new-data-name\> \<data-name1\> \<data-name2\>
+    // sum-data \<new-data\> \<data1\> \<data2\>
     if (cmd_units.size() < 4) throw std::runtime_error("sum: insufficient parameters.");
 
-    gctl::array<std::string> copy_str(3, "null");
-    for (size_t i = 1; i <= GCTL_MIN(cmd_units.size() - 1, 3); i++)
-    {
-        copy_str[i - 1] = cmd_units[i];
-    }
-
-    rg.sum(copy_str[0], copy_str[1], copy_str[2]);
+    rg.sum(cmd_units[1], cmd_units[2], cmd_units[3]);
     return;
 }
 
 void diff_data(const std::vector<std::string> &cmd_units)
 {
-    // sum-data \<new-data-name\> \<data-name1\> \<data-name2\>
+    // diff \<new-data\> \<data1\> \<data2\>
     if (cmd_units.size() < 4) throw std::runtime_error("diff: insufficient parameters.");
 
-    gctl::array<std::string> copy_str(3, "null");
-    for (size_t i = 1; i <= GCTL_MIN(cmd_units.size() - 1, 3); i++)
-    {
-        copy_str[i - 1] = cmd_units[i];
-    }
-
-    rg.diff(copy_str[0], copy_str[1], copy_str[2]);
+    rg.diff(cmd_units[1], cmd_units[2], cmd_units[3]);
     return;
 }
 
@@ -706,7 +676,7 @@ void get_stats(const std::vector<std::string> &cmd_units)
 
 void get_profile(const std::vector<std::string> &cmd_units)
 {
-    // profile \<data-name\> \<out-file\> {\<x\>,\<y\> \<x\>,\<y\> \<num\>}|{\<xy-file\> [\<file-format\>]}
+    // profile \<data-name\> \<out-file\> {\<x\>,\<y\> \<x\>,\<y\> \<num\>}|{\<xy-file\>}
     if (cmd_units.size() < 4) throw std::runtime_error("profile: insufficient parameters.");
 
     gctl::array<std::string> copy_str(5, "null");
@@ -726,29 +696,9 @@ void get_profile(const std::vector<std::string> &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<double> 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);
     }

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;

tool/gridmanager/readme_gridmanager.md

@@ -22,8 +22,8 @@ Mathematic expresstions could be used to specify the grid ranges.
 4. `save gmsh <file>` Save the grid using the Gmsh legacy format. This command will write all data that are allowed for outputing.
 5. `save text <file>` Save the grid data to text file.
 
-#### cloud \<new-data-name\> node|cell \<data-could-file\> \<search-x\> \<search-y\> \<search-deg\>  [\<text_descriptor\>]
-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. \<data-could-file\> is the input data file. \<search-x\>,\<search-y\> and \<search-deg\> defines an oval of which the command is used to interpolate the could data to grid points. Additional options for reading the \<data-could-file\> could be given by the [\<text_descriptor\>] argument which has a format of \<order_str\>/\<delimiter\>/\<annotate\>/\<head_record\>. 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 \<new-data-name\> node|cell \<data-could-file\> \<search-x\> \<search-y\> \<search-deg\>
+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. \<data-could-file\> is the input data file. \<search-x\>,\<search-y\> and \<search-deg\> defines an oval of which the command is used to interpolate the could data to grid points.
 
 #### gradient \<data-name\> \<new-data-name\> dx|dy \<order\>
 Calculate directional gradient(s) of the selected data. The new gradient data will have the same
@@ -32,10 +32,10 @@ data type of the used data. 'dx' and 'dy' give the calculation directions. highe
 #### wavelet \<data-name\> \<wavelet-name\> \<order\> [\<show-summary\>]
 Using a wavelet \<wavelet-name\> to decompose the selecte data. Avaiable wavelets include 'bd2'... .The order of decomposition could be set using \<order\>. And set \<show-summary\> to yes to see the summary.
 
-#### sum \<new-data-name\> \<data-name1\> \<data-name2\>
+#### sum \<new-data\> \<data1\> \<data2\>
 Calculate the sum of the two data. The arguments' format should be clear enough.
 
-#### diff \<new-data-name\> \<data-name1\> \<data-name2\>
+#### diff \<new-data\> \<data1\> \<data2\>
 Calculate the difference of the two data. The arguments' format should be clear enough.
 
 #### bool \<new-data-name\> \<data-name1> \<bool-data\> [reverse]