Merge pull request 'dev_yi' (#1) from dev_yi into main

Reviewed-on: #1

README.md

@@ -1,2 +1,245 @@
-# gctl_mesh
+![logo](GCTL_logo.jpg)
+
+# GCTL Mesh Library
+
+GCTL Mesh是地球物理计算工具库(Geophysical Computational Tools & Library)的网格处理模块，提供了多种类型的网格实现和操作功能。
+
+## 功能特点
+
+- 支持多种网格类型的创建和操作
+- 提供2D和3D网格的实现
+- 包含规则网格和非规则网格的处理
+- 支持球面网格的处理
+- 支持网格数据的导入导出（包括Gmsh格式和二进制格式）
+- 支持标量、矢量和张量数据的处理
+
+## 代码库结构
+
+### 核心库 (lib/)
+```
+lib/
+├── mesh/               # 核心网格实现
+│   ├── mesh.h         # 网格基类
+│   ├── meshdata.h     # 网格数据结构
+│   ├── regular_grid.h # 规则网格实现
+│   ├── linear_mesh_2d.h  # 2D线性网格
+│   ├── linear_mesh_3d.h  # 3D线性网格
+│   ├── tri_mesh.h     # 三角形网格
+│   ├── tet_mesh.h     # 四面体网格
+│   └── regular_mesh_sph_3d.h  # 球面网格
+├── io/                 # 输入输出处理
+│   ├── gmsh_io.h      # Gmsh格式支持
+│   ├── netcdf_io.h    # NetCDF格式支持
+│   └── binary_io.h    # 二进制格式支持
+└── utils/             # 工具函数
+    ├── math_utils.h   # 数学工具
+    └── geo_utils.h    # 几何工具
+```
+
+### 依赖关系
+- `mesh.h` → 基类，定义网格接口
+- `meshdata.h` → 依赖 mesh.h，实现数据管理
+- 具体网格类 → 继承自 mesh.h，实现特定网格功能
+
+## 实现细节
+
+### 数据结构
+
+#### 网格节点 (Node)
+```cpp
+struct Node {
+    int id;            // 节点ID
+    double x, y, z;    // 节点坐标
+    array<int> elems;  // 相邻单元列表
+};
+```
+
+#### 网格单元 (Element)
+```cpp
+struct Element {
+    int id;                 // 单元ID
+    array<Node*> vertices;  // 顶点列表
+    array<Element*> neighbors;  // 相邻单元
+};
+```
+
+#### 数据管理 (MeshData)
+```cpp
+class MeshData {
+    string name;            // 数据名称
+    DataType type;         // 数据类型
+    array<double> values;  // 数据值
+    bool output_enabled;   // 输出控制
+};
+```
+
+### 核心功能实现
+
+#### 网格生成
+- 规则网格：基于给定维度和间距生成
+- 非规则网格：支持从外部文件导入或程序化生成
+- 球面网格：基于经纬度和径向参数生成
+
+#### 数据操作
+- 插值：支持线性插值和高阶插值
+- 梯度计算：支持一阶和二阶导数
+- 数据过滤：支持平滑和降噪操作
+
+#### 文件格式支持
+- Gmsh格式：`.msh`文件，支持网格和数据的导入导出
+- NetCDF格式：`.nc`文件，支持结构化数据的存储
+- 二进制格式：`.bin`文件，用于高效的数据存储和读取
+
+### 性能优化
+
+#### 数据结构优化
+- 使用连续内存存储提高访问效率
+- 采用邻接表存储拓扑关系
+- 使用哈希表加速查找操作
+
+#### 算法优化
+- 网格生成采用增量构建策略
+- 查找操作使用空间分区技术
+- 并行计算支持（可选）
+
+## 文档
+
+### 基础类
+- [`mesh.h`](doc/mesh.md): 网格基类，定义了基本的网格接口和功能
+- [`meshdata.h`](doc/meshdata.md): 网格数据结构的定义和处理
+
+### 规则网格
+- [`regular_grid.h`](doc/regular_grid.md): 通用规则网格实现，支持2D和3D
+- [`regular_mesh_2d.h`](doc/regular_mesh_2d.md): 二维规则网格，专注于2D结构化网格
+- [`regular_mesh_3d.h`](doc/regular_mesh_3d.md): 三维规则网格，专注于3D结构化网格
+- [`regular_mesh_sph_3d.h`](doc/regular_mesh_sph_3d.md): 三维球面规则网格，适用于地球物理建模
+
+### 线性网格
+- [`linear_mesh_2d.h`](doc/linear_mesh_2d.md): 二维线性网格，用于表示2D空间中的线段
+- [`linear_mesh_3d.h`](doc/linear_mesh_3d.md): 三维线性网格，用于表示3D空间中的曲线
+
+### 非规则网格
+- [`tri_mesh.h`](doc/tri_mesh.md): 三角形网格，用于复杂2D区域的离散化
+- [`tet_mesh.h`](doc/tet_mesh.md): 四面体网格，用于复杂3D区域的离散化
+
+## 编译和安装
+
+### 系统要求
+- C++11或更高版本
+- CMake 3.10或更高版本
+- 可选依赖：
+  - NetCDF (用于NetCDF格式支持)
+  - MathGL (用于可视化支持)
+  - GMT (用于地理数据处理)
+
+### 编译步骤
+```bash
+mkdir build
+cd build
+cmake ..
+make
+make install
+```
+
+### 配置选项
+- `GCTL_ENABLE_NETCDF`: 启用NetCDF支持
+- `GCTL_ENABLE_MATHGL`: 启用MathGL可视化
+- `GCTL_ENABLE_GMT`: 启用GMT支持
+- `GCTL_BUILD_EXAMPLES`: 构建示例程序
+- `GCTL_BUILD_TESTS`: 构建测试程序
+
+## 应用场景
+
+- **规则网格**：适用于简单几何形状的离散化，如矩形/立方体区域
+- **线性网格**：
+  - 2D：用于表示边界、裂缝等线性特征
+  - 3D：用于表示输电线路、管道系统、裂缝网络等
+- **三角形网格**：适用于复杂2D区域的有限元分析、流体动力学等
+- **四面体网格**：适用于复杂3D区域的有限元分析、计算流体力学等
+- **球面网格**：特别适用于地球物理学、大气科学等全球尺度的模拟
+
+## 性能考虑
+
+### 内存管理
+- 大规模网格处理时注意内存使用
+- 支持数据流式处理
+- 提供内存使用估算功能
+
+### 计算效率
+- 网格操作采用高效算法
+- 支持并行计算（可选）
+- 提供性能分析工具
+
+## API示例
+
+详细的API使用示例请参考各个网格类型的文档。基本使用模式如下：
+
+```cpp
+// 创建网格对象
+gctl::regular_grid rgd;
+rgd.init("grid-1", "test grid", 4, 4, 0.0, 0.0, 1.0, 1.0);
+
+// 添加数据
+rgd.add_data(gctl::NodeData, gctl::Scalar, "temperature", 25.0);
+rgd.add_data(gctl::ElemData, gctl::Vector, "velocity", velocity_array);
+
+// 保存网格
+rgd.save_gmsh("output.msh", gctl::OverWrite);
+```
+
+## 许可证
+
+该项目采用双重许可证方案：
+1. GNU Lesser General Public License v2或更高版本
+2. 商业许可证（需单独申请）
+
+## 联系方式
+
+作者：Yi Zhang (yizhang-geo@zju.edu.cn)
+
+## 版权声明
+
+Copyright (c) 2023 Yi Zhang
+
+## 使用示例
+
+库提供了多个完整的示例程序，展示了不同功能的使用方法：
+
+### 基础网格操作
+```cpp
+// mesh_ex1.cpp - 基础网格创建和数据操作
+gctl::regular_grid rgd;
+rgd.init("grid-1", "null", 4, 4, 0.0, 0.0, 1.0, 1.0);  // 创建4x4网格
+rgd.add_data(gctl::NodeData, gctl::Scalar, "data-1", 2.5);  // 添加节点数据
+rgd.add_data(gctl::ElemData, gctl::Scalar, "data-2", 1.2);  // 添加单元数据
+```
+
+### 文件导入导出
+```cpp
+// mesh_ex5.cpp - NetCDF格式文件操作
+gctl::regular_grid rgd;
+rgd.load_netcdf_grid("data/input.nc", gctl::ElemData, "x", "y");  // 读取网格
+rgd.save_netcdf_grid("data/output.nc", "example");  // 保存网格
+```
+
+### 球面网格处理
+```cpp
+// mesh_ex10.cpp - 3D球面网格创建
+gctl::regular_mesh_sph_3d rm_3ds;
+rm_3ds.init("mesh-1", "null", 30.25, 30.25, 2005, 0.5, 0.5, 10, 40, 40, 50);
+rm_3ds.add_data(gctl::ElemData, gctl::Scalar, "data-1", 2.5);
+rm_3ds.save_gmsh("mesh_sample10", gctl::OverWrite, gctl::NotPacked);
+```
+
+更多示例可以在`example`目录下找到：
+- `mesh_ex1.cpp`: 基础网格创建和数据操作
+- `mesh_ex2.cpp`: 网格数据的导入导出
+- `mesh_ex3.cpp`: 网格变换和操作
+- `mesh_ex4.cpp`: 非结构化网格处理
+- `mesh_ex5.cpp`: NetCDF格式文件操作
+- `mesh_ex6.cpp`: 网格数据插值
+- `mesh_ex7.cpp`: 网格数据统计分析
+- `mesh_ex8.cpp`: 二维规则网格操作
+- `mesh_ex9.cpp`: 三维规则网格操作
+- `mesh_ex10.cpp`: 球面网格处理
 

backup/mesh.cpp

@@ -0,0 +1,864 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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 "mesh.h"
+
+gctl::base_mesh::base_mesh()
+{
+	initialized = false;
+	//std::clog << "A new mesh object is created." << std::endl;
+}
+
+gctl::base_mesh::~base_mesh()
+{
+	clear();
+	//std::clog << "A mesh object is destroyed." << std::endl;
+}
+
+void gctl::base_mesh::clear()
+{
+	meshdata *data_ptr;
+	if (!saved_data.empty())
+	{
+		for (iter = saved_data.begin(); iter != saved_data.end(); ++iter)
+		{
+			data_ptr = *iter;
+			meshdata::destroy(data_ptr);
+		}
+		saved_data.clear();
+	}
+
+	initialized = false;
+	return;
+}
+
+bool gctl::base_mesh::initiated() const
+{
+	return initialized;
+}
+
+bool gctl::base_mesh::saved(std::string datname) const
+{
+	if (saved_data.empty())
+	{
+		return false;
+	}
+	else
+	{
+		meshdata *data_ptr = nullptr;
+		// 这里我们需要使用常量迭代器
+		std::list<meshdata*>::const_iterator c_iter;
+		for (c_iter = saved_data.begin(); c_iter != saved_data.end(); ++c_iter)
+		{
+			data_ptr = *c_iter;
+			if (data_ptr->get_datname() == datname)
+			{
+				return true;
+			}
+		}
+		return false;
+	}
+}
+
+gctl::meshdata *gctl::base_mesh::get_data(std::string datname) const
+{
+	if (saved_data.empty())
+	{
+		throw std::runtime_error("[gctl::base_mesh] No data saved.");
+	}
+
+	meshdata *curr_data = nullptr;
+	std::list<meshdata*>::const_iterator c_iter;
+	for (c_iter = saved_data.begin(); c_iter != saved_data.end(); ++c_iter)
+	{
+		curr_data = *c_iter;
+		if (curr_data->get_datname() == datname)
+		{
+			return curr_data;
+		}
+	}
+
+	throw gctl::runtime_error("[gctl::base_mesh] No data found by the name: " + datname);
+}
+
+void gctl::base_mesh::get_all_data(array<meshdata*>& out_list) const
+{
+	if (saved_data.empty())
+	{
+		throw runtime_error("[gctl::base_mesh] No data saved.");
+	}
+
+	int c_count = 0;
+	out_list.resize(saved_data.size());
+
+	meshdata *curr_data = nullptr;
+	std::list<meshdata*>::const_iterator c_iter;
+	for (c_iter = saved_data.begin(); c_iter != saved_data.end(); ++c_iter)
+	{
+		curr_data = *c_iter;
+		out_list[c_count] = curr_data;
+		c_count++;
+	}
+	return;
+}
+
+void *gctl::base_mesh::get_datval(std::string datname) const
+{
+	meshdata *curr_data = get_data(datname);
+	return curr_data->get_datval_ptr();
+}
+
+void gctl::base_mesh::remove_data(std::string datname)
+{
+	if (saved_data.empty())
+	{
+		throw runtime_error("[gctl::base_mesh] No data saved.");
+	}
+
+	meshdata *curr_data;
+	for (iter = saved_data.begin(); iter != saved_data.end(); ++iter)
+	{
+		curr_data = *iter;
+		if (curr_data->get_datname() == datname)
+		{
+			meshdata::destroy(curr_data);
+			iter = saved_data.erase(iter);
+			//std::clog << "Meshdata: " << datname << " is destroyed." << std::endl;
+			break;
+		}
+	}
+	return;
+}
+
+void gctl::base_mesh::show_info(std::ostream &os) const
+{
+	if (meshtype == REGULAR_MESH) os << "Mesh: Regular | ";
+	if (meshtype == LINEAR_MESH)  os << "Mesh: Linear | ";
+	if (meshtype == TRI_TET_MESH) os << "Mesh: Unstructured | ";
+	if (meshtype == REGULAR_MESH_SPH) os << "Mesh: Regular (spherical) | ";
+	if (meshtype == LINEAR_MESH_SPH)  os << "Mesh: Linear (spherical) | ";
+	if (meshtype == TRI_TET_MESH_SPH) os << "Mesh: Unstructured (spherical) | ";
+	if (meshtype == REGULAR_GRID) os << "Mesh: Grid | ";
+
+	if (meshdim == MESH_2D) os << "2D" << std::endl;
+	else if (meshdim == MESH_3D) os << "3D" << std::endl;
+
+	os << "Name: " << meshname << std::endl;
+	os << "Info: " << meshinfo << std::endl;
+
+	show_mesh_dimension(os);
+
+	meshdata *curr_data;
+	std::list<meshdata*>::const_iterator c_iter;
+	for (c_iter = saved_data.begin(); c_iter != saved_data.end(); ++c_iter)
+	{
+		curr_data = *c_iter;
+		curr_data->show_info();
+	}
+	return;
+}
+
+void gctl::base_mesh::rename_data(std::string oldname, std::string newname)
+{
+	meshdata *curr_data = get_data(oldname);
+	curr_data->set_datname(newname);
+	return;
+}
+
+gctl::mesh_type_e gctl::base_mesh::get_meshtype() const
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	return meshtype;
+}
+
+gctl::mesh_dim_e gctl::base_mesh::get_meshdim() const
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	return meshdim;
+}
+
+int gctl::base_mesh::get_nodenum() const
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	return node_num;
+}
+
+int gctl::base_mesh::get_elenum() const
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	return ele_num;
+}
+
+int gctl::base_mesh::get_datanum() const
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	return saved_data.size();
+}
+
+std::string gctl::base_mesh::get_meshname() const
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	return meshname;
+}
+
+void gctl::base_mesh::set_meshname(std::string in_name)
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	if (in_name.empty())
+	{
+		throw std::runtime_error("[gctl::base_mesh] The input name is empty.");
+	}
+
+	meshname = in_name;
+	return;
+}
+
+std::string gctl::base_mesh::get_meshinfo() const
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	return meshinfo;
+}
+
+void gctl::base_mesh::set_meshinfo(std::string in_info)
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	if (in_info == "")
+	{
+		throw runtime_error("[gctl::base_mesh] The input info. is empty.");
+	}
+
+	meshinfo = in_info;
+	return;
+}
+
+gctl::meshdata *gctl::base_mesh::add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, double init_val)
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	meshdata *data_ptr;
+	if(saved(in_name))
+	{
+		data_ptr = get_data(in_name);
+		if (data_ptr->get_valtype() != Scalar || data_ptr->get_dattype() != in_type)
+		{
+			throw gctl::runtime_error("[gctl::base_mesh] A data object is already created with a different data or value type by the name:" + in_name);
+		}
+
+		// 存在一个同名 同数据类型 同赋值类型的数据 则将其数据设置为初始值
+		array<double>* val_ptr = (array<double>*) data_ptr->get_datval_ptr();
+		val_ptr->assign_all(init_val);
+		return data_ptr;
+	}
+
+	if (in_type == NodeData)
+	{
+		data_ptr = meshdata_scalar::create(in_name, in_type, node_num, if_output, init_val);
+	}
+	else if (in_type == ElemData || in_type == ElemData2D || in_type == ElemData3D)
+	{
+		data_ptr = meshdata_scalar::create(in_name, in_type, ele_num, if_output, init_val);
+	}
+
+	saved_data.push_back(data_ptr);
+	return data_ptr;
+}
+
+gctl::meshdata *gctl::base_mesh::add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, gctl::point3dc init_val)
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	meshdata *data_ptr;
+	if(saved(in_name))
+	{
+		data_ptr = get_data(in_name);
+		if (data_ptr->get_valtype() != Scalar || data_ptr->get_dattype() != in_type)
+		{
+			throw gctl::runtime_error("[gctl::base_mesh] A data object is already created with a different data or value type by the name:" + in_name);
+		}
+
+		// 存在一个同名 同数据类型 同赋值类型的数据 则将其数据设置为初始值
+		array<point3dc>* val_ptr = (array<point3dc>*) data_ptr->get_datval_ptr();
+		val_ptr->assign_all(init_val);
+		return data_ptr;
+	}
+
+	if (in_type == NodeData)
+	{
+		data_ptr = meshdata_vector::create(in_name, in_type, node_num, if_output, init_val);
+	}
+	else if (in_type == ElemData || in_type == ElemData2D || in_type == ElemData3D)
+	{
+		data_ptr = meshdata_vector::create(in_name, in_type, ele_num, if_output, init_val);
+	}
+
+	saved_data.push_back(data_ptr);
+	return data_ptr;
+}
+
+gctl::meshdata *gctl::base_mesh::add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, gctl::tensor init_val)
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	meshdata *data_ptr;
+	if(saved(in_name))
+	{
+		data_ptr = get_data(in_name);
+		if (data_ptr->get_valtype() != Scalar || data_ptr->get_dattype() != in_type)
+		{
+			throw gctl::runtime_error("[gctl::base_mesh] A data object is already created with a different data or value type by the name:" + in_name);
+		}
+
+		// 存在一个同名 同数据类型 同赋值类型的数据 则将其数据设置为初始值
+		array<tensor>* val_ptr = (array<tensor>*) data_ptr->get_datval_ptr();
+		val_ptr->assign_all(init_val);
+		return data_ptr;
+	}
+
+	if (in_type == NodeData)
+	{
+		data_ptr = meshdata_tensor::create(in_name, in_type, node_num, if_output, init_val);
+	}
+	else if (in_type == ElemData || in_type == ElemData2D || in_type == ElemData3D)
+	{
+		data_ptr = meshdata_tensor::create(in_name, in_type, ele_num, if_output, init_val);
+	}
+
+	saved_data.push_back(data_ptr);
+	return data_ptr;
+}
+
+gctl::meshdata *gctl::base_mesh::add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, mesh_data_value_e val_type)
+{
+	if (!initialized)
+	{
+		throw std::runtime_error("[gctl::base_mesh] Mesh not initialized.");
+	}
+
+	meshdata *data_ptr;
+	if(saved(in_name))
+	{
+		data_ptr = get_data(in_name);
+		if (data_ptr->get_valtype() != Scalar || data_ptr->get_dattype() != in_type)
+		{
+			throw gctl::runtime_error("[gctl::base_mesh] A data object is already created with a different data or value type by the name:" + in_name);
+		}
+
+		return data_ptr;
+	}
+
+	if (val_type == Scalar && in_type == NodeData)
+	{
+		data_ptr = meshdata_scalar::create(in_name, in_type, node_num, if_output, 0.0);
+	}
+	else if (val_type == Scalar && in_type == ElemData)
+	{
+		data_ptr = meshdata_scalar::create(in_name, in_type, ele_num, if_output, 0.0);
+	}
+	else if (val_type == Vector && in_type == NodeData)
+	{
+		point3dc init_val = {0.0, 0.0, 0.0};
+		data_ptr = meshdata_vector::create(in_name, in_type, node_num, if_output, init_val);
+	}
+	else if (val_type == Vector && in_type == ElemData)
+	{
+		point3dc init_val = {0.0, 0.0, 0.0};
+		data_ptr = meshdata_vector::create(in_name, in_type, ele_num, if_output, init_val);
+	}
+	else if (val_type == Tensor && in_type == NodeData)
+	{
+		tensor init_val = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+		data_ptr = meshdata_tensor::create(in_name, in_type, node_num, if_output, init_val);
+	}
+	else if (val_type == Tensor && in_type == ElemData)
+	{
+		tensor init_val = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+		data_ptr = meshdata_tensor::create(in_name, in_type, ele_num, if_output, init_val);
+	}
+
+
+	saved_data.push_back(data_ptr);
+	return data_ptr;
+}
+
+void gctl::base_mesh::init(std::string in_name, std::string in_info, int xnum, int ynum, 
+	double xmin, double ymin, double dx, double dy)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::init(std::string in_name, std::string in_info, int xbnum, int ybnum, int zbnum, 
+	double xmin, double ymin, double zmin, double xsize, double ysize, double zsize)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::init(std::string in_name, std::string in_info, double lon_min, double lat_min, 
+	double rad_min, double lon_size, double lat_size, double rad_size, int lon_bnum, int lat_bnum, int rad_bnum)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::init(std::string in_name, std::string in_info, double xmin, double ymin, 
+	const gctl::array<double> &xsizes, const gctl::array<double> &ysizes)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::init(std::string in_name, std::string in_info, double xmin, double ymin, 
+	double zmin, const gctl::array<double> &xsizes, const gctl::array<double> &ysizes, 
+	const gctl::array<double> &zsizes)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::init(std::string in_name, std::string in_info, const gctl::array<gctl::vertex2dc> &in_nodes, 
+	const gctl::array<gctl::triangle2d> &in_triangles)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::init(std::string in_name, std::string in_info, const gctl::array<gctl::vertex3dc> &in_nodes, 
+	const gctl::array<gctl::tetrahedron> &in_tets)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::show_mesh_dimension(std::ostream &os) const
+{
+	return;
+}
+
+void gctl::base_mesh::save_gmsh(std::string filename, index_packed_e packed)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
+{
+	if (out_mode == OverWrite) save_gmsh(filename, packed);
+
+	std::ofstream outfile;
+	gctl::open_outfile(outfile, filename, ".msh", std::ios::out|std::ios::app);
+
+	meshdata *curr_data;
+	for (iter = saved_data.begin(); iter != saved_data.end(); ++iter)
+	{
+		curr_data = *iter;
+		if (curr_data->get_dattype() == d_type && d_type == NodeData && 
+			curr_data->get_valtype() == Scalar && curr_data->get_output())
+		{
+			gctl::array<double> *data_ptr = (gctl::array<double>*) curr_data->get_datval_ptr();
+			gctl::save_gmsh_data(outfile, curr_data->get_datname(), *data_ptr, gctl::NodeData, packed);
+		}
+		else if (curr_data->get_dattype() == d_type && d_type == ElemData && 
+			curr_data->get_valtype() == Scalar && curr_data->get_output())
+		{
+			gctl::array<double> *data_ptr = (gctl::array<double>*) curr_data->get_datval_ptr();
+			gctl::save_gmsh_data(outfile, curr_data->get_datname(), *data_ptr, gctl::ElemData, packed);
+		}
+		else if (curr_data->get_dattype() == d_type && d_type == NodeData && 
+			curr_data->get_valtype() == Vector && curr_data->get_output())
+		{
+			gctl::array<gctl::point3dc> *data_ptr = (gctl::array<gctl::point3dc>*) curr_data->get_datval_ptr();
+			gctl::save_gmsh_data(outfile, curr_data->get_datname(), *data_ptr, gctl::NodeData, packed);
+		}
+		else if (curr_data->get_dattype() == d_type && d_type == ElemData && 
+			curr_data->get_valtype() == Vector && curr_data->get_output())
+		{
+			gctl::array<gctl::point3dc> *data_ptr = (gctl::array<gctl::point3dc>*) curr_data->get_datval_ptr();
+			gctl::save_gmsh_data(outfile, curr_data->get_datname(), *data_ptr, gctl::ElemData, packed);
+		}
+	}
+
+	outfile.close();
+	return;
+}
+
+void gctl::base_mesh::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
+{
+	meshdata *curr_data = get_data(datname);
+
+	if (!curr_data->get_output())
+	{
+		throw gctl::runtime_error("[gctl::base_mesh] Output is disabled for the data:" + datname);
+	}
+
+	if (out_mode == OverWrite) save_gmsh(filename, packed);
+
+	std::ofstream outfile;
+	gctl::open_outfile(outfile, filename, ".msh", std::ios::out|std::ios::app);
+
+	if (curr_data->get_valtype() == Scalar)
+	{
+		gctl::array<double> *data_ptr = (gctl::array<double>*) curr_data->get_datval_ptr();
+		if (curr_data->get_dattype() == NodeData)
+		{
+			gctl::save_gmsh_data(outfile, curr_data->get_datname(), *data_ptr, gctl::NodeData, packed);
+		}
+		else if (curr_data->get_dattype() == ElemData)
+		{
+			gctl::save_gmsh_data(outfile, curr_data->get_datname(), *data_ptr, gctl::ElemData, packed);
+		}
+	}
+	else if (curr_data->get_valtype() == Vector)
+	{
+		gctl::array<gctl::point3dc> *data_ptr = (gctl::array<gctl::point3dc>*) curr_data->get_datval_ptr();
+		if (curr_data->get_dattype() == NodeData)
+		{
+			gctl::save_gmsh_data(outfile, curr_data->get_datname(), *data_ptr, gctl::NodeData, packed);
+		}
+		else if (curr_data->get_dattype() == ElemData)
+		{
+			gctl::save_gmsh_data(outfile, curr_data->get_datname(), *data_ptr, gctl::ElemData, packed);
+		}
+	}
+	else if (curr_data->get_valtype() == Tensor)
+	{
+		gctl::array<gctl::tensor> *data_ptr = (gctl::array<gctl::tensor>*) curr_data->get_datval_ptr();
+		if (curr_data->get_dattype() == NodeData)
+		{
+			array<double> v(data_ptr->size());
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[0][0];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Txx", v, gctl::NodeData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[0][1];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Txy", v, gctl::NodeData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[0][2];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Txz", v, gctl::NodeData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[1][1];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Tyy", v, gctl::NodeData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[1][2];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Tyz", v, gctl::NodeData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[2][2];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Tzz", v, gctl::NodeData, packed);
+		}
+		else if (curr_data->get_dattype() == ElemData)
+		{
+			array<double> v(data_ptr->size());
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[0][0];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Txx", v, gctl::ElemData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[0][1];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Txy", v, gctl::ElemData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[0][2];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Txz", v, gctl::ElemData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[1][1];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Tyy", v, gctl::ElemData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[1][2];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Tyz", v, gctl::ElemData, packed);
+
+			for (size_t i = 0; i < v.size(); i++)
+			{
+				v[i] = data_ptr->at(i).val[2][2];
+			}
+			gctl::save_gmsh_data(outfile, curr_data->get_datname() + "_Tzz", v, gctl::ElemData, packed);
+		}
+	}
+
+	outfile.close();
+	return;
+}
+
+void gctl::base_mesh::load_data_cloud(const array<point2dc> &in_posi, const array<double> &in_val, 
+	double search_xlen, double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::load_data_cloud(const array<point3dc> &in_posi, const array<double> &in_val, 
+	double search_xlen, double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::extract_points(std::string datname, const array<point2dc> &in_posi, array<double> &out_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::extract_points(std::string datname, const array<point3dc> &in_posi, array<double> &out_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::extract_profile(std::string datname, const point2dc &start_p, const point2dc &end_p, int size_p, 
+	array<point2dc> &out_posi, array<double> &out_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::extract_profile(std::string datname, const point3dc &start_p, const point3dc &end_p, int size_p, 
+	double dh, array<point3dc> &out_posi, array<double> &out_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, double in_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, point3dc in_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, tensor in_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::purge_data(std::string datname, double in_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::purge_data(std::string datname, point3dc in_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+void gctl::base_mesh::purge_data(std::string datname, tensor in_val)
+{
+	throw runtime_error("[gctl::base_mesh] Invalid mesh type for calling this function.");
+	return;
+}
+
+/**
+ * 以下是类的私有函数，可以简单一些
+ */
+
+void gctl::base_mesh::init(mesh_type_e in_type, mesh_dim_e in_dim, std::string in_name, std::string in_info)
+{
+	if (in_name == "")
+	{
+		throw std::runtime_error("[gctl::base_mesh] The input name is empty.");
+	}
+
+	if (in_info == "")
+	{
+		throw std::runtime_error("[gctl::base_mesh] The input info. is empty.");
+	}
+
+	meshname = in_name;
+	meshinfo = in_info;
+	meshtype = in_type;
+	meshdim  = in_dim;
+	return;
+}
+
+void gctl::base_mesh::load_headinfo(std::ifstream &infile, mesh_type_e expected_type, mesh_dim_e expected_dim)
+{
+	// 读入网格头信息
+	infile.read((char*)&meshtype, sizeof(int));
+	infile.read((char*)&meshdim, sizeof(int));
+	if (meshdim != expected_dim || meshtype != expected_type)
+	{
+		infile.close();
+		throw std::runtime_error("[gctl::base_mesh] Invalid input mesh type.");
+	}
+
+	int info_size;
+	infile.read((char*)&info_size, sizeof(int));
+	meshname.resize(info_size);
+	infile.read((char*)meshname.c_str(), info_size);
+
+	infile.read((char*)&info_size, sizeof(int));
+	meshinfo.resize(info_size);
+	infile.read((char*)meshinfo.c_str(), info_size);
+	return;
+}
+
+void gctl::base_mesh::load_datablock(std::ifstream &infile)
+{
+	meshdata *new_data;
+	int in_num, info_size;
+	mesh_data_type_e in_dattype;
+	mesh_data_value_e in_valtype;
+	std::string in_name;
+
+	infile.read((char*)&in_num, sizeof(int));
+	for (int i = 0; i < in_num; i++)
+	{
+		// 首先读入三个整形和数据名称
+		infile.read((char*)&in_dattype, sizeof(int));
+		infile.read((char*)&in_valtype, sizeof(int));
+		infile.read((char*)&info_size, sizeof(int));
+
+		in_name.resize(info_size);
+		infile.read((char*)in_name.c_str(), info_size);
+
+		new_data = add_data(in_name, in_dattype, true, in_valtype);
+		new_data->load_binary(infile);
+	}
+	return;
+}
+
+void gctl::base_mesh::save_headinfo(std::ofstream &outfile)
+{
+	// 首先输出网格的类型和维度
+	outfile.write((char*)&meshtype, sizeof(int));
+	outfile.write((char*)&meshdim, sizeof(int));
+	// 输出网格名称与信息
+	int info_size = meshname.size();
+	outfile.write((char*)&info_size, sizeof(int));
+	outfile.write((char*)meshname.c_str(), info_size);
+	info_size = meshinfo.size();
+	outfile.write((char*)&info_size, sizeof(int));
+	outfile.write((char*)meshinfo.c_str(), info_size);
+	return;
+}
+
+void gctl::base_mesh::save_datablock(std::ofstream &outfile)
+{
+	// 统计输出的模型数量
+	int out_num = 0;
+	meshdata *curr_data = nullptr;
+	for (iter = saved_data.begin(); iter != saved_data.end(); ++iter)
+	{
+		curr_data = *iter;
+		if (curr_data->get_output())
+		{
+			out_num++;
+		}
+	}
+	outfile.write((char*)&out_num, sizeof(int));
+
+	for (iter = saved_data.begin(); iter != saved_data.end(); ++iter)
+	{
+		curr_data = *iter;
+		if (curr_data->get_output())
+		{
+			curr_data->save_binary(outfile);
+		}
+	}
+	return;
+}

backup/mesh.h

@@ -0,0 +1,579 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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_BASE_MESH_H
+#define _GCTL_BASE_MESH_H
+
+#include "list"
+
+#include "new_meshdata.h"
+//#include "meshdata.h"
+//#include "meshdata_scalar.h"
+//#include "meshdata_vector.h"
+//#include "meshdata_tensor.h"
+
+#include "gctl/io.h"
+#include "gctl/algorithm.h"
+
+namespace gctl
+{
+	enum mesh_type_e
+	{
+		REGULAR_MESH,
+		LINEAR_MESH,
+		TRI_TET_MESH,
+		REGULAR_MESH_SPH,
+		LINEAR_MESH_SPH,
+		TRI_TET_MESH_SPH,
+		REGULAR_GRID,
+	};
+
+	enum mesh_dim_e
+	{
+		MESH_2D,
+		MESH_3D,
+	};
+
+	/**
+	 * @brief      网格对象（纯虚类）。此类为所有具体网格类型的基类
+	 */
+	class base_mesh
+	{
+	public:
+		base_mesh();
+		virtual ~base_mesh();
+		
+		/**
+		 * @brief 清除所有网格数据。
+		 * 
+		 */
+		void clear();
+
+		/**
+		 * @brief 检查网格是否已经初始化。
+		 * 
+		 */
+		bool initiated() const;
+
+		/**
+		 * @brief 检查是否存在对应名称的数据。
+		 * 
+		 * @param datname 数据名
+		 */
+		bool saved(std::string datname) const;
+
+		/**
+		 * @brief 返回对应名称的数据对象指针。
+		 * 
+		 * @param datname 数据名
+		 * @return 数据指针
+		 */
+		meshdata *get_data(std::string datname) const;
+
+		/**
+		 * @brief 返回所有数据指针。
+		 * 
+		 * @param out_list 数据指针列表
+		 */
+		void get_all_data(array<meshdata*> &out_list) const;
+
+		/**
+		 * @brief 返回对应名称的数据对象的数组指针。
+		 * 
+		 * @param datname 数据名
+		 * @return 数组指针（注意需要转换为对应的数组指针类型）
+		 */
+		void *get_datval(std::string datname) const;
+
+		/**
+		 * @brief 删除对应名称的数据对象。
+		 * 
+		 * @param datname 数据名
+		 */
+		void remove_data(std::string datname);
+
+		/**
+		 * @brief 显示网格与数据信息。
+		 * 
+		 * @param os 输出流
+		 */
+		void show_info(std::ostream &os = std::clog) const;
+
+		/**
+		 * @brief 重命名数据对象。
+		 * 
+		 * @param oldname 旧名称
+		 * @param newname 新名称
+		 */
+		void rename_data(std::string oldname, std::string newname);
+
+		/**
+		 * @brief 返回网格类型
+		 * 
+		 * @return 网格类型 
+		 */
+		mesh_type_e get_meshtype() const;
+
+		/**
+		 * @brief 返回网格维度类型
+		 * 
+		 * @return 网格维度
+		 */
+		mesh_dim_e get_meshdim() const;
+
+		/**
+		 * @brief 返回顶点数量
+		 * 
+		 * @return 顶点数量
+		 */
+		int get_nodenum() const;
+
+		/**
+		 * @brief 返回网格单元数量
+		 * 
+		 * @return 单元数量
+		 */
+		int get_elenum() const;
+
+		/**
+		 * @brief 返回数据对象数量
+		 * 
+		 * @return 数据对象数量
+		 */
+		int get_datanum() const;
+
+		/**
+		 * @brief 返回网格名称
+		 * 
+		 * @return 网格名称
+		 */
+		std::string get_meshname() const;
+
+		/**
+		 * @brief 设置网格名称
+		 * 
+		 * @param in_name 网格名称
+		 */
+		void set_meshname(std::string in_name);
+
+		/**
+		 * @brief 返回网格说明信息
+		 * 
+		 * @return 说明信息
+		 */
+		std::string get_meshinfo() const;
+
+		/**
+		 * @brief 设置网格说明信息
+		 * 
+		 * @param in_info 说明信息
+		 */
+		void set_meshinfo(std::string in_info);
+
+		/**
+		 * @brief 添加标量网格数据对象并赋初始值，如果对象已经存在则赋初始值
+		 * 
+		 * @param in_name 数据对象名称
+		 * @param in_type 数据对象类型
+		 * @param if_output 数据对象是否可输出
+		 * @param init_val 初始值
+		 * @return 新建的数据对象指针
+		 */
+		meshdata *add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, double init_val);
+		
+		/**
+		 * @brief 添加适量网格数据对象并赋初始值，如果对象已经存在则赋初始值
+		 * 
+		 * @param in_name 数据对象名称
+		 * @param in_type 数据对象类型
+		 * @param if_output 数据对象是否可输出
+		 * @param init_val 初始值
+		 * @return 新建的数据对象指针
+		 */
+		meshdata *add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, gctl::point3dc init_val);
+		
+		/**
+		 * @brief 添加张量网格数据对象并赋初始值，如果对象已经存在则赋初始值
+		 * 
+		 * @param in_name 数据对象名称
+		 * @param in_type 数据对象类型
+		 * @param if_output 数据对象是否可输出
+		 * @param init_val 初始值
+		 * @return 新建的数据对象指针
+		 */
+		meshdata *add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, gctl::tensor init_val);
+		
+		/**
+		 * @brief 添加张量网格数据对象并初始为0值
+		 * 
+		 * @param in_name 数据对象名称
+		 * @param in_type 数据对象类型
+		 * @param if_output 数据对象是否可输出
+		 * @param init_val 初始值
+		 * @return 新建的数据对象指针
+		 */
+		meshdata *add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, mesh_data_value_e val_type);
+
+		/**
+		 * @brief （基类虚函数）初始化规则网格
+		 * 
+		 * @param in_name 网格名称
+		 * @param in_info 网格说明信息
+		 * @param xnum 网格x轴数量
+		 * @param ynum 网格y轴数量
+		 * @param xmin 网格x轴最小值
+		 * @param ymin 网格y轴最小值
+		 * @param dx 网格x轴间隔
+		 * @param dy 网格y轴间隔
+		 */
+		virtual void init(std::string in_name, std::string in_info, int xnum, int ynum, 
+			double xmin, double ymin, double dx, double dy);
+
+		/**
+		 * @brief 
+		 * 
+		 * @param in_name 
+		 * @param in_info 
+		 * @param xbnum 
+		 * @param ybnum 
+		 * @param zbnum 
+		 * @param xmin 
+		 * @param ymin 
+		 * @param zmin 
+		 * @param xsize 
+		 * @param ysize 
+		 * @param zsize 
+		 */
+		virtual void init(std::string in_name, std::string in_info, int xbnum, int ybnum, int zbnum, 
+			double xmin, double ymin, double zmin, double xsize, double ysize, double zsize);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param in_name 
+		 * @param in_info 
+		 * @param lon_min 
+		 * @param lat_min 
+		 * @param rad_min 
+		 * @param lon_size 
+		 * @param lat_size 
+		 * @param rad_size 
+		 * @param lon_bnum 
+		 * @param lat_bnum 
+		 * @param rad_bnum 
+		 */
+		virtual void init(std::string in_name, std::string in_info, double lon_min, double lat_min, 
+			double rad_min, double lon_size, double lat_size, double rad_size, int lon_bnum, int lat_bnum, int rad_bnum);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param in_name 
+		 * @param in_info 
+		 * @param xmin 
+		 * @param ymin 
+		 * @param xsizes 
+		 * @param ysizes 
+		 */
+		virtual void init(std::string in_name, std::string in_info, double xmin, double ymin, 
+			const gctl::array<double> &xsizes, const gctl::array<double> &ysizes);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param in_name 
+		 * @param in_info 
+		 * @param xmin 
+		 * @param ymin 
+		 * @param zmin 
+		 * @param xsizes 
+		 * @param ysizes 
+		 * @param zsizes 
+		 */
+		virtual void init(std::string in_name, std::string in_info, double xmin, double ymin, 
+			double zmin, const gctl::array<double> &xsizes, const gctl::array<double> &ysizes, 
+			const gctl::array<double> &zsizes);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param in_name 
+		 * @param in_info 
+		 * @param in_nodes 
+		 * @param in_triangles 
+		 */
+		virtual void init(std::string in_name, std::string in_info, const gctl::array<gctl::vertex2dc> &in_nodes, 
+			const gctl::array<gctl::triangle2d> &in_triangles);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param in_name 
+		 * @param in_info 
+		 * @param in_nodes 
+		 * @param in_tets 
+		 */
+		virtual void init(std::string in_name, std::string in_info, const gctl::array<gctl::vertex3dc> &in_nodes, 
+			const gctl::array<gctl::tetrahedron> &in_tets);
+
+		/**
+		 * @brief 显示网格的维度信息
+		 * 
+		 */
+		virtual void show_mesh_dimension(std::ostream &os) const;
+
+		/**
+		 * @brief （基类纯虚函数）读入二进制网格文件
+		 * 
+		 * @param filename 文件名
+		 */
+		virtual void load_binary(std::string filename) = 0;
+
+		/**
+		 * @brief （基类纯虚函数）写入二进制网格文件
+		 * 
+		 * @param filename 文件名
+		 */
+		virtual void save_binary(std::string filename) = 0;
+
+		/**
+		 * @brief （基类虚函数）保存网格到Gmsh文件。
+		 * 
+		 * @param filename 文件名
+		 * @param packed 索引是否为紧凑（从0开始）
+		 */
+		virtual void save_gmsh(std::string filename, index_packed_e packed = Packed);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param filename 
+		 * @param d_type 
+		 * @param out_mode 
+		 * @param packed 
+		 */
+		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param filename 
+		 * @param datname 
+		 * @param out_mode 
+		 * @param packed 
+		 */
+		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
+
+		/**
+		 * @brief      将输入的自由点位上的浮点数据插值为规则网络的顶点或者块体数据
+		 *
+		 * @param      in_posi      输入的自由点位的坐标
+		 * @param      in_val       输入的自由点位的数值
+		 * @param[in]  search_xlen  插值搜索区域的x半径
+		 * @param[in]  search_ylen  插值搜索区域的y半径
+		 * @param[in]  search_deg   插值搜索区域的x半径绕x轴正方向逆时针旋转的角度（弧度）
+		 * @param[in]  datname      网格后的数据名称
+		 * @param[in]  d_type       网格后的数据类型
+		 */
+		virtual void load_data_cloud(const array<point2dc> &in_posi, const array<double> &in_val, double search_xlen, 
+			double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type = NodeData);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param in_posi 
+		 * @param in_val 
+		 * @param search_xlen 
+		 * @param search_ylen 
+		 * @param search_deg 
+		 * @param datname 
+		 * @param d_type 
+		 */
+		virtual void load_data_cloud(const array<point3dc> &in_posi, const array<double> &in_val, double search_xlen, 
+			double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type = NodeData);
+
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param in_posi 
+		 * @param out_val 
+		 */
+		virtual void extract_points(std::string datname, const array<point2dc> &in_posi, array<double> &out_val);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param in_posi 
+		 * @param out_val 
+		 */
+		virtual void extract_points(std::string datname, const array<point3dc> &in_posi, array<double> &out_val);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param start_p 
+		 * @param end_p 
+		 * @param size_p 
+		 * @param out_posi 
+		 * @param out_val 
+		 */
+		virtual void extract_profile(std::string datname, const point2dc &start_p, const point2dc &end_p, int size_p, 
+			array<point2dc> &out_posi, array<double> &out_val);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param start_p 
+		 * @param end_p 
+		 * @param size_p 
+		 * @param dh 
+		 * @param out_posi 
+		 * @param out_val 
+		 */
+		virtual void extract_profile(std::string datname, const point3dc &start_p, const point3dc &end_p, int size_p, 
+			double dh, array<point3dc> &out_posi, array<double> &out_val);
+
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param p_type 
+		 * @param v_type 
+		 * @param para_str 
+		 * @param in_val 
+		 */
+		virtual void edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, double in_val);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param p_type 
+		 * @param v_type 
+		 * @param para_str 
+		 * @param in_val 
+		 */
+		virtual void edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, point3dc in_val);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param p_type 
+		 * @param v_type 
+		 * @param para_str 
+		 * @param in_val 
+		 */
+		virtual void edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, tensor in_val);
+
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param in_val 
+		 */
+		virtual void purge_data(std::string datname, double in_val);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param in_val 
+		 */
+		virtual void purge_data(std::string datname, point3dc in_val);
+		
+		/**
+		 * @brief 
+		 * 
+		 * @param datname 
+		 * @param in_val 
+		 */
+		virtual void purge_data(std::string datname, tensor in_val);
+
+	protected:
+		mesh_type_e meshtype;
+		mesh_dim_e  meshdim;
+		std::string meshname;
+		std::string meshinfo;
+
+		int node_num, ele_num;
+		bool initialized;
+
+		std::list<meshdata*> saved_data;
+		std::list<meshdata*>::iterator iter;
+
+		/**
+		 * 以下为类的私有函数 只能被类的公共函数调用
+		 */
+		
+		/**
+		 * @brief 初始化网格
+		 * 
+		 * @param in_type 
+		 * @param in_dim 
+		 * @param in_name 
+		 * @param in_info 
+		 */
+		void init(mesh_type_e in_type, mesh_dim_e in_dim, std::string in_name, std::string in_info);
+
+		/**
+		 * @brief 读入二进制网格头信息
+		 * 
+		 * @param infile 输入流
+		 * @param expected_type 预期获取的网格类型 
+		 * @param expected_dim 预期获得的网格维度类型
+		 */
+		void load_headinfo(std::ifstream &infile, mesh_type_e expected_type, mesh_dim_e expected_dim);
+		
+		/**
+		 * @brief 读入二进制数据块
+		 * 
+		 * @param infile 输入流
+		 */
+		void load_datablock(std::ifstream &infile);
+		
+		/**
+		 * @brief 写入二进制网格头信息
+		 * 
+		 * @param outfile 输出流
+		 */
+		void save_headinfo(std::ofstream &outfile);
+		
+		/**
+		 * @brief 写入二进制数据块
+		 * 
+		 * @param outfile 输出流
+		 */
+		void save_datablock(std::ofstream &outfile);
+	};
+}
+
+#endif //_GCTL_BASE_MESH_H

backup/meshdata.cpp

@@ -0,0 +1,94 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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 "meshdata.h"
+
+gctl::meshdata::meshdata(std::string in_name, mesh_data_type_e in_type, bool if_output)
+{
+	if (in_name == "")
+	{
+		throw std::runtime_error("[gctl::meshdata] The input name is empty.");
+	}
+
+	datname = in_name;
+	dattype = in_type;
+	output_ok = if_output;
+}
+
+gctl::meshdata::~meshdata(){}
+
+void gctl::meshdata::set_datname(std::string in_name)
+{
+	if (in_name == "")
+	{
+		throw std::runtime_error("[gctl::meshdata] The input name is empty.");
+	}
+
+	datname = in_name;
+	return;
+}
+
+std::string gctl::meshdata::get_datname()
+{
+	return datname;
+}
+
+gctl::mesh_data_type_e gctl::meshdata::get_dattype()
+{
+	return dattype;
+}
+
+gctl::mesh_data_value_e gctl::meshdata::get_valtype()
+{
+	return valtype;
+}
+
+void gctl::meshdata::set_output(bool if_output)
+{
+	output_ok = if_output;
+	return;
+}
+
+bool gctl::meshdata::get_output()
+{
+	return output_ok;
+}
+
+void gctl::meshdata::show_info(std::ostream &os)
+{
+	os << "Data: " << datname << " | ";
+	if (dattype == NodeData) os << "Type: Node data | ";
+	if (dattype == ElemData) os << "Type: Element data | ";
+	if (dattype == ElemData2D) os << "Type: 2D element data | ";
+	if (dattype == ElemData3D) os << "Type: 3D element data | ";
+	if (valtype == Scalar) os << "Value: Scalar | ";
+	if (valtype == Vector) os << "Value: Vector | ";
+	if (valtype == Tensor) os << "Value: Tensor | ";
+	if (output_ok) os << "Output: Yes" << std::endl;
+	else os<< "Output: No" << std::endl;
+	return;
+}

backup/meshdata.h

@@ -0,0 +1,158 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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_MESHDATA_H
+#define _GCTL_MESHDATA_H
+
+#include "gctl_mesh_config.h"
+#include "gctl/core.h"
+#include "gctl/geometry.h"
+
+namespace gctl
+{
+	/**
+	 * @brief      网格数据的数值类型
+	 */
+	enum mesh_data_value_e
+	{
+		Scalar, ///< 标量数据
+		Vector, ///< 矢量数据
+		Tensor, ///< 张量数据
+	};
+
+	/**
+	 * @brief      网格文件中的数据对象（纯虚类）。此类为具体类型的数据对象的基类
+	 */
+	class meshdata
+	{
+	protected:
+		std::string datname; // 数据的名称
+		mesh_data_type_e dattype; // 数据的赋值属性 顶点或是元素（设置后不可更改）
+		mesh_data_value_e valtype; // 数据的类型（设置后不可更改）
+		bool output_ok; // 是否可输出数据
+
+	public:
+		// 构造函数
+		meshdata(std::string in_name, mesh_data_type_e in_type, bool if_output);
+
+		// 析构函数
+		virtual ~meshdata();
+
+		/**
+		 * @brief      设置数据名称
+		 *
+		 * @param[in]  in_name  名称字符串
+		 */
+		void set_datname(std::string in_name);
+
+		/**
+		 * @brief      返回数据的名称
+		 *
+		 * @return     名称字符串
+		 */
+		std::string get_datname();
+
+		/**
+		 * @brief      返回数据的赋值类型
+		 *
+		 * @return     数据赋值类型的枚举
+		 */
+		mesh_data_type_e get_dattype();
+
+		/**
+		 * @brief      返回数据值的类型
+		 *
+		 * @return     数据值类型的枚举
+		 */
+		mesh_data_value_e get_valtype();
+
+		/**
+		 * @brief      设置数据对象输出状态的设置情况
+		 *
+		 * @param[in]  if_output  是否可输出此对象
+		 */
+		void set_output(bool if_output);
+
+		/**
+		 * @brief      返回数据对象输出状态的设置情况
+		 *
+		 * @return     是否可输出此对象
+		 */
+		bool get_output();
+
+		/**
+		 * @brief      显示数据对象的头信息
+		 */
+		void show_info(std::ostream &os = std::clog);
+
+		/**
+		 * @brief     显示数据统计参数
+		 */
+		virtual void show_stats(std::ostream &os = std::clog) = 0;
+
+		/**
+		 * @brief      返回数据块的指针 所有子类的数据块均为gctl::array类型，这是一个通用的接口。
+		 * 
+		 * @note       注意因为具体的返回类型未知 因此需在调用后按实际数据类型转换指针
+		 *
+		 * @return     数据块的指针
+		 */
+		virtual void *get_datval_ptr() = 0;
+
+		/**
+		 * @brief      载入二进制文件
+		 * 
+		 * @warning    不要直接调用此函数，应该在网格类型的相应函数中调用
+		 *
+		 * @param      infile  输入文件的流（以二进制方式打开）
+		 */
+		virtual void load_binary(std::ifstream &infile) = 0;
+
+		/**
+		 * @brief      保存二进制文件
+		 * 
+		 * @warning    不要直接调用此函数，应该在网格类型的相应函数中调用
+		 *
+		 * @param      outfile  输出的文件流（以二进制方式打开）
+		 */
+		virtual void save_binary(std::ofstream &outfile) = 0;
+
+		/**
+		 * @brief      销毁指针的对象 可销毁此基类或子类指针指向的对象
+		 *
+		 * @param      obj_ptr  网格数据指针
+		 */
+		static void destroy(meshdata *obj_ptr)
+		{
+			delete obj_ptr;
+			obj_ptr = nullptr;
+			return;
+		}
+	};
+}
+
+#endif //_GCTL_MESHDATA_H

doc/linear_mesh_2d.md

@@ -0,0 +1,109 @@
+# GCTL Linear Mesh 2D 文档
+
+## 简介
+
+`linear_mesh_2d.h` 定义了GCTL网格库中的二维线性网格类 `linear_mesh_2d`，用于创建和管理二维非结构化线性网格。该类支持不规则分布的网格点，适用于复杂几何形状的离散化。
+
+## 类继承
+
+`linear_mesh_2d` 继承自 `base_mesh` 基类，实现了所有虚函数，并添加了二维线性网格特有的功能。
+
+## 主要API
+
+### 初始化
+```cpp
+// 初始化二维线性网格
+void init(std::string name,                // 网格名称
+          std::string info,                // 网格信息
+          const array<point2dc> &nodes,    // 节点坐标数组
+          const array<int2> &elements);    // 单元连接关系数组
+```
+
+### 网格信息
+```cpp
+// 获取节点坐标
+void get_node_coord(int node_idx,          // 节点索引
+                   double &x, double &y) const;
+
+// 获取单元节点索引
+void get_elem_nodes(int elem_idx,          // 单元索引
+                   int &n1, int &n2) const;
+
+// 获取单元中心坐标
+void get_elem_center(int elem_idx,         // 单元索引
+                    double &x, double &y) const;
+
+// 获取单元长度
+double get_elem_length(int elem_idx) const;
+```
+
+### 网格查询
+```cpp
+// 查找包含指定点的单元
+bool find_element(double x, double y,      // 目标点坐标
+                 int &elem_idx) const;     // 返回单元索引
+
+// 获取节点相邻的单元
+array<int> get_node_elements(int node_idx) const;
+
+// 获取单元相邻的单元
+array<int> get_elem_neighbors(int elem_idx) const;
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) override;
+
+// 保存为二进制格式
+virtual void save_binary(std::string filename) override;
+
+// 读取二进制格式
+virtual void load_binary(std::string filename) override;
+```
+
+## 使用示例
+
+```cpp
+// 创建节点坐标数组
+array<point2dc> nodes;
+nodes.push_back(point2dc(0.0, 0.0));
+nodes.push_back(point2dc(1.0, 0.0));
+nodes.push_back(point2dc(1.0, 1.0));
+nodes.push_back(point2dc(0.0, 1.0));
+
+// 创建单元连接关系
+array<int2> elements;
+elements.push_back(int2(0, 1));  // 第一条边
+elements.push_back(int2(1, 2));  // 第二条边
+elements.push_back(int2(2, 3));  // 第三条边
+elements.push_back(int2(3, 0));  // 第四条边
+
+// 创建线性网格
+gctl::linear_mesh_2d mesh2d;
+mesh2d.init("boundary", "2D boundary mesh", nodes, elements);
+
+// 获取单元信息
+double length = mesh2d.get_elem_length(0);  // 获取第一条边的长度
+
+// 查找包含特定点的单元
+int elem_idx;
+if (mesh2d.find_element(0.5, 0.5, elem_idx)) {
+    // 找到包含点(0.5,0.5)的单元
+}
+
+// 添加数据并保存
+mesh2d.add_data(gctl::NodeData, gctl::Scalar, "displacement", 0.0);
+mesh2d.save_gmsh("boundary.msh", gctl::Packed);
+```
+
+## 注意事项
+
+1. 线性网格主要用于表示一维结构（如边界）在二维空间中的分布
+2. 每个单元由两个节点定义，表示一条线段
+3. 节点和单元的编号从0开始
+4. 支持不规则分布的节点，适合复杂几何形状
+5. 提供了网格拓扑关系查询功能（节点相邻单元、单元相邻单元等）
+6. 支持Gmsh和二进制格式的文件操作
+7. 适用于边界元素法、有限元法等数值计算 

doc/linear_mesh_3d.md

@@ -0,0 +1,121 @@
+# GCTL Linear Mesh 3D 文档
+
+## 简介
+
+`linear_mesh_3d.h` 定义了GCTL网格库中的三维线性网格类 `linear_mesh_3d`，用于创建和管理三维空间中的非结构化线性网格。该类支持不规则分布的网格点，适用于复杂三维几何结构中的线性特征离散化。
+
+## 类继承
+
+`linear_mesh_3d` 继承自 `base_mesh` 基类，实现了所有虚函数，并添加了三维线性网格特有的功能。
+
+## 主要API
+
+### 初始化
+```cpp
+// 初始化三维线性网格
+void init(std::string name,                // 网格名称
+          std::string info,                // 网格信息
+          const array<point3dc> &nodes,    // 节点坐标数组
+          const array<int2> &elements);    // 单元连接关系数组
+```
+
+### 网格信息
+```cpp
+// 获取节点坐标
+void get_node_coord(int node_idx,          // 节点索引
+                   double &x, double &y, double &z) const;
+
+// 获取单元节点索引
+void get_elem_nodes(int elem_idx,          // 单元索引
+                   int &n1, int &n2) const;
+
+// 获取单元中心坐标
+void get_elem_center(int elem_idx,         // 单元索引
+                    double &x, double &y, double &z) const;
+
+// 获取单元长度
+double get_elem_length(int elem_idx) const;
+
+// 获取单元方向向量
+void get_elem_direction(int elem_idx,      // 单元索引
+                       double &dx, double &dy, double &dz) const;
+```
+
+### 网格查询
+```cpp
+// 查找最近的单元
+bool find_nearest_element(double x, double y, double z,  // 目标点坐标
+                         int &elem_idx,                  // 返回单元索引
+                         double &distance) const;        // 返回最短距离
+
+// 获取节点相邻的单元
+array<int> get_node_elements(int node_idx) const;
+
+// 获取单元相邻的单元
+array<int> get_elem_neighbors(int elem_idx) const;
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) override;
+
+// 保存为二进制格式
+virtual void save_binary(std::string filename) override;
+
+// 读取二进制格式
+virtual void load_binary(std::string filename) override;
+```
+
+## 使用示例
+
+```cpp
+// 创建节点坐标数组
+array<point3dc> nodes;
+nodes.push_back(point3dc(0.0, 0.0, 0.0));
+nodes.push_back(point3dc(1.0, 0.0, 0.0));
+nodes.push_back(point3dc(1.0, 1.0, 1.0));
+nodes.push_back(point3dc(0.0, 1.0, 1.0));
+
+// 创建单元连接关系
+array<int2> elements;
+elements.push_back(int2(0, 1));  // 第一条边
+elements.push_back(int2(1, 2));  // 第二条边
+elements.push_back(int2(2, 3));  // 第三条边
+
+// 创建线性网格
+gctl::linear_mesh_3d mesh3d;
+mesh3d.init("curve3d", "3D curve mesh", nodes, elements);
+
+// 获取单元信息
+double length = mesh3d.get_elem_length(0);  // 获取第一条边的长度
+
+// 获取单元方向
+double dx, dy, dz;
+mesh3d.get_elem_direction(0, dx, dy, dz);  // 获取第一条边的方向向量
+
+// 查找最近的单元
+int elem_idx;
+double distance;
+if (mesh3d.find_nearest_element(0.5, 0.5, 0.5, elem_idx, distance)) {
+    // 找到距离点(0.5,0.5,0.5)最近的单元
+}
+
+// 添加数据并保存
+mesh3d.add_data(gctl::NodeData, gctl::Vector, "displacement", 
+                array<point3dc>(nodes.size(), point3dc(0,0,0)));
+mesh3d.save_gmsh("curve3d.msh", gctl::Packed);
+```
+
+## 注意事项
+
+1. 线性网格主要用于表示一维结构（如曲线、管线）在三维空间中的分布
+2. 每个单元由两个节点定义，表示一条空间线段
+3. 节点和单元的编号从0开始
+4. 支持不规则分布的节点，适合复杂三维曲线的离散化
+5. 提供了网格拓扑关系查询功能
+6. 支持计算单元的方向向量，便于处理定向问题
+7. 查找最近单元功能适用于路径规划和碰撞检测
+8. 支持Gmsh和二进制格式的文件操作
+9. 适用于输电线路、管道系统、裂缝网络等模拟 

doc/mesh.md

@@ -0,0 +1,116 @@
+# GCTL Mesh 基类文档
+
+## 简介
+
+`mesh.h` 定义了GCTL网格库的基础类 `base_mesh`，这是所有具体网格类型的基类。它提供了网格操作的基本接口和功能实现。
+
+## 类型定义
+
+### 网格类型 (mesh_type_e)
+```cpp
+enum mesh_type_e {
+    UNDEFINED,           // 未定义网格
+    REGULAR_MESH,        // 规则网格
+    LINEAR_MESH,         // 线性网格
+    TRI_TET_MESH,       // 三角形/四面体网格
+    REGULAR_MESH_SPH,    // 球面规则网格
+    LINEAR_MESH_SPH,     // 球面线性网格
+    TRI_TET_MESH_SPH,    // 球面三角形/四面体网格
+    REGULAR_GRID,        // 规则网格
+};
+```
+
+### 网格维度 (mesh_dim_e)
+```cpp
+enum mesh_dim_e {
+    MESH_0D,    // 0维网格
+    MESH_2D,    // 2维网格
+    MESH_3D     // 3维网格
+};
+```
+
+## 主要API
+
+### 基本操作
+```cpp
+void clear();                    // 清除所有网格数据
+bool initiated() const;          // 检查网格是否已初始化
+void show_info();               // 显示网格与数据信息
+```
+
+### 网格信息获取
+```cpp
+mesh_type_e get_meshtype() const;    // 获取网格类型
+mesh_dim_e get_meshdim() const;      // 获取网格维度
+int get_nodenum() const;             // 获取顶点数量
+int get_elenum() const;              // 获取网格单元数量
+int get_datanum() const;             // 获取数据对象数量
+```
+
+### 网格属性设置
+```cpp
+std::string get_meshname() const;             // 获取网格名称
+void set_meshname(std::string in_name);       // 设置网格名称
+std::string get_meshinfo() const;             // 获取网格说明信息
+void set_meshinfo(std::string in_info);       // 设置网格说明信息
+```
+
+### 数据操作
+```cpp
+// 添加网格数据（标量初始值）
+meshdata &add_data(mesh_data_type_e in_loctype,    // 数据位置类型
+                  mesh_data_value_e in_valtype,     // 数据值类型
+                  std::string name,                 // 数据名称
+                  double init_val,                  // 初始值
+                  bool if_output = true);           // 是否可输出
+
+// 添加网格数据（数组初始值）
+meshdata &add_data(mesh_data_type_e in_loctype,    // 数据位置类型
+                  std::string name,                 // 数据名称
+                  const array<double> &init_arr,    // 初始值数组
+                  bool if_output = true);           // 是否可输出
+
+// 获取数据对象
+meshdata &get_data(std::string datname);
+
+// 删除数据对象
+void remove_data(std::string datname);
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式（含数据）
+void save_gmsh_withdata(std::string filename,           // 文件名
+                       output_type_e out_mode,          // 输出模式
+                       index_packed_e packed = Packed); // 索引打包方式
+
+// 保存为Gmsh格式（指定数据）
+void save_gmsh_withdata(std::string filename,           // 文件名
+                       std::string datname,             // 数据名称
+                       output_type_e out_mode,          // 输出模式
+                       index_packed_e packed = Packed); // 索引打包方式
+
+// 纯虚函数 - 需要由派生类实现
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) = 0;
+virtual void load_binary(std::string filename) = 0;
+virtual void save_binary(std::string filename) = 0;
+```
+
+## 使用示例
+
+```cpp
+// 创建网格对象（以regular_grid为例）
+gctl::regular_grid rgd;
+rgd.init("grid-1", "test grid", 4, 4, 0.0, 0.0, 1.0, 1.0);
+
+// 添加数据
+rgd.add_data(gctl::NodeData, gctl::Scalar, "temperature", 25.0);
+rgd.add_data(gctl::ElemData, gctl::Vector, "velocity", velocity_array);
+
+// 显示网格信息
+rgd.show_info();
+
+// 保存网格
+rgd.save_gmsh_withdata("output.msh", gctl::OverWrite);
+``` 

doc/meshdata.md

@@ -0,0 +1,116 @@
+# GCTL MeshData 文档
+
+## 简介
+
+`meshdata.h` 定义了GCTL网格库中的数据对象类 `meshdata`，用于在网格上存储和管理各种类型的数据（标量、矢量和张量）。
+
+## 类型定义
+
+### 数据值类型 (mesh_data_value_e)
+```cpp
+enum mesh_data_value_e {
+    Scalar,     // 标量数据
+    Vector,     // 矢量数据
+    Tensor      // 张量数据
+};
+```
+
+## 类成员
+
+### 属性
+```cpp
+std::string name_;           // 数据的名称
+mesh_data_type_e loctype_;   // 数据的赋值位置属性（顶点或元素）
+mesh_data_value_e valtype_;  // 数据的类型（标量/矢量/张量）
+bool output_ok_;             // 是否可输出数据
+array<double> datval_;       // 数据值数组
+double nan_val_;             // 无效数据的标记值
+```
+
+## 主要API
+
+### 构造和初始化
+```cpp
+// 默认构造函数
+meshdata();
+
+// 带参数构造函数
+meshdata(mesh_data_type_e in_loctype,     // 数据位置类型
+         mesh_data_value_e in_valtype,     // 数据值类型
+         size_t size,                      // 数据大小
+         std::string name = "untitled",    // 数据名称
+         bool if_output = true,            // 是否可输出
+         double nan_val = GCTL_BDL_MAX);   // 无效值标记
+
+// 创建数据对象
+void create(mesh_data_type_e in_loctype,   // 数据位置类型
+           mesh_data_value_e in_valtype,   // 数据值类型
+           size_t size,                    // 数据大小
+           std::string name = "untitled",  // 数据名称
+           bool if_output = true,          // 是否可输出
+           double nan_val = GCTL_BDL_MAX); // 无效值标记
+```
+
+### 数据操作
+```cpp
+// 清空数据对象
+void clear();
+
+// 导出矢量数据
+array<point3dc> export_vector() const;
+
+// 导出张量数据
+array<tensor> export_tensor() const;
+```
+
+### 信息显示
+```cpp
+// 显示数据对象的头信息
+void show_info(std::ostream &os = std::clog) const;
+
+// 显示数据统计参数
+void show_stats(std::ostream &os = std::clog) const;
+```
+
+### 文件操作
+```cpp
+// 载入二进制文件
+void load_binary(std::ifstream &infile);
+
+// 保存二进制文件
+void save_binary(std::ofstream &outfile);
+```
+
+## 使用示例
+
+```cpp
+// 创建数据对象
+gctl::meshdata data(gctl::NodeData,        // 节点数据
+                   gctl::Scalar,           // 标量类型
+                   100,                    // 数据大小
+                   "temperature",          // 数据名称
+                   true);                  // 可输出
+
+// 显示数据信息
+data.show_info();
+
+// 显示数据统计信息
+data.show_stats();
+
+// 导出矢量数据
+if (data.valtype_ == gctl::Vector) {
+    array<point3dc> vectors = data.export_vector();
+}
+
+// 导出张量数据
+if (data.valtype_ == gctl::Tensor) {
+    array<tensor> tensors = data.export_tensor();
+}
+```
+
+## 注意事项
+
+1. 数据对象支持三种类型：标量、矢量和张量
+2. 数据可以定义在网格的顶点或元素上
+3. 数据导出时会根据数据类型自动转换为对应的格式
+4. 二进制文件操作通常由网格类调用，不建议直接使用 

doc/regular_grid.md

@@ -0,0 +1,114 @@
+# GCTL Regular Grid 文档
+
+## 简介
+
+`regular_grid.h` 定义了GCTL网格库中的规则网格类 `regular_grid`，用于创建和管理规则网格。规则网格是最基本的网格类型，具有规则的网格点分布和单元结构。
+
+## 类继承
+
+`regular_grid` 继承自 `base_mesh` 基类，实现了所有虚函数，并添加了规则网格特有的功能。
+
+## 主要API
+
+### 初始化
+```cpp
+// 初始化二维规则网格
+void init(std::string name,          // 网格名称
+          std::string info,          // 网格信息
+          int nx, int ny,            // x和y方向的网格数
+          double x0, double y0,      // 起始点坐标
+          double x1, double y1);     // 终止点坐标
+
+// 初始化三维规则网格
+void init(std::string name,          // 网格名称
+          std::string info,          // 网格信息
+          int nx, int ny, int nz,    // x、y和z方向的网格数
+          double x0, double y0, double z0,  // 起始点坐标
+          double x1, double y1, double z1); // 终止点坐标
+```
+
+### 网格信息
+```cpp
+// 获取网格维度
+int get_dimension() const;
+
+// 获取网格点数
+void get_grid_num(int &nx, int &ny) const;          // 2D
+void get_grid_num(int &nx, int &ny, int &nz) const; // 3D
+
+// 获取网格范围
+void get_grid_range(double &x0, double &y0,         // 2D起点
+                   double &x1, double &y1) const;    // 2D终点
+void get_grid_range(double &x0, double &y0, double &z0,  // 3D起点
+                   double &x1, double &y1, double &z1) const; // 3D终点
+```
+
+### 网格操作
+```cpp
+// 获取指定位置的网格点坐标
+void get_node_coord(int i, int j,             // 2D索引
+                   double &x, double &y) const;
+void get_node_coord(int i, int j, int k,      // 3D索引
+                   double &x, double &y, double &z) const;
+
+// 获取指定位置的单元中心坐标
+void get_elem_center(int i, int j,            // 2D索引
+                    double &x, double &y) const;
+void get_elem_center(int i, int j, int k,     // 3D索引
+                    double &x, double &y, double &z) const;
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) override;
+
+// 保存为二进制格式
+virtual void save_binary(std::string filename) override;
+
+// 读取二进制格式
+virtual void load_binary(std::string filename) override;
+
+// NetCDF格式操作
+void load_netcdf_grid(std::string filename,
+                     mesh_data_type_e loctype,
+                     std::string xname,
+                     std::string yname);
+
+void save_netcdf_grid(std::string filename,
+                     std::string title = "untitled");
+```
+
+## 使用示例
+
+```cpp
+// 创建2D规则网格
+gctl::regular_grid grid2d;
+grid2d.init("grid2d", "2D regular grid",
+            100, 100,           // 100x100网格
+            0.0, 0.0,          // 起点(0,0)
+            1.0, 1.0);         // 终点(1,1)
+
+// 创建3D规则网格
+gctl::regular_grid grid3d;
+grid3d.init("grid3d", "3D regular grid",
+            50, 50, 50,        // 50x50x50网格
+            0.0, 0.0, 0.0,     // 起点(0,0,0)
+            1.0, 1.0, 1.0);    // 终点(1,1,1)
+
+// 获取网格信息
+int nx, ny, nz;
+grid3d.get_grid_num(nx, ny, nz);
+
+// 添加数据并保存
+grid3d.add_data(gctl::NodeData, gctl::Scalar, "temperature", 25.0);
+grid3d.save_gmsh("output.msh", gctl::Packed);
+```
+
+## 注意事项
+
+1. 规则网格支持2D和3D两种维度
+2. 网格点和单元的编号是规则的，可以通过索引直接访问
+3. 支持多种文件格式的导入导出（Gmsh、二进制、NetCDF）
+4. 网格范围由起点和终点坐标定义，网格间距自动计算 

doc/regular_mesh_2d.md

@@ -0,0 +1,98 @@
+# GCTL Regular Mesh 2D 文档
+
+## 简介
+
+`regular_mesh_2d.h` 定义了GCTL网格库中的二维规则网格类 `regular_mesh_2d`，用于创建和管理二维结构化网格。与`regular_grid`相比，该类更专注于二维网格的处理，提供了更多针对二维网格的专门操作。
+
+## 类继承
+
+`regular_mesh_2d` 继承自 `base_mesh` 基类，实现了所有虚函数，并添加了二维规则网格特有的功能。
+
+## 主要API
+
+### 初始化
+```cpp
+// 初始化二维规则网格
+void init(std::string name,          // 网格名称
+          std::string info,          // 网格信息
+          int nx, int ny,            // x和y方向的网格数
+          double x0, double y0,      // 起始点坐标
+          double x1, double y1,      // 终止点坐标
+          double dx = 0.0,           // x方向网格间距（0表示自动计算）
+          double dy = 0.0);          // y方向网格间距（0表示自动计算）
+```
+
+### 网格信息
+```cpp
+// 获取网格点数
+void get_mesh_num(int &nx, int &ny) const;
+
+// 获取网格范围
+void get_mesh_range(double &x0, double &y0,    // 起点坐标
+                   double &x1, double &y1) const; // 终点坐标
+
+// 获取网格间距
+void get_mesh_step(double &dx, double &dy) const;
+```
+
+### 网格操作
+```cpp
+// 获取指定位置的网格点坐标
+void get_node_coord(int i, int j,          // 节点索引
+                   double &x, double &y) const;
+
+// 获取指定位置的单元中心坐标
+void get_elem_center(int i, int j,         // 单元索引
+                    double &x, double &y) const;
+
+// 获取指定坐标所在的网格单元索引
+bool get_elem_index(double x, double y,    // 坐标
+                   int &i, int &j) const;  // 返回单元索引
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) override;
+
+// 保存为二进制格式
+virtual void save_binary(std::string filename) override;
+
+// 读取二进制格式
+virtual void load_binary(std::string filename) override;
+```
+
+## 使用示例
+
+```cpp
+// 创建二维规则网格
+gctl::regular_mesh_2d mesh2d;
+mesh2d.init("mesh2d", "2D regular mesh",
+            100, 100,           // 100x100网格
+            0.0, 0.0,          // 起点(0,0)
+            1.0, 1.0);         // 终点(1,1)
+
+// 获取网格信息
+int nx, ny;
+mesh2d.get_mesh_num(nx, ny);
+
+double dx, dy;
+mesh2d.get_mesh_step(dx, dy);
+
+// 获取特定位置的坐标
+double x, y;
+mesh2d.get_node_coord(50, 50, x, y);  // 获取中心点坐标
+
+// 添加数据并保存
+mesh2d.add_data(gctl::NodeData, gctl::Scalar, "elevation", 0.0);
+mesh2d.save_gmsh("terrain.msh", gctl::Packed);
+```
+
+## 注意事项
+
+1. 网格间距可以自动计算或手动指定
+2. 提供了坐标与网格索引之间的转换功能
+3. 所有索引都是从0开始计数
+4. 网格点和单元的编号是规则的，遵循行优先顺序
+5. 支持Gmsh和二进制格式的文件操作 

doc/regular_mesh_3d.md

@@ -0,0 +1,106 @@
+# GCTL Regular Mesh 3D 文档
+
+## 简介
+
+`regular_mesh_3d.h` 定义了GCTL网格库中的三维规则网格类 `regular_mesh_3d`，用于创建和管理三维结构化网格。该类专门处理三维空间中的规则网格，提供了完整的三维网格操作功能。
+
+## 类继承
+
+`regular_mesh_3d` 继承自 `base_mesh` 基类，实现了所有虚函数，并添加了三维规则网格特有的功能。
+
+## 主要API
+
+### 初始化
+```cpp
+// 初始化三维规则网格
+void init(std::string name,          // 网格名称
+          std::string info,          // 网格信息
+          int nx, int ny, int nz,    // x、y和z方向的网格数
+          double x0, double y0, double z0,  // 起始点坐标
+          double x1, double y1, double z1,  // 终止点坐标
+          double dx = 0.0,           // x方向网格间距（0表示自动计算）
+          double dy = 0.0,           // y方向网格间距（0表示自动计算）
+          double dz = 0.0);          // z方向网格间距（0表示自动计算）
+```
+
+### 网格信息
+```cpp
+// 获取网格点数
+void get_mesh_num(int &nx, int &ny, int &nz) const;
+
+// 获取网格范围
+void get_mesh_range(double &x0, double &y0, double &z0,    // 起点坐标
+                   double &x1, double &y1, double &z1) const; // 终点坐标
+
+// 获取网格间距
+void get_mesh_step(double &dx, double &dy, double &dz) const;
+```
+
+### 网格操作
+```cpp
+// 获取指定位置的网格点坐标
+void get_node_coord(int i, int j, int k,          // 节点索引
+                   double &x, double &y, double &z) const;
+
+// 获取指定位置的单元中心坐标
+void get_elem_center(int i, int j, int k,         // 单元索引
+                    double &x, double &y, double &z) const;
+
+// 获取指定坐标所在的网格单元索引
+bool get_elem_index(double x, double y, double z,  // 坐标
+                   int &i, int &j, int &k) const;  // 返回单元索引
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) override;
+
+// 保存为二进制格式
+virtual void save_binary(std::string filename) override;
+
+// 读取二进制格式
+virtual void load_binary(std::string filename) override;
+```
+
+## 使用示例
+
+```cpp
+// 创建三维规则网格
+gctl::regular_mesh_3d mesh3d;
+mesh3d.init("mesh3d", "3D regular mesh",
+            50, 50, 50,           // 50x50x50网格
+            0.0, 0.0, 0.0,        // 起点(0,0,0)
+            1.0, 1.0, 1.0);       // 终点(1,1,1)
+
+// 获取网格信息
+int nx, ny, nz;
+mesh3d.get_mesh_num(nx, ny, nz);
+
+double dx, dy, dz;
+mesh3d.get_mesh_step(dx, dy, dz);
+
+// 获取特定位置的坐标
+double x, y, z;
+mesh3d.get_node_coord(25, 25, 25, x, y, z);  // 获取中心点坐标
+
+// 查找特定坐标所在的单元
+int i, j, k;
+if (mesh3d.get_elem_index(0.5, 0.5, 0.5, i, j, k)) {
+    // 找到包含点(0.5,0.5,0.5)的单元
+}
+
+// 添加数据并保存
+mesh3d.add_data(gctl::NodeData, gctl::Scalar, "pressure", 0.0);
+mesh3d.save_gmsh("volume.msh", gctl::Packed);
+```
+
+## 注意事项
+
+1. 网格间距可以自动计算或手动指定
+2. 提供了三维空间中坐标与网格索引之间的转换功能
+3. 所有索引都是从0开始计数
+4. 网格点和单元的编号是规则的，遵循x-y-z优先顺序
+5. 支持Gmsh和二进制格式的文件操作
+6. 内存占用随网格规模的立方增长，使用大规模网格时需注意内存管理 

doc/regular_mesh_sph_3d.md

@@ -0,0 +1,115 @@
+# GCTL Regular Mesh Spherical 3D 文档
+
+## 简介
+
+`regular_mesh_sph_3d.h` 定义了GCTL网格库中的三维球面规则网格类 `regular_mesh_sph_3d`，用于创建和管理球面坐标系下的三维规则网格。该类特别适用于地球物理学中的球面网格计算。
+
+## 类继承
+
+`regular_mesh_sph_3d` 继承自 `base_mesh` 基类，实现了所有虚函数，并添加了球面网格特有的功能。
+
+## 主要API
+
+### 初始化
+```cpp
+// 初始化三维球面规则网格
+void init(std::string name,          // 网格名称
+          std::string info,          // 网格信息
+          double lat0,               // 起始纬度
+          double lon0,               // 起始经度
+          double r0,                 // 起始半径
+          double dlat,               // 纬度间隔
+          double dlon,               // 经度间隔
+          double dr,                 // 半径间隔
+          int nlat,                  // 纬度方向网格数
+          int nlon,                  // 经度方向网格数
+          int nr);                   // 径向网格数
+```
+
+### 网格信息
+```cpp
+// 获取网格点数
+void get_mesh_num(int &nlat, int &nlon, int &nr) const;
+
+// 获取网格范围
+void get_mesh_range(double &lat0, double &lon0, double &r0,    // 起点
+                   double &lat1, double &lon1, double &r1) const; // 终点
+
+// 获取网格间距
+void get_mesh_step(double &dlat, double &dlon, double &dr) const;
+```
+
+### 网格操作
+```cpp
+// 获取指定位置的网格点球面坐标
+void get_node_coord_sph(int i, int j, int k,              // 节点索引
+                       double &lat, double &lon, double &r) const;
+
+// 获取指定位置的网格点笛卡尔坐标
+void get_node_coord(int i, int j, int k,                  // 节点索引
+                   double &x, double &y, double &z) const;
+
+// 获取指定位置的单元中心球面坐标
+void get_elem_center_sph(int i, int j, int k,            // 单元索引
+                        double &lat, double &lon, double &r) const;
+
+// 获取指定位置的单元中心笛卡尔坐标
+void get_elem_center(int i, int j, int k,                // 单元索引
+                    double &x, double &y, double &z) const;
+
+// 获取指定球面坐标所在的网格单元索引
+bool get_elem_index(double lat, double lon, double r,    // 球面坐标
+                   int &i, int &j, int &k) const;        // 返回单元索引
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) override;
+
+// 保存为二进制格式
+virtual void save_binary(std::string filename) override;
+
+// 读取二进制格式
+virtual void load_binary(std::string filename) override;
+```
+
+## 使用示例
+
+```cpp
+// 创建三维球面规则网格
+gctl::regular_mesh_sph_3d mesh_sph;
+mesh_sph.init("earth_mesh", "Earth model mesh",
+              -90.0, -180.0, 6371.0,     // 起始点：南极点，经度-180度，地球半径
+              1.0, 1.0, 10.0,            // 间隔：1度纬度，1度经度，10km深度
+              180, 360, 100);            // 网格数：180纬度，360经度，100层
+
+// 获取网格信息
+int nlat, nlon, nr;
+mesh_sph.get_mesh_num(nlat, nlon, nr);
+
+// 获取特定位置的坐标（球面坐标）
+double lat, lon, r;
+mesh_sph.get_node_coord_sph(90, 180, 0, lat, lon, r);  // 获取中心点坐标
+
+// 获取特定位置的坐标（笛卡尔坐标）
+double x, y, z;
+mesh_sph.get_node_coord(90, 180, 0, x, y, z);
+
+// 添加数据并保存
+mesh_sph.add_data(gctl::NodeData, gctl::Scalar, "temperature", 0.0);
+mesh_sph.save_gmsh("earth_model.msh", gctl::Packed);
+```
+
+## 注意事项
+
+1. 坐标系采用地理坐标系约定：
+   - 纬度范围：-90度（南极）到90度（北极）
+   - 经度范围：-180度到180度
+   - 半径：从地心开始计算
+2. 提供了球面坐标和笛卡尔坐标的双重表示
+3. 网格在极点附近可能出现变形，需要特别处理
+4. 经度方向是周期性的，可以处理跨越180度经线的情况
+5. 支持Gmsh和二进制格式的文件操作
+6. 适用于地球物理建模、大气科学等领域 

doc/tet_mesh.md

@@ -0,0 +1,152 @@
+# GCTL Tetrahedral Mesh 文档
+
+## 简介
+
+`tet_mesh.h` 定义了GCTL网格库中的四面体网格类 `tet_mesh`，用于创建和管理三维空间中的非结构化四面体网格。该类特别适用于复杂三维区域的离散化，支持任意形状的体积域。
+
+## 类继承
+
+`tet_mesh` 继承自 `base_mesh` 基类，实现了所有虚函数，并添加了四面体网格特有的功能。
+
+## 主要API
+
+### 初始化
+```cpp
+// 初始化四面体网格
+void init(std::string name,                // 网格名称
+          std::string info,                // 网格信息
+          const array<point3dc> &nodes,    // 节点坐标数组
+          const array<int4> &elements);    // 单元连接关系数组
+```
+
+### 网格信息
+```cpp
+// 获取节点坐标
+void get_node_coord(int node_idx,          // 节点索引
+                   double &x, double &y, double &z) const;
+
+// 获取单元节点索引
+void get_elem_nodes(int elem_idx,          // 单元索引
+                   int &n1, int &n2, int &n3, int &n4) const;
+
+// 获取单元中心坐标
+void get_elem_center(int elem_idx,         // 单元索引
+                    double &x, double &y, double &z) const;
+
+// 获取单元体积
+double get_elem_volume(int elem_idx) const;
+
+// 获取单元外接球信息
+void get_elem_circumsphere(int elem_idx,   // 单元索引
+                          double &x, double &y, double &z,  // 球心坐标
+                          double &radius) const;           // 半径
+```
+
+### 网格查询
+```cpp
+// 查找包含指定点的单元
+bool find_element(double x, double y, double z,  // 目标点坐标
+                 int &elem_idx) const;           // 返回单元索引
+
+// 获取节点相邻的单元
+array<int> get_node_elements(int node_idx) const;
+
+// 获取单元相邻的单元
+array<int> get_elem_neighbors(int elem_idx) const;
+
+// 获取边界面
+array<int3> get_boundary_faces() const;
+```
+
+### 网格质量评估
+```cpp
+// 获取单元质量参数
+double get_elem_quality(int elem_idx) const;
+
+// 获取最小二面角
+double get_elem_min_dihedral_angle(int elem_idx) const;
+
+// 获取最大二面角
+double get_elem_max_dihedral_angle(int elem_idx) const;
+
+// 获取纵横比
+double get_elem_aspect_ratio(int elem_idx) const;
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) override;
+
+// 保存为二进制格式
+virtual void save_binary(std::string filename) override;
+
+// 读取二进制格式
+virtual void load_binary(std::string filename) override;
+```
+
+## 使用示例
+
+```cpp
+// 创建节点坐标数组
+array<point3dc> nodes;
+nodes.push_back(point3dc(0.0, 0.0, 0.0));  // 底面三角形的顶点
+nodes.push_back(point3dc(1.0, 0.0, 0.0));
+nodes.push_back(point3dc(0.0, 1.0, 0.0));
+nodes.push_back(point3dc(0.0, 0.0, 1.0));  // 顶点
+
+// 创建单元连接关系
+array<int4> elements;
+elements.push_back(int4(0, 1, 2, 3));  // 一个四面体
+
+// 创建四面体网格
+gctl::tet_mesh mesh;
+mesh.init("volume", "3D tetrahedral mesh", nodes, elements);
+
+// 获取单元信息
+double volume = mesh.get_elem_volume(0);  // 获取四面体的体积
+double quality = mesh.get_elem_quality(0); // 获取四面体的质量参数
+
+// 获取外接球信息
+double cx, cy, cz, radius;
+mesh.get_elem_circumsphere(0, cx, cy, cz, radius);
+
+// 查找包含特定点的单元
+int elem_idx;
+if (mesh.find_element(0.25, 0.25, 0.25, elem_idx)) {
+    // 找到包含点(0.25,0.25,0.25)的四面体
+}
+
+// 获取边界
+array<int3> boundary = mesh.get_boundary_faces();
+
+// 添加数据并保存
+mesh.add_data(gctl::NodeData, gctl::Vector, "displacement",
+              array<point3dc>(nodes.size(), point3dc(0,0,0)));
+mesh.save_gmsh("volume.msh", gctl::Packed);
+```
+
+## 注意事项
+
+1. 四面体网格适用于三维区域的离散化
+2. 每个单元由四个节点定义，遵循右手法则定义方向
+3. 节点和单元的编号从0开始
+4. 提供了完整的网格质量评估功能：
+   - 体积
+   - 二面角
+   - 纵横比
+   - 质量参数
+5. 支持复杂的拓扑关系查询
+6. 可以自动识别和提取边界面
+7. 支持Gmsh和二进制格式的文件操作
+8. 适用于：
+   - 有限元分析
+   - 计算流体动力学
+   - 结构力学
+   - 电磁场分析
+9. 网格质量对计算精度和收敛性有重要影响，建议：
+   - 控制最小二面角不小于10度
+   - 控制最大二面角不大于160度
+   - 保持适当的纵横比
+</rewritten_file> 

doc/tri_mesh.md

@@ -0,0 +1,138 @@
+# GCTL Triangle Mesh 文档
+
+## 简介
+
+`tri_mesh.h` 定义了GCTL网格库中的三角形网格类 `tri_mesh`，用于创建和管理二维平面上的非结构化三角形网格。该类特别适用于复杂二维区域的离散化，支持任意形状的边界。
+
+## 类继承
+
+`tri_mesh` 继承自 `base_mesh` 基类，实现了所有虚函数，并添加了三角形网格特有的功能。
+
+## 主要API
+
+### 初始化
+```cpp
+// 初始化三角形网格
+void init(std::string name,                // 网格名称
+          std::string info,                // 网格信息
+          const array<point2dc> &nodes,    // 节点坐标数组
+          const array<int3> &elements);    // 单元连接关系数组
+```
+
+### 网格信息
+```cpp
+// 获取节点坐标
+void get_node_coord(int node_idx,          // 节点索引
+                   double &x, double &y) const;
+
+// 获取单元节点索引
+void get_elem_nodes(int elem_idx,          // 单元索引
+                   int &n1, int &n2, int &n3) const;
+
+// 获取单元中心坐标
+void get_elem_center(int elem_idx,         // 单元索引
+                    double &x, double &y) const;
+
+// 获取单元面积
+double get_elem_area(int elem_idx) const;
+
+// 获取单元外接圆信息
+void get_elem_circumcircle(int elem_idx,   // 单元索引
+                          double &x, double &y,  // 圆心坐标
+                          double &radius) const; // 半径
+```
+
+### 网格查询
+```cpp
+// 查找包含指定点的单元
+bool find_element(double x, double y,      // 目标点坐标
+                 int &elem_idx) const;     // 返回单元索引
+
+// 获取节点相邻的单元
+array<int> get_node_elements(int node_idx) const;
+
+// 获取单元相邻的单元
+array<int> get_elem_neighbors(int elem_idx) const;
+
+// 获取边界边
+array<int2> get_boundary_edges() const;
+```
+
+### 网格质量评估
+```cpp
+// 获取单元质量参数
+double get_elem_quality(int elem_idx) const;
+
+// 获取最小内角
+double get_elem_min_angle(int elem_idx) const;
+
+// 获取最大内角
+double get_elem_max_angle(int elem_idx) const;
+```
+
+### 文件操作
+```cpp
+// 保存为Gmsh格式
+virtual void save_gmsh(std::string filename,
+                      index_packed_e packed = Packed) override;
+
+// 保存为二进制格式
+virtual void save_binary(std::string filename) override;
+
+// 读取二进制格式
+virtual void load_binary(std::string filename) override;
+```
+
+## 使用示例
+
+```cpp
+// 创建节点坐标数组
+array<point2dc> nodes;
+nodes.push_back(point2dc(0.0, 0.0));
+nodes.push_back(point2dc(1.0, 0.0));
+nodes.push_back(point2dc(0.5, 0.866));
+nodes.push_back(point2dc(0.5, 0.289));
+
+// 创建单元连接关系
+array<int3> elements;
+elements.push_back(int3(0, 1, 3));  // 第一个三角形
+elements.push_back(int3(0, 3, 2));  // 第二个三角形
+elements.push_back(int3(1, 2, 3));  // 第三个三角形
+
+// 创建三角形网格
+gctl::tri_mesh mesh;
+mesh.init("domain", "2D triangle mesh", nodes, elements);
+
+// 获取单元信息
+double area = mesh.get_elem_area(0);  // 获取第一个三角形的面积
+double quality = mesh.get_elem_quality(0);  // 获取第一个三角形的质量参数
+
+// 获取外接圆信息
+double cx, cy, radius;
+mesh.get_elem_circumcircle(0, cx, cy, radius);
+
+// 查找包含特定点的单元
+int elem_idx;
+if (mesh.find_element(0.5, 0.5, elem_idx)) {
+    // 找到包含点(0.5,0.5)的三角形
+}
+
+// 获取边界
+array<int2> boundary = mesh.get_boundary_edges();
+
+// 添加数据并保存
+mesh.add_data(gctl::NodeData, gctl::Scalar, "pressure", 0.0);
+mesh.save_gmsh("domain.msh", gctl::Packed);
+```
+
+## 注意事项
+
+1. 三角形网格适用于二维区域的离散化
+2. 每个单元由三个节点定义，按逆时针顺序排列
+3. 节点和单元的编号从0开始
+4. 提供了完整的网格质量评估功能
+5. 支持复杂的拓扑关系查询
+6. 可以自动识别和提取边界边
+7. 支持Gmsh和二进制格式的文件操作
+8. 适用于有限元分析、流体动力学等数值模拟
+9. 网格质量对计算结果有重要影响，建议进行质量检查 

example/CMakeLists.txt

@@ -21,7 +21,7 @@ add_example(mesh_ex1 ON)
 add_example(mesh_ex2 ON)
 add_example(mesh_ex3 ON)
 add_example(mesh_ex4 ON)
-add_example(mesh_ex5 ON)
+add_example(mesh_ex5 OFF)
 add_example(mesh_ex6 ON)
 add_example(mesh_ex7 ON)
 add_example(mesh_ex8 ON)

example/mesh_ex1.cpp

@@ -32,26 +32,20 @@ int main(int argc, char *argv[])
 	gctl::regular_grid rgd;
 	rgd.init("grid-1", "null", 4, 4, 0.0, 0.0, 1.0, 1.0);
 
-	gctl::meshdata *data = rgd.add_data("data-1", gctl::NodeData, true, 2.5);
+	rgd.add_data(gctl::NodeData, gctl::Scalar, "data-1", 2.5);
-	gctl::meshdata *data2= rgd.add_data("data-2", gctl::ElemData, false, 1.2);
+	rgd.add_data(gctl::ElemData, gctl::Scalar, "data-2", 1.2, false);
-	gctl::meshdata *data3= rgd.add_data("data-3", gctl::NodeData, true, gctl::Scalar);
+	rgd.add_data(gctl::NodeData, gctl::Scalar, "data-3", 1.0);
+	rgd.show_info();
 
-	std::cout << "data name: " << data2->get_datname() << std::endl;
+	gctl::meshdata data2 = rgd.get_data("data-2");
-	//gctl::array<double> *data_ptr = (gctl::array<double>*) rgd.get_datval(data2->get_name());
+	std::cout << "data name: " << data2.name_ << std::endl;
-	gctl::array<double> *data_ptr = (gctl::array<double>*) data2->get_datval_ptr();
+	data2.datval_.show();
-	for (int i = 0; i < data_ptr->size(); i++)
-	{
-		std::cout << data_ptr->at(i) << std::endl;
-	}
 
-	std::cout << "data name: " << data3->get_datname() << std::endl;
+	gctl::meshdata data3 = rgd.get_data("data-3");
-	gctl::array<int> *data_ptr2 = (gctl::array<int>*) rgd.get_datval(data3->get_datname());
+	std::cout << "data name: " << data3.name_ << std::endl;
-	for (int i = 0; i < data_ptr2->size(); i++)
+	data3.datval_.show();
-	{
-		std::cout << data_ptr2->at(i) << std::endl;
-	}
 
 	rgd.remove_data("data-1");
-
+	rgd.show_info();
 	return 0;
 }

example/mesh_ex10.cpp

@@ -31,7 +31,7 @@ int main(int argc, char *argv[])
 {
 	gctl::regular_mesh_sph_3d rm_3ds;
 	rm_3ds.init("mesh-1", "null", 30.25, 30.25, 2005, 0.5, 0.5, 10, 40, 40, 50);
-	rm_3ds.add_data("data-1", gctl::ElemData, true, 2.5);
+	rm_3ds.add_data(gctl::ElemData, gctl::Scalar, "data-1", 2.5);
-	rm_3ds.save_gmsh("data/out/mesh_sample10", gctl::ElemData, gctl::OverWrite, gctl::NotPacked);
+	rm_3ds.save_gmsh_withdata("mesh_sample10",gctl::OverWrite, gctl::NotPacked);
 	return 0;
 }

example/mesh_ex2.cpp

@@ -34,32 +34,29 @@ int main(int argc, char *argv[])
 
 	gctl::regular_grid rgd;
 	rgd.init("grid-1", "null", 15, 10, 0.0, 0.0, 1.0, 1.0);
-	rgd.add_data(dname, gctl::NodeData, true, gctl::Scalar);
 
-	gctl::array<double> *data_ptr = (gctl::array<double>*) rgd.get_datval(dname);
+	gctl::meshdata &data = rgd.add_data(gctl::NodeData, gctl::Scalar, dname, 0.0);
 	for (int j = 0; j < rgd.get_ydim(); j++)
 	{
 		for (int i = 0; i < rgd.get_xdim(); i++)
 		{
-			data_ptr->at(i + j*rgd.get_xdim()) = i;
+			data.datval_[i + j*rgd.get_xdim()] = i;
 		}
 	}
 
-	rgd.add_data(dname2, gctl::NodeData, true, gctl::Scalar);
+	gctl::meshdata &data2 = rgd.add_data(gctl::NodeData, gctl::Scalar, dname2, 0.0);
-	data_ptr = (gctl::array<double>*) rgd.get_datval(dname2);
 	for (int j = 0; j < rgd.get_ydim(); j++)
 	{
 		for (int i = 0; i < rgd.get_xdim(); i++)
 		{
-			data_ptr->at(i + j*rgd.get_xdim()) = j;
+			data2.datval_[i + j*rgd.get_xdim()] = j;
 		}
 	}
 
 	// disable the output of data object named dname2
-	gctl::meshdata *data2 = rgd.get_data(dname2);
+	//data2.output_ok_ = false;
-	data2->set_output(false);
-
-	rgd.save_netcdf_grid("data/out/sample2-out", gctl::NodeData);
 
+	rgd.show_info();
+	rgd.save_netcdf_grid("sample2-out", gctl::NodeData);
 	return 0;
 }

example/mesh_ex3.cpp

@@ -27,62 +27,56 @@
 
 #include "../lib/mesh.h"
 
-int main(int argc, char *argv[])
+int main(int argc, char *argv[]) try
 {
-	try
+	std::string dname = "example";
+	std::string dname2= "another-example";
+
+	gctl::regular_grid rgd;
+	rgd.init("grid-1", "this is a test message.", 15, 10, 0.0, 0.0, 1.0, 1.0);
+
+	gctl::meshdata &data = rgd.add_data(gctl::NodeData, gctl::Scalar, dname, 0.0);
+	for (int j = 0; j < rgd.get_ydim(); j++)
 	{
-		std::string dname = "example";
+		for (int i = 0; i < rgd.get_xdim(); i++)
-		std::string dname2= "another-example";
-
-		gctl::regular_grid rgd;
-		rgd.init("grid-1", "this is a test message.", 15, 10, 0.0, 0.0, 1.0, 1.0);
-		rgd.add_data(dname, gctl::NodeData, true, gctl::Scalar);
-
-		gctl::array<double> *data_ptr = (gctl::array<double>*) rgd.get_datval(dname);
-		for (int j = 0; j < rgd.get_ydim(); j++)
 		{
-			for (int i = 0; i < rgd.get_xdim(); i++)
+			data.datval_[i + j*rgd.get_xdim()] = i;
-			{
-				data_ptr->at(i + j*rgd.get_xdim()) = i;
-			}
 		}
-
-		rgd.add_data(dname2, gctl::ElemData, true, gctl::Scalar);
-		data_ptr = (gctl::array<double>*) rgd.get_datval(dname2);
-		for (int j = 0; j < rgd.get_ydim()-1; j++)
-		{
-			for (int i = 0; i < rgd.get_xdim()-1; i++)
-			{
-				data_ptr->at(i + j*(rgd.get_xdim()-1)) = j;
-			}
-		}
-
-		rgd.save_binary("data/out/sample3-out");
-		rgd.save_netcdf_grid("data/out/sample3-out1", dname);
-		rgd.save_netcdf_grid("data/out/sample3-out2", dname2);
-
-		gctl::regular_grid rgd2;
-		rgd2.load_binary("data/out/sample3-out");
-
-		rgd2.show_info();
-
-		data_ptr = (gctl::array<double>*) rgd2.get_datval(dname);
-		for (int i = 0; i < data_ptr->size(); i++)
-		{
-			std::cout << data_ptr->at(i) << " ";
-		}
-		std::cout << std::endl;
-
-		data_ptr = (gctl::array<double>*) rgd2.get_datval(dname2);
-		for (int i = 0; i < data_ptr->size(); i++)
-		{
-			std::cout << data_ptr->at(i) << " ";
-		}
-		std::cout << std::endl;
 	}
-	catch(std::exception &e)
+
+	gctl::meshdata &data2 = rgd.add_data(gctl::ElemData, gctl::Scalar, dname2, 0.0);
+	for (int j = 0; j < rgd.get_ydim()-1; j++)
 	{
-		GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+		for (int i = 0; i < rgd.get_xdim()-1; i++)
+		{
+			data2.datval_[i + j*(rgd.get_xdim()-1)] = j;
+		}
 	}
+
+	rgd.save_binary("sample3-out");
+	rgd.save_netcdf_grid("sample3-out1", dname);
+	rgd.save_netcdf_grid("sample3-out2", dname2);
+
+	gctl::regular_grid rgd2;
+	rgd2.load_binary("sample3-out");
+	rgd2.show_info();
+
+	gctl::meshdata &data3 = rgd2.get_data(dname);
+	for (int i = 0; i < data3.datval_.size(); i++)
+	{
+		std::cout << data3.datval_[i] << " ";
+	}
+	std::cout << std::endl;
+
+	gctl::meshdata &data4 = rgd2.get_data(dname2);
+	for (int i = 0; i < data4.datval_.size(); i++)
+	{
+		std::cout << data4.datval_[i] << " ";
+	}
+	std::cout << std::endl;
 	return 0;
 }
+catch(std::exception &e)
+{
+	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}

example/mesh_ex4.cpp

@@ -27,37 +27,34 @@
 
 #include "../lib/mesh.h"
 
-int main(int argc, char *argv[])
+int main(int argc, char *argv[]) try
 {
-	try
+	std::string dname = "vector-data";
+
+	gctl::regular_grid rgd;
+	rgd.init("grid-example", "this grid can store vectors. Yeah!", 15, 10, 0.0, 0.0, 1.0, 1.0);
+
+	gctl::array<gctl::point3dc> init_data(rgd.get_ydim()*rgd.get_xdim());
+	for (int j = 0; j < rgd.get_ydim(); j++)
 	{
-		std::string dname = "vector-data";
+		for (int i = 0; i < rgd.get_xdim(); i++)
-
-		gctl::regular_grid rgd;
-		rgd.init("grid-example", "this grid can store vectors. Yeah!", 15, 10, 0.0, 0.0, 1.0, 1.0);
-		rgd.add_data(dname, gctl::NodeData, true, gctl::Vector);
-
-		gctl::array<gctl::point3dc> *data_ptr = (gctl::array<gctl::point3dc>*) rgd.get_datval(dname);
-		for (int j = 0; j < rgd.get_ydim(); j++)
 		{
-			for (int i = 0; i < rgd.get_xdim(); i++)
+			init_data[i + j*rgd.get_xdim()].x = i;
-			{
+			init_data[i + j*rgd.get_xdim()].y = j;
-				data_ptr->at(i + j*rgd.get_xdim()).x = i;
+			init_data[i + j*rgd.get_xdim()].z = i+j;
-				data_ptr->at(i + j*rgd.get_xdim()).y = j;
-				data_ptr->at(i + j*rgd.get_xdim()).z = i+j;
-			}
 		}
-
-		rgd.save_binary("data/out/sample4-out");
-
-		gctl::regular_grid rgd2;
-		rgd2.load_binary("data/out/sample4-out");
-
-		rgd2.show_info();
-	}
-	catch(std::exception &e)
-	{
-		GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
 	}
+
+	gctl::meshdata data = rgd.add_data(gctl::NodeData, dname, init_data);
+
+	rgd.save_binary("sample4-out");
+
+	gctl::regular_grid rgd2;
+	rgd2.load_binary("sample4-out");
+	rgd2.show_info();
 	return 0;
 }
+catch(std::exception &e)
+{
+	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}

example/mesh_ex5.cpp

@@ -27,20 +27,17 @@
 
 #include "../lib/mesh.h"
 
-int main(int argc, char *argv[])
+int main(int argc, char *argv[]) try
 {
-	try
+	gctl::regular_grid rgd;
-	{
+	rgd.load_netcdf_grid("data/out/sample3-out1", gctl::ElemData, "x", "y");
-		gctl::regular_grid rgd;
-		rgd.load_netcdf_grid("data/out/sample3-out1", gctl::ElemData, "x", "y");
 
-		rgd.show_info();
+	rgd.show_info();
 
-		rgd.save_netcdf_grid("data/out/sample5-out", "example");
+	rgd.save_netcdf_grid("data/out/sample5-out", "example");
-	}
-	catch(std::exception &e)
-	{
-		GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
-	}
 	return 0;
 }
+catch(std::exception &e)
+{
+	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}

example/mesh_ex6.cpp

@@ -27,46 +27,40 @@
 
 #include "../lib/mesh.h"
 
-int main(int argc, char *argv[])
+int main(int argc, char *argv[]) try
 {
-	try
+	gctl::regular_mesh_2d rg_mesh;
+	rg_mesh.init("mesh-example", "this mesh can store vectors. Yeah!", 15, 10, 0.5, 0.5, 1.0, 1.0);
+
+	int ybnum = rg_mesh.get_ybnum();
+	int xbnum = rg_mesh.get_xbnum();
+
+	gctl::array<gctl::point3dc> init_data(ybnum*xbnum);
+	for (int j = 0; j < ybnum; j++)
 	{
-		gctl::regular_mesh_2d rg_mesh;
+		for (int i = 0; i < xbnum; i++)
-		rg_mesh.init("mesh-example", "this mesh can store vectors. Yeah!", 15, 10, 0.5, 0.5, 1.0, 1.0);
-
-		int ybnum = rg_mesh.get_ybnum();
-		int xbnum = rg_mesh.get_xbnum();
-
-		rg_mesh.add_data("vector-data", gctl::ElemData, true, gctl::Vector);
-		gctl::array<gctl::point3dc> *data_ptr = (gctl::array<gctl::point3dc>*) rg_mesh.get_datval("vector-data");
-		for (int j = 0; j < ybnum; j++)
 		{
-			for (int i = 0; i < xbnum; i++)
+			init_data[i + j*xbnum].x = i;
-			{
+			init_data[i + j*xbnum].y = j;
-				data_ptr->at(i + j*xbnum).x = i;
+			init_data[i + j*xbnum].z = i+j;
-				data_ptr->at(i + j*xbnum).y = j;
-				data_ptr->at(i + j*xbnum).z = i+j;
-			}
 		}
-
-		rg_mesh.add_data("double-data", gctl::NodeData, true, gctl::Scalar);
-		gctl::array<double> *data_ptr2 = (gctl::array<double>*) rg_mesh.get_datval("double-data");
-		for (int j = 0; j < ybnum+1; j++)
-		{
-			for (int i = 0; i < xbnum+1; i++)
-			{
-				data_ptr2->at(i + j*(xbnum+1)) = i+j;
-			}
-		}
-
-		//rg_mesh.save_gmsh("data/sample6-out", "vector-data");
-		//rg_mesh.save_gmsh("data/sample6-out", "double-data", Append);
-
-		rg_mesh.save_gmsh("data/out/sample6-out", gctl::ElemData, gctl::OverWrite);
 	}
-	catch(std::exception &e)
+
+	rg_mesh.add_data(gctl::ElemData, "vector-data", init_data);
+	
+	gctl::meshdata &data2= rg_mesh.add_data(gctl::NodeData, gctl::Scalar, "double-data", 0.0);
+	for (int j = 0; j < ybnum+1; j++)
 	{
-		GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+		for (int i = 0; i < xbnum+1; i++)
+		{
+			data2.datval_[i + j*(xbnum+1)] = i+j;
+		}
 	}
+
+	rg_mesh.save_gmsh_withdata("sample6-out", gctl::OverWrite);
 	return 0;
 }
+catch(std::exception &e)
+{
+	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}

example/mesh_ex7.cpp

@@ -27,35 +27,31 @@
 
 #include "../lib/mesh.h"
 
-int main(int argc, char *argv[])
+int main(int argc, char *argv[]) try
 {
-	try
+	gctl::regular_mesh_3d rg_mesh;
+	rg_mesh.init("mesh-example", "this mesh can store vectors. Yeah!", 15, 10, 10, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0);
+
+	int ybnum = rg_mesh.get_ybnum();
+	int xbnum = rg_mesh.get_xbnum();
+	int zbnum = rg_mesh.get_zbnum();
+
+	gctl::meshdata &data = rg_mesh.add_data(gctl::NodeData, gctl::Scalar, "double-data", 0.0);
+	for (int k = 0; k < zbnum+1; k++)
 	{
-		gctl::regular_mesh_3d rg_mesh;
+		for (int j = 0; j < ybnum+1; j++)
-		rg_mesh.init("mesh-example", "this mesh can store vectors. Yeah!", 15, 10, 10, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0);
-
-		int ybnum = rg_mesh.get_ybnum();
-		int xbnum = rg_mesh.get_xbnum();
-		int zbnum = rg_mesh.get_zbnum();
-
-		rg_mesh.add_data("double-data", gctl::NodeData, gctl::Packed, gctl::Scalar);
-		gctl::array<double> *data_ptr = (gctl::array<double>*) rg_mesh.get_datval("double-data");
-		for (int k = 0; k < zbnum+1; k++)
 		{
-			for (int j = 0; j < ybnum+1; j++)
+			for (int i = 0; i < xbnum+1; i++)
 			{
-				for (int i = 0; i < xbnum+1; i++)
+				data.datval_[i + j*(xbnum+1) + k*(xbnum+1)*(ybnum+1)] = i+j+k;
-				{
-					data_ptr->at(i + j*(xbnum+1) + k*(xbnum+1)*(ybnum+1)) = i+j+k;
-				}
 			}
 		}
+	}
 
-		rg_mesh.save_gmsh("data/out/sample7-out", "double-data", gctl::OverWrite);
+	rg_mesh.save_gmsh_withdata("sample7-out", "double-data", gctl::OverWrite);
-	}
-	catch(std::exception &e)
-	{
-		GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
-	}
 	return 0;
 }
+catch(std::exception &e)
+{
+	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}

example/mesh_ex8.cpp

@@ -27,24 +27,21 @@
 
 #include "../lib/mesh.h"
 
-int main(int argc, char *argv[])
+int main(int argc, char *argv[]) try
 {
-	try
+	gctl::triangle2d_mesh t_mesh;
-	{
+	t_mesh.load_triangle("sample8", gctl::Packed);
-		gctl::triangle_mesh t_mesh;
+	t_mesh.add_data(gctl::ElemData, gctl::Scalar, "example", 1.0);
-		t_mesh.load_triangle("data/out/sample8", gctl::Packed);
-		t_mesh.add_data("example", gctl::ElemData, gctl::Packed, 1.0);
 
-		t_mesh.save_gmsh("data/out/sample8-out", "example", gctl::OverWrite, gctl::NotPacked);
+	t_mesh.save_gmsh_withdata("sample8-out", "example", gctl::OverWrite, gctl::NotPacked);
-		t_mesh.save_binary("data/out/sample8-out");
+	t_mesh.save_binary("sample8-out");
 
-		gctl::triangle_mesh t2_mesh;
+	gctl::triangle2d_mesh t2_mesh;
-		t2_mesh.load_binary("data/out/sample8-out");
+	t2_mesh.load_binary("sample8-out");
-		t2_mesh.show_info();
+	t2_mesh.show_info();
-	}
-	catch(std::exception &e)
-	{
-		GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
-	}
 	return 0;
 }
+catch(std::exception &e)
+{
+	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}

example/mesh_ex9.cpp

@@ -27,43 +27,39 @@
 
 #include "../lib/mesh.h"
 
-int main(int argc, char *argv[])
+int main(int argc, char *argv[]) try
 {
-	try
+	gctl::array<double> xs(15);
+	for (int i = 0; i < 15; i++)
 	{
-		gctl::array<double> xs(15);
+		xs[i] = i+1;
-		for (int i = 0; i < 15; i++)
-		{
-			xs[i] = i+1;
-		}
-
-		gctl::array<double> ys(10);
-		for (int i = 0; i < 10; i++)
-		{
-			ys[i] = i+1;
-		}
-
-		gctl::linear_mesh_2d l2d_mesh;
-		l2d_mesh.init("mesh-example", "This is a linear mesh", 0.5, 0.5, xs, ys);
-
-		int ybnum = l2d_mesh.get_ybnum();
-		int xbnum = l2d_mesh.get_xbnum();
-
-		l2d_mesh.add_data("double-data", gctl::ElemData, gctl::Packed, gctl::Scalar);
-		gctl::array<double> *data_ptr2 = (gctl::array<double>*) l2d_mesh.get_datval("double-data");
-		for (int j = 0; j < ybnum; j++)
-		{
-			for (int i = 0; i < xbnum; i++)
-			{
-				data_ptr2->at(i + j*xbnum) = i+j;
-			}
-		}
-
-		l2d_mesh.save_gmsh("data/out/sample9-out", gctl::ElemData, gctl::OverWrite, gctl::NotPacked);
 	}
-	catch(std::exception &e)
+
+	gctl::array<double> ys(10);
+	for (int i = 0; i < 10; i++)
 	{
-		GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+		ys[i] = i+1;
 	}
+
+	gctl::linear_mesh_2d l2d_mesh;
+	l2d_mesh.init("mesh-example", "This is a linear mesh", 0.5, 0.5, xs, ys);
+
+	int ybnum = l2d_mesh.get_ybnum();
+	int xbnum = l2d_mesh.get_xbnum();
+
+	gctl::meshdata &data = l2d_mesh.add_data(gctl::ElemData, gctl::Scalar, "double-data", 0.0);
+	for (int j = 0; j < ybnum; j++)
+	{
+		for (int i = 0; i < xbnum; i++)
+		{
+			data.datval_[i + j*xbnum] = i+j;
+		}
+	}
+
+	l2d_mesh.save_gmsh_withdata("sample9-out", gctl::OverWrite, gctl::NotPacked);
 	return 0;
 }
+catch(std::exception &e)
+{
+	GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}

lib/mesh.h

@@ -35,7 +35,7 @@
 #include "mesh/regular_mesh_3d.h"
 #include "mesh/linear_mesh_2d.h"
 #include "mesh/linear_mesh_3d.h"
-#include "mesh/tri_mesh.h"
+#include "mesh/tri2d_mesh.h"
 #include "mesh/tet_mesh.h"
 #include "mesh/regular_mesh_sph_3d.h"
 

lib/mesh/linear_mesh_2d.cpp

@@ -30,11 +30,7 @@
 void gctl::linear_mesh_2d::init(std::string in_name, std::string in_info, double xmin, double ymin, 
 	const gctl::array<double> &xsizes, const gctl::array<double> &ysizes)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is already initialized.");
-	}
-
 	base_mesh::init(LINEAR_MESH, MESH_2D, in_name, in_info);
 
 	lm_xmin  = xmin;
@@ -44,8 +40,8 @@ void gctl::linear_mesh_2d::init(std::string in_name, std::string in_info, double
 	lm_xsizes.resize(lm_xbnum);
 	lm_ysizes.resize(lm_ybnum);
 
-	node_num = (lm_xbnum+1) * (lm_ybnum+1);
+	node_num_ = (lm_xbnum+1) * (lm_ybnum+1);
-	ele_num  = lm_xbnum * lm_ybnum;
+	ele_num_  = lm_xbnum * lm_ybnum;
 
 	for (int i = 0; i < lm_xbnum; i++)
 	{
@@ -57,8 +53,8 @@ void gctl::linear_mesh_2d::init(std::string in_name, std::string in_info, double
 		lm_ysizes[i] = ysizes[i];
 	}
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	double x_p, y_p = lm_ymin - 0.5*lm_ysizes[0];
 	for (int j = 0; j < lm_ybnum+1; j++)
@@ -88,16 +84,23 @@ void gctl::linear_mesh_2d::init(std::string in_name, std::string in_info, double
 		}
 	}
 
-	initialized = true;
+	initialized_ = true;
+	return;
+}
+
+void gctl::linear_mesh_2d::show_mesh_dimension(std::ostream &os) const
+{
+	os << "node num: " << node_num_ << std::endl;
+	os << "elem num: " << ele_num_ << std::endl;
+	os << "x-range: " << get_xmin() << " -> " << get_xmax() << "\n";
+	os << "y-range: " << get_ymin() << " -> " << get_ymax() << "\n";
+	os << "dimension: " << lm_xbnum << ", " << lm_ybnum << "\n";
 	return;
 }
 
 void gctl::linear_mesh_2d::load_binary(std::string filename)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is already initialized.");
-	}
 
 	std::ifstream infile;
 	gctl::open_infile(infile, filename, ".2m", std::ios::in|std::ios::binary);
@@ -113,8 +116,8 @@ void gctl::linear_mesh_2d::load_binary(std::string filename)
 
 	lm_xsizes.resize(lm_xbnum);
 	lm_ysizes.resize(lm_ybnum);
-	node_num = (lm_xbnum+1) * (lm_ybnum+1);
+	node_num_ = (lm_xbnum+1) * (lm_ybnum+1);
-	ele_num  = lm_xbnum * lm_ybnum;
+	ele_num_  = lm_xbnum * lm_ybnum;
 
 	for (int i = 0; i < lm_xbnum; i++)
 	{
@@ -126,8 +129,8 @@ void gctl::linear_mesh_2d::load_binary(std::string filename)
 		infile.read((char*)lm_ysizes.get(i), sizeof(double));
 	}
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	double x_p, y_p = lm_ymin - 0.5*lm_ysizes[0];
 	for (int j = 0; j < lm_ybnum+1; j++)
@@ -157,7 +160,7 @@ void gctl::linear_mesh_2d::load_binary(std::string filename)
 		}
 	}
 
-	initialized = true;
+	initialized_ = true;
 
 	// 读入模型数据单元
 	load_datablock(infile);
@@ -168,10 +171,7 @@ void gctl::linear_mesh_2d::load_binary(std::string filename)
 
 void gctl::linear_mesh_2d::save_binary(std::string filename)
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is not initialized.");
-	}
 
 	std::ofstream outfile;
 	gctl::open_outfile(outfile, filename, ".2m", std::ios::out|std::ios::binary);
@@ -214,61 +214,63 @@ gctl::linear_mesh_2d::~linear_mesh_2d(){}
 
 int gctl::linear_mesh_2d::get_xbnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is not initialized.");
-	}
-
 	return lm_xbnum;
 }
 
 int gctl::linear_mesh_2d::get_ybnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is not initialized.");
-	}
-
 	return lm_ybnum;
 }
 
 double gctl::linear_mesh_2d::get_xmin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is not initialized.");
-	}
-
 	return lm_xmin;
 }
 
 double gctl::linear_mesh_2d::get_ymin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is not initialized.");
-	}
-
 	return lm_ymin;
 }
 
+double gctl::linear_mesh_2d::get_xmax() const
+{
+	check_initiated(); // 检查网格是否已经初始化
+
+	double xmax = lm_xmin - 0.5*lm_xsizes[0];
+	for (size_t i = 0; i < lm_xsizes.size(); i++)
+	{
+		xmax += lm_xsizes[i];
+	}
+	xmax -= 0.5*lm_xsizes.back();
+	return xmax;
+}
+
+double gctl::linear_mesh_2d::get_ymax() const
+{
+	check_initiated(); // 检查网格是否已经初始化
+
+	double ymax = lm_ymin - 0.5*lm_ysizes[0];
+	for (size_t i = 0; i < lm_ysizes.size(); i++)
+	{
+		ymax += lm_ysizes[i];
+	}
+	ymax -= 0.5*lm_ysizes.back();
+	return ymax;
+}
+
 const gctl::array<double> *gctl::linear_mesh_2d::get_xsizes() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is not initialized.");
-	}
-
 	return &lm_xsizes;
 }
 
 const gctl::array<double> *gctl::linear_mesh_2d::get_ysizes() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_2d] The mesh is not initialized.");
-	}
-
 	return &lm_ysizes;
 }
 
@@ -280,16 +282,3 @@ void gctl::linear_mesh_2d::save_gmsh(std::string filename, index_packed_e packed
 	outfile.close();
 	return;
 }
-
-
-void gctl::linear_mesh_2d::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
-	return;
-}
-
-void gctl::linear_mesh_2d::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, datname, out_mode, packed);
-	return;
-}

lib/mesh/linear_mesh_2d.h

@@ -40,12 +40,14 @@ namespace gctl
 		 * 以下为mesh类型的虚函数实现
 		 */
 
-		void init(std::string in_name, std::string in_info, double xmin, double ymin, 
+		void load_gmsh(std::string filename, index_packed_e packed = Packed){}
-			const array<double> &xsizes, const array<double> &ysizes);
+		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;
+
 		/**
 		 * 以下为regular_mesh_2d的专有函数
 		 */
@@ -55,17 +57,18 @@ namespace gctl
 			const array<double> &xsizes, const array<double> &ysizes);
 		virtual ~linear_mesh_2d();
 
+		void init(std::string in_name, std::string in_info, double xmin, double ymin, 
+			const array<double> &xsizes, const array<double> &ysizes);
+
 		int get_xbnum() const;
 		int get_ybnum() const;
 		double get_xmin() const;
 		double get_ymin() const;
+		double get_xmax() const;
+		double get_ymax() const;
 		const array<double> *get_xsizes() const;
 		const array<double> *get_ysizes() const;
 
-		void save_gmsh(std::string filename, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
-
 	protected:
 		int lm_xbnum, lm_ybnum;
 		double lm_xmin, lm_ymin;

lib/mesh/linear_mesh_3d.cpp

@@ -30,11 +30,7 @@
 void gctl::linear_mesh_3d::init(std::string in_name, std::string in_info, double xmin, double ymin, 
 	double zmin, const array<double> &xsizes, const array<double> &ysizes, const array<double> &zsizes)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is already initialized.");
-	}
-
 	base_mesh::init(LINEAR_MESH, MESH_3D, in_name, in_info);
 
 	lm_xmin  = xmin;
@@ -47,8 +43,8 @@ void gctl::linear_mesh_3d::init(std::string in_name, std::string in_info, double
 	lm_ysizes.resize(lm_ybnum);
 	lm_zsizes.resize(lm_zbnum);
 
-	node_num = (lm_xbnum+1) * (lm_ybnum+1) * (lm_zbnum+1);
+	node_num_ = (lm_xbnum+1) * (lm_ybnum+1) * (lm_zbnum+1);
-	ele_num  = lm_xbnum * lm_ybnum * lm_zbnum;
+	ele_num_  = lm_xbnum * lm_ybnum * lm_zbnum;
 
 	for (int i = 0; i < lm_xbnum; i++)
 	{
@@ -65,8 +61,8 @@ void gctl::linear_mesh_3d::init(std::string in_name, std::string in_info, double
 		lm_zsizes[i] = zsizes[i];
 	}
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	int tmp_index;
 	double x_p, y_p, z_p = lm_zmin - 0.5*lm_zsizes[0];
@@ -112,16 +108,24 @@ void gctl::linear_mesh_3d::init(std::string in_name, std::string in_info, double
 		}
 	}
 
-	initialized = true;
+	initialized_ = true;
+	return;
+}
+
+void gctl::linear_mesh_3d::show_mesh_dimension(std::ostream &os) const
+{
+	os << "node num: " << node_num_ << std::endl;
+	os << "elem num: " << ele_num_ << std::endl;
+	os << "x-range: " << get_xmin() << " -> " << get_xmax() << "\n";
+	os << "y-range: " << get_ymin() << " -> " << get_ymax() << "\n";
+	os << "z-range: " << get_zmin() << " -> " << get_zmax() << "\n";
+	os << "dimension: " << lm_xbnum << ", " << lm_ybnum << ", " << lm_zbnum << "\n";
 	return;
 }
 
 void gctl::linear_mesh_3d::load_binary(std::string filename)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is already initialized.");
-	}
 
 	std::ifstream infile;
 	gctl::open_infile(infile, filename, ".2m", 
@@ -141,8 +145,8 @@ void gctl::linear_mesh_3d::load_binary(std::string filename)
 	lm_xsizes.resize(lm_xbnum);
 	lm_ysizes.resize(lm_ybnum);
 	lm_zsizes.resize(lm_zbnum);
-	node_num = (lm_xbnum+1) * (lm_ybnum+1) * (lm_zbnum+1);
+	node_num_ = (lm_xbnum+1) * (lm_ybnum+1) * (lm_zbnum+1);
-	ele_num  = lm_xbnum * lm_ybnum * lm_zbnum;
+	ele_num_  = lm_xbnum * lm_ybnum * lm_zbnum;
 
 	for (int i = 0; i < lm_xbnum; i++)
 	{
@@ -159,8 +163,8 @@ void gctl::linear_mesh_3d::load_binary(std::string filename)
 		infile.read((char*)lm_zsizes.get(i), sizeof(double));
 	}
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	int tmp_index;
 	double x_p, y_p, z_p = lm_zmin - 0.5*lm_zsizes[0];
@@ -206,7 +210,7 @@ void gctl::linear_mesh_3d::load_binary(std::string filename)
 		}
 	}
 
-	initialized = true;
+	initialized_ = true;
 
 	// 读入模型数据单元
 	load_datablock(infile);
@@ -217,10 +221,7 @@ void gctl::linear_mesh_3d::load_binary(std::string filename)
 
 void gctl::linear_mesh_3d::save_binary(std::string filename)
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
 
 	std::ofstream outfile;
 	gctl::open_outfile(outfile, filename, ".2m", 
@@ -272,91 +273,94 @@ gctl::linear_mesh_3d::~linear_mesh_3d(){}
 
 int gctl::linear_mesh_3d::get_xbnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return lm_xbnum;
 }
 
 int gctl::linear_mesh_3d::get_ybnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return lm_ybnum;
 }
 
 int gctl::linear_mesh_3d::get_zbnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return lm_zbnum;
 }
 
 double gctl::linear_mesh_3d::get_xmin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return lm_xmin;
 }
 
 double gctl::linear_mesh_3d::get_ymin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return lm_ymin;
 }
 
 double gctl::linear_mesh_3d::get_zmin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return lm_zmin;
 }
 
+double gctl::linear_mesh_3d::get_xmax() const
+{
+	check_initiated(); // 检查网格是否已经初始化
+
+	double xmax = lm_xmin - 0.5*lm_xsizes[0];
+	for (size_t i = 0; i < lm_xsizes.size(); i++)
+	{
+		xmax += lm_xsizes[i];
+	}
+	xmax -= 0.5*lm_xsizes.back();
+	return xmax;
+}
+
+double gctl::linear_mesh_3d::get_ymax() const
+{
+	check_initiated(); // 检查网格是否已经初始化
+
+	double ymax = lm_ymin - 0.5*lm_ysizes[0];
+	for (size_t i = 0; i < lm_ysizes.size(); i++)
+	{
+		ymax += lm_ysizes[i];
+	}
+	ymax -= 0.5*lm_ysizes.back();
+	return ymax;
+}
+
+double gctl::linear_mesh_3d::get_zmax() const
+{
+	check_initiated(); // 检查网格是否已经初始化
+
+	double zmax = lm_zmin - 0.5*lm_zsizes[0];
+	for (size_t i = 0; i < lm_zsizes.size(); i++)
+	{
+		zmax += lm_zsizes[i];
+	}
+	zmax -= 0.5*lm_zsizes.back();
+	return zmax;
+}
+
 const gctl::array<double> *gctl::linear_mesh_3d::get_xsizes() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return &lm_xsizes;
 }
 
 const gctl::array<double> *gctl::linear_mesh_3d::get_ysizes() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return &lm_ysizes;
 }
 
 const gctl::array<double> *gctl::linear_mesh_3d::get_zsizes() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::linear_mesh_3d] The mesh is not initialized.");
-	}
-
 	return &lm_zsizes;
 }
 
@@ -368,15 +372,3 @@ void gctl::linear_mesh_3d::save_gmsh(std::string filename, index_packed_e packed
 	outfile.close();
 	return;
 }
-
-void gctl::linear_mesh_3d::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
-	return;
-}
-
-void gctl::linear_mesh_3d::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, datname, out_mode, packed);
-	return;
-}

lib/mesh/linear_mesh_3d.h

@@ -40,13 +40,14 @@ namespace gctl
 		 * 以下为mesh类型的虚函数实现
 		 */
 
-		void init(std::string in_name, std::string in_info, double xmin, double ymin, 
+		void load_gmsh(std::string filename, index_packed_e packed = Packed){}
-			double zmin, const array<double> &xsizes, const array<double> &ysizes, 
+		void save_gmsh(std::string filename, index_packed_e packed = Packed);
-			const array<double> &zsizes);
 
 		void load_binary(std::string filename);
 		void save_binary(std::string filename);
 
+		void show_mesh_dimension(std::ostream &os) const;
+
 		/**
 		 * 以下为regular_mesh_2d的专有函数
 		 */
@@ -57,20 +58,23 @@ namespace gctl
 			const array<double> &zsizes);
 		virtual ~linear_mesh_3d();
 
+		void init(std::string in_name, std::string in_info, double xmin, double ymin, 
+			double zmin, const array<double> &xsizes, const array<double> &ysizes, 
+			const array<double> &zsizes);
+
 		int get_xbnum() const;
 		int get_ybnum() const;
 		int get_zbnum() const;
 		double get_xmin() const;
 		double get_ymin() const;
 		double get_zmin() const;
+		double get_xmax() const;
+		double get_ymax() const;
+		double get_zmax() const;
 		const array<double> *get_xsizes() const;
 		const array<double> *get_ysizes() const;
 		const array<double> *get_zsizes() const;
 
-		void save_gmsh(std::string filename, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
-
 	protected:
 		int lm_xbnum, lm_ybnum, lm_zbnum;
 		double lm_xmin, lm_ymin, lm_zmin;

lib/mesh/mesh.h

@@ -1,4 +1,4 @@
-/********************************************************
+/********************************************************
  *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
  * ██║  ███╗██║        ██║   ██║
@@ -28,13 +28,7 @@
 #ifndef _GCTL_BASE_MESH_H
 #define _GCTL_BASE_MESH_H
 
-#include "list"
-
 #include "meshdata.h"
-#include "meshdata_scalar.h"
-#include "meshdata_vector.h"
-#include "meshdata_tensor.h"
-
 #include "gctl/io.h"
 #include "gctl/algorithm.h"
 
@@ -42,6 +36,7 @@ namespace gctl
 {
 	enum mesh_type_e
 	{
+		UNDEFINED,
 		REGULAR_MESH,
 		LINEAR_MESH,
 		TRI_TET_MESH,
@@ -53,6 +48,7 @@ namespace gctl
 
 	enum mesh_dim_e
 	{
+		MESH_0D,
 		MESH_2D,
 		MESH_3D,
 	};
@@ -74,7 +70,6 @@ namespace gctl
 
 		/**
 		 * @brief 检查网格是否已经初始化。
-		 * 
 		 */
 		bool initiated() const;
 
@@ -85,36 +80,6 @@ namespace gctl
 		 */
 		bool saved(std::string datname) const;
 
-		/**
-		 * @brief 返回对应名称的数据对象指针。
-		 * 
-		 * @param datname 数据名
-		 * @return 数据指针
-		 */
-		meshdata *get_data(std::string datname) const;
-
-		/**
-		 * @brief 返回所有数据指针。
-		 * 
-		 * @param out_list 数据指针列表
-		 */
-		void get_all_data(array<meshdata*> &out_list) const;
-
-		/**
-		 * @brief 返回对应名称的数据对象的数组指针。
-		 * 
-		 * @param datname 数据名
-		 * @return 数组指针（注意需要转换为对应的数组指针类型）
-		 */
-		void *get_datval(std::string datname) const;
-
-		/**
-		 * @brief 删除对应名称的数据对象。
-		 * 
-		 * @param datname 数据名
-		 */
-		void remove_data(std::string datname);
-
 		/**
 		 * @brief 显示网格与数据信息。
 		 * 
@@ -122,14 +87,6 @@ namespace gctl
 		 */
 		void show_info(std::ostream &os = std::clog) const;
 
-		/**
-		 * @brief 重命名数据对象。
-		 * 
-		 * @param oldname 旧名称
-		 * @param newname 新名称
-		 */
-		void rename_data(std::string oldname, std::string newname);
-
 		/**
 		 * @brief 返回网格类型
 		 * 
@@ -194,156 +151,143 @@ namespace gctl
 		void set_meshinfo(std::string in_info);
 
 		/**
-		 * @brief 添加标量网格数据对象并赋初始值，如果对象已经存在则赋初始值
+		 * @brief 添加网格数据，并返回数据对象的引用。
 		 * 
-		 * @param in_name 数据对象名称
+		 * @param in_loctype 网格数据类型（块体或顶点数据）
-		 * @param in_type 数据对象类型
+		 * @param in_valtype 网格数据值类型
-		 * @param if_output 数据对象是否可输出
+		 * @param name 数据名称
-		 * @param init_val 初始值
+		 * @param init_val 数据初始值
-		 * @return 新建的数据对象指针
+		 * @param if_output 数据是否可用于输出
+		 * @param nan_val 数据无效值标记
+		 * 
+		 * @return 数据对象的引用（非常量的）
 		 */
-		meshdata *add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, double init_val);
+		meshdata &add_data(mesh_data_type_e in_loctype, mesh_data_value_e in_valtype, std::string name, 
+			double init_val, bool if_output = true, double nan_val = GCTL_BDL_MAX);
 
 		/**
-		 * @brief 添加适量网格数据对象并赋初始值，如果对象已经存在则赋初始值
+		 * @brief 添加网格数据，并返回数据对象的引用。
 		 *
-		 * @param in_name 数据对象名称
+		 * @param in_loctype 网格数据类型（块体或顶点数据）
-		 * @param in_type 数据对象类型
+		 * @param name 数据名称
-		 * @param if_output 数据对象是否可输出
+		 * @param init_arr 数据初始值数组
-		 * @param init_val 初始值
+		 * @param if_output 数据是否可用于输出
-		 * @return 新建的数据对象指针
+		 * @param nan_val 数据无效值标记
+		 * 
+		 * @return 数据对象的引用（非常量的）
 		 */
-		meshdata *add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, gctl::point3dc init_val);
+		meshdata &add_data(mesh_data_type_e in_loctype, std::string name, const array<double> &init_arr, 
+			bool if_output = true, double nan_val = GCTL_BDL_MAX);
 
 		/**
-		 * @brief 添加张量网格数据对象并赋初始值，如果对象已经存在则赋初始值
+		 * @brief 添加网格数据，并返回数据对象的引用。
 		 *
-		 * @param in_name 数据对象名称
+		 * @param in_loctype 网格数据类型（块体或顶点数据）
-		 * @param in_type 数据对象类型
+		 * @param name 数据名称
-		 * @param if_output 数据对象是否可输出
+		 * @param init_arr 数据初始值数组
-		 * @param init_val 初始值
+		 * @param if_output 数据是否可用于输出
-		 * @return 新建的数据对象指针
+		 * @param nan_val 数据无效值标记
+		 * 
+		 * @return 数据对象的引用（非常量的）
 		 */
-		meshdata *add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, gctl::tensor init_val);
+		meshdata &add_data(mesh_data_type_e in_loctype, std::string name, const array<point3dc> &init_arr, 
+			bool if_output = true, double nan_val = GCTL_BDL_MAX);
 
 		/**
-		 * @brief 添加张量网格数据对象并初始为0值
+		 * @brief 添加网格数据，并返回数据对象的引用。
 		 *
-		 * @param in_name 数据对象名称
+		 * @param in_loctype 网格数据类型（块体或顶点数据）
-		 * @param in_type 数据对象类型
+		 * @param name 数据名称
-		 * @param if_output 数据对象是否可输出
+		 * @param init_arr 数据初始值数组
-		 * @param init_val 初始值
+		 * @param if_output 数据是否可用于输出
-		 * @return 新建的数据对象指针
+		 * @param nan_val 数据无效值标记
+		 * 
+		 * @return 数据对象的引用（非常量的）
 		 */
-		meshdata *add_data(std::string in_name, mesh_data_type_e in_type, bool if_output, mesh_data_value_e val_type);
+		meshdata &add_data(mesh_data_type_e in_loctype, std::string name, const array<tensor> &init_arr, 
+			bool if_output = true, double nan_val = GCTL_BDL_MAX);
 
 		/**
-		 * @brief （基类虚函数）初始化规则网格
+		 * @brief 返回对应名称的数据对象引用。
 		 * 
-		 * @param in_name 网格名称
+		 * @param datname 数据名
-		 * @param in_info 网格说明信息
+		 * @return 数据对象的引用（非常量的）
-		 * @param xnum 网格x轴数量
-		 * @param ynum 网格y轴数量
-		 * @param xmin 网格x轴最小值
-		 * @param ymin 网格y轴最小值
-		 * @param dx 网格x轴间隔
-		 * @param dy 网格y轴间隔
 		 */
-		virtual void init(std::string in_name, std::string in_info, int xnum, int ynum, 
+		meshdata &get_data(std::string datname);
-			double xmin, double ymin, double dx, double dy);
 
 		/**
-		 * @brief 
+		 * @brief 返回对应名称的数据对象常量引用。
 		 * 
-		 * @param in_name 
+		 * @param datname 数据名
-		 * @param in_info 
+		 * @return 数据对象的引用（常量的）
-		 * @param xbnum 
-		 * @param ybnum 
-		 * @param zbnum 
-		 * @param xmin 
-		 * @param ymin 
-		 * @param zmin 
-		 * @param xsize 
-		 * @param ysize 
-		 * @param zsize 
 		 */
-		virtual void init(std::string in_name, std::string in_info, int xbnum, int ybnum, int zbnum, 
+		const meshdata &get_data(std::string datname) const;
-			double xmin, double ymin, double zmin, double xsize, double ysize, double zsize);
 
 		/**
-		 * @brief 
+		 * @brief 返回对应名称的数据对象指针。
 		 *
-		 * @param in_name 
+		 * @param datname 数据名
-		 * @param in_info 
+		 * @return 数据对象的指针
-		 * @param lon_min 
-		 * @param lat_min 
-		 * @param rad_min 
-		 * @param lon_size 
-		 * @param lat_size 
-		 * @param rad_size 
-		 * @param lon_bnum 
-		 * @param lat_bnum 
-		 * @param rad_bnum 
 		 */
-		virtual void init(std::string in_name, std::string in_info, double lon_min, double lat_min, 
+		meshdata *get_data_ptr(std::string datname);
-			double rad_min, double lon_size, double lat_size, double rad_size, int lon_bnum, int lat_bnum, int rad_bnum);
 
 		/**
-		 * @brief 
+		 * @brief 删除对应名称的数据对象。
 		 * 
-		 * @param in_name 
+		 * @param datname 数据名
-		 * @param in_info 
-		 * @param xmin 
-		 * @param ymin 
-		 * @param xsizes 
-		 * @param ysizes 
 		 */
-		virtual void init(std::string in_name, std::string in_info, double xmin, double ymin, 
+		void remove_data(std::string datname);
-			const gctl::array<double> &xsizes, const gctl::array<double> &ysizes);
 
 		/**
-		 * @brief 
+		 * @brief 保存网格到Gmsh文件，并保存所有数据。
 		 *
-		 * @param in_name 
+		 * @param filename 文件名
-		 * @param in_info 
+		 * @param out_mode 输出模式
-		 * @param xmin 
+		 * @param packed 索引是否为紧凑（从0开始）
-		 * @param ymin 
-		 * @param zmin 
-		 * @param xsizes 
-		 * @param ysizes 
-		 * @param zsizes 
 		 */
-		virtual void init(std::string in_name, std::string in_info, double xmin, double ymin, 
+		void save_gmsh_withdata(std::string filename, output_type_e out_mode, index_packed_e packed = Packed);
-			double zmin, const gctl::array<double> &xsizes, const gctl::array<double> &ysizes, 
-			const gctl::array<double> &zsizes);
 		
 		/**
-		 * @brief 
+		 * @brief 保存网格到Gmsh文件，并保存指定名称的数据。
 		 *
-		 * @param in_name 
+		 * @param filename 文件名
-		 * @param in_info 
+		 * @param datname 数据名
-		 * @param in_nodes 
+		 * @param out_mode 输出模式
-		 * @param in_triangles 
+		 * @param packed 索引是否为紧凑（从0开始）
 		 */
-		virtual void init(std::string in_name, std::string in_info, const gctl::array<gctl::vertex2dc> &in_nodes, 
+		void save_gmsh_withdata(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
-			const gctl::array<gctl::triangle2d> &in_triangles);
 
 		/**
-		 * @brief 
+		 * @brief 保存网格到Gmsh文件，并读入所有数据。
 		 *
-		 * @param in_name 
+		 * @param filename 文件名
-		 * @param in_info 
+		 * @param packed 索引是否为紧凑（从0开始）
-		 * @param in_nodes 
-		 * @param in_tets 
 		 */
-		virtual void init(std::string in_name, std::string in_info, const gctl::array<gctl::vertex3dc> &in_nodes, 
+		void load_gmsh_withdata(std::string filename, index_packed_e packed = Packed);
-			const gctl::array<gctl::tetrahedron> &in_tets);
 
 		/**
-		 * @brief 显示网格的维度信息
+		 * @brief 保存网格到Gmsh文件，并读入特定数据。
 		 *
+		 * @param filename 文件名
+		 * @param packed 索引是否为紧凑（从0开始）
 		 */
-		virtual void show_mesh_dimension(std::ostream &os) const;
+		void load_gmsh_withdata(std::string filename, std::string datname, mesh_data_type_e type, index_packed_e packed = Packed);
+
+		/**
+		 * @brief （基类纯虚函数）保存网格到Gmsh文件。
+		 * 
+		 * @param filename 文件名
+		 * @param packed 索引是否为紧凑（从0开始）
+		 */
+		virtual void save_gmsh(std::string filename, index_packed_e packed = Packed) = 0;
+
+		/**
+		 * @brief （基类纯虚函数）读入Gmsh文件。
+		 *
+		 * @param filename 文件名
+		 * @param packed 索引是否为紧凑（从0开始）
+		 */
+		virtual void load_gmsh(std::string filename, index_packed_e packed = Packed) = 0;
 
 		/**
 		 * @brief （基类纯虚函数）读入二进制网格文件
@@ -360,189 +304,49 @@ namespace gctl
 		virtual void save_binary(std::string filename) = 0;
 
 		/**
-		 * @brief （基类虚函数）保存网格到Gmsh文件。
+		 * @brief （基类纯虚函数）显示网格的维度信息
-		 * 
-		 * @param filename 文件名
-		 * @param packed 索引是否为紧凑（从0开始）
 		 */
-		virtual void save_gmsh(std::string filename, index_packed_e packed = Packed);
+		virtual void show_mesh_dimension(std::ostream &os) const = 0;
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param filename 
-		 * @param d_type 
-		 * @param out_mode 
-		 * @param packed 
-		 */
-		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param filename 
-		 * @param datname 
-		 * @param out_mode 
-		 * @param packed 
-		 */
-		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
-
-		/**
-		 * @brief      将输入的自由点位上的浮点数据插值为规则网络的顶点或者块体数据
-		 *
-		 * @param      in_posi      输入的自由点位的坐标
-		 * @param      in_val       输入的自由点位的数值
-		 * @param[in]  search_xlen  插值搜索区域的x半径
-		 * @param[in]  search_ylen  插值搜索区域的y半径
-		 * @param[in]  search_deg   插值搜索区域的x半径绕x轴正方向逆时针旋转的角度（弧度）
-		 * @param[in]  datname      网格后的数据名称
-		 * @param[in]  d_type       网格后的数据类型
-		 */
-		virtual void load_data_cloud(const array<point2dc> &in_posi, const array<double> &in_val, double search_xlen, 
-			double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type = NodeData);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param in_posi 
-		 * @param in_val 
-		 * @param search_xlen 
-		 * @param search_ylen 
-		 * @param search_deg 
-		 * @param datname 
-		 * @param d_type 
-		 */
-		virtual void load_data_cloud(const array<point3dc> &in_posi, const array<double> &in_val, double search_xlen, 
-			double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type = NodeData);
-
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param in_posi 
-		 * @param out_val 
-		 */
-		virtual void extract_points(std::string datname, const array<point2dc> &in_posi, array<double> &out_val);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param in_posi 
-		 * @param out_val 
-		 */
-		virtual void extract_points(std::string datname, const array<point3dc> &in_posi, array<double> &out_val);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param start_p 
-		 * @param end_p 
-		 * @param size_p 
-		 * @param out_posi 
-		 * @param out_val 
-		 */
-		virtual void extract_profile(std::string datname, const point2dc &start_p, const point2dc &end_p, int size_p, 
-			array<point2dc> &out_posi, array<double> &out_val);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param start_p 
-		 * @param end_p 
-		 * @param size_p 
-		 * @param dh 
-		 * @param out_posi 
-		 * @param out_val 
-		 */
-		virtual void extract_profile(std::string datname, const point3dc &start_p, const point3dc &end_p, int size_p, 
-			double dh, array<point3dc> &out_posi, array<double> &out_val);
-
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param p_type 
-		 * @param v_type 
-		 * @param para_str 
-		 * @param in_val 
-		 */
-		virtual void edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, double in_val);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param p_type 
-		 * @param v_type 
-		 * @param para_str 
-		 * @param in_val 
-		 */
-		virtual void edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, point3dc in_val);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param p_type 
-		 * @param v_type 
-		 * @param para_str 
-		 * @param in_val 
-		 */
-		virtual void edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, tensor in_val);
-
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param in_val 
-		 */
-		virtual void purge_data(std::string datname, double in_val);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param in_val 
-		 */
-		virtual void purge_data(std::string datname, point3dc in_val);
-		
-		/**
-		 * @brief 
-		 * 
-		 * @param datname 
-		 * @param in_val 
-		 */
-		virtual void purge_data(std::string datname, tensor in_val);
 
 	protected:
-		mesh_type_e meshtype;
+		mesh_type_e meshtype_;
-		mesh_dim_e  meshdim;
+		mesh_dim_e  meshdim_;
-		std::string meshname;
+		std::string meshname_;
-		std::string meshinfo;
+		std::string meshinfo_;
+		gmshio meshio_;
 
-		int node_num, ele_num;
+		int node_num_, ele_num_;
-		bool initialized;
+		bool initialized_;
 
-		std::list<meshdata*> saved_data;
+		std::vector<meshdata> datalist_;
-		std::list<meshdata*>::iterator iter;
 
 		/**
 		 * 以下为类的私有函数 只能被类的公共函数调用
 		 */
 
+		/**
+		 * @brief 检查网格是否已经初始化。
+		 * 
+		 */
+		void check_initiated(bool inverse = false) const;
+		
 		/**
 		 * @brief 初始化网格
 		 *
-		 * @param in_type 
+		 * @param in_type 网格类型
-		 * @param in_dim 
+		 * @param in_dim 网格维度
-		 * @param in_name 
+		 * @param in_name 网格名称
-		 * @param in_info 
+		 * @param in_info 网格说明
 		 */
 		void init(mesh_type_e in_type, mesh_dim_e in_dim, std::string in_name, std::string in_info);
 		
+		/**
+		 * @brief 返回数据索引并检查数据对象是否存在
+		 *
+		 * @param datname 数据名
+		 */
+		int data_index(std::string datname) const;
+		
 		/**
 		 * @brief 读入二进制网格头信息
 		 * 

lib/mesh/meshdata.cpp

@@ -1,4 +1,4 @@
-/********************************************************
+/********************************************************
  *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
  * ██║  ███╗██║        ██║   ██║
@@ -26,69 +26,239 @@
  ******************************************************/
 
 #include "meshdata.h"
+#include "fstream"
 
-gctl::meshdata::meshdata(std::string in_name, mesh_data_type_e in_type, bool if_output)
+gctl::meshdata::meshdata()
 {
-	if (in_name == "")
+    name_ = "untitled";
-	{
+    loctype_ = NodeData;
-		throw std::runtime_error("[gctl::meshdata] The input name is empty.");
+    valtype_ = Scalar;
-	}
+    output_ok_ = true;
-
+    nan_val_ = GCTL_BDL_MAX;
-	datname = in_name;
-	dattype = in_type;
-	output_ok = if_output;
 }
 
-gctl::meshdata::~meshdata(){}
+gctl::meshdata::meshdata(mesh_data_type_e in_loctype, 
-
+    mesh_data_value_e in_valtype, size_t size, 
-void gctl::meshdata::set_datname(std::string in_name)
+    std::string name, bool if_output, double nan_val)
 {
-	if (in_name == "")
+    create(in_loctype, in_valtype, size, name, if_output, nan_val);
-	{
+}
-		throw std::runtime_error("[gctl::meshdata] The input name is empty.");
-	}
 
-	datname = in_name;
+gctl::meshdata::~meshdata()
+{
+    clear();
+}
+
+void gctl::meshdata::clear()
+{
+    name_ = "untitled";
+    loctype_ = NodeData;
+    valtype_ = Scalar;
+    output_ok_ = true;
+    nan_val_ = GCTL_BDL_MAX;
+    datval_.clear();
+}
+
+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)
+{
+    name_ = name;
+    loctype_ = in_loctype;
+    valtype_ = in_valtype;
+    output_ok_ = if_output;
+    nan_val_ = nan_val;
+    if (size != 0) datval_.resize(size, 0.0);
+    return;
+}
+
+gctl::array<gctl::point3dc> gctl::meshdata::export_vector() const
+{
+    if (valtype_ != Vector) throw std::runtime_error("[gctl::meshdata::export_vector] Invalid data type.");
+
+    size_t n = 0;
+    array<point3dc> tmp(datval_.size()/3);
+    for (size_t i = 0; i < datval_.size(); i += 3)
+    {
+        tmp[n] = point3dc(datval_[3*i], datval_[3*i + 1], datval_[3*i + 2]);
+        n++;
+    }
+    return tmp;
+}
+
+gctl::array<gctl::tensor> gctl::meshdata::export_tensor() const
+{
+    if (valtype_ != Tensor) throw std::runtime_error("[gctl::meshdata::export_tensor] Invalid data type.");
+
+    size_t n = 0;
+    array<tensor> tmp(datval_.size()/9);
+    for (size_t i = 0; i < datval_.size(); i += 9)
+    {
+        tmp[n] = tensor(datval_[9*i], datval_[9*i + 1], datval_[9*i + 2],
+            datval_[9*i + 3], datval_[9*i + 4], datval_[9*i + 5],
+            datval_[9*i + 6], datval_[9*i + 7], datval_[9*i + 8]);
+        n++;
+    }
+    return tmp;
+}
+
+void gctl::meshdata::show_info(std::ostream &os) const
+{
+    os << "Data: " << name_ << " | ";
+	if (loctype_ == NodeData) os << "node data | ";
+	if (loctype_ == ElemData) os << "element data | ";
+	if (loctype_ == ElemData2D) os << "2D element data | ";
+	if (loctype_ == ElemData3D) os << "3D element data | ";
+	if (valtype_ == Scalar) os << "scalar | ";
+	if (valtype_ == Vector) os << "vector | ";
+	if (valtype_ == Tensor) os << "tensor | ";
+	if (output_ok_) os << "output" << std::endl;
+	else os<< "no-output" << std::endl;
 	return;
 }
 
-std::string gctl::meshdata::get_datname()
+void gctl::meshdata::show_stats(std::ostream &os) const
 {
-	return datname;
+    int vn = 0;
-}
+    if (valtype_ == Vector)
+    {
+        for (size_t i = 0; i < datval_.size(); i += 3)
+        {
+            if (fabs(datval_[3*i] - nan_val_) > 1e-10) vn++;
+        }
 
-gctl::mesh_data_type_e gctl::meshdata::get_dattype()
+        array<double> tmp_x(vn), tmp_y(vn), tmp_z(vn);
-{
-	return dattype;
-}
 
-gctl::mesh_data_value_e gctl::meshdata::get_valtype()
+        vn = 0;
-{
+        for (size_t i = 0; i < datval_.size(); i += 3)
-	return valtype;
+        {
-}
+            if (fabs(datval_[3*i] - nan_val_) > 1e-10)
+            {
+                tmp_x[vn] = datval_[3*i];
+                tmp_y[vn] = datval_[3*i + 1];
+                tmp_z[vn] = datval_[3*i + 2];
+                vn++;
+            }
+        }
         
-void gctl::meshdata::set_output(bool if_output)
+        os  << "max  = (" << tmp_x.max() << ", " << tmp_y.max() << ", " << tmp_z.max() << ")\n"
-{
+            << "min  = (" << tmp_x.min() << ", " << tmp_y.min() << ", " << tmp_z.min() << ")\n"
-	output_ok = if_output;
+            << "mean = (" << tmp_x.mean() << ", " << tmp_y.mean() << ", " << tmp_z.mean() << ")\n"
+            << "std  = (" << tmp_x.std() << ", " << tmp_y.std() << ", " << tmp_z.std() << ")\n"
+            << "rms  = (" << tmp_x.rms() << ", " << tmp_y.rms() << ", " << tmp_z.rms() << ")\n";
+    }
+    else if (valtype_ == Tensor)
+    {
+        for (size_t i = 0; i < datval_.size(); i += 9)
+        {
+            if (fabs(datval_[9*i] - nan_val_) > 1e-10) vn++;
+        }
+
+        array<double> tmp_xx(vn), tmp_xy(vn), tmp_xz(vn);
+        array<double> tmp_yx(vn), tmp_yy(vn), tmp_yz(vn);
+        array<double> tmp_zx(vn), tmp_zy(vn), tmp_zz(vn);
+
+        vn = 0;
+        for (size_t i = 0; i < datval_.size(); i += 9)
+        {
+            if (fabs(datval_[9*i] - nan_val_) > 1e-10)
+            {
+                tmp_xx[vn] = datval_[9*i];
+                tmp_xy[vn] = datval_[9*i + 1];
+                tmp_xz[vn] = datval_[9*i + 2];
+                tmp_yx[vn] = datval_[9*i + 3];
+                tmp_yy[vn] = datval_[9*i + 4];
+                tmp_yz[vn] = datval_[9*i + 5];
+                tmp_zx[vn] = datval_[9*i + 6];
+                tmp_zy[vn] = datval_[9*i + 7];
+                tmp_zz[vn] = datval_[9*i + 8];
+                vn++;
+            }
+        }
+
+        os << "max = \n"
+        << tmp_xx.max() << ", " << tmp_xy.max() << ", " << tmp_xz.max() << "\n"
+        << tmp_yx.max() << ", " << tmp_yy.max() << ", " << tmp_yz.max() << "\n"
+        << tmp_zx.max() << ", " << tmp_zy.max() << ", " << tmp_zz.max() << "\n";
+
+        os << "min = \n"
+        << tmp_xx.min() << ", " << tmp_xy.min() << ", " << tmp_xz.min() << "\n"
+        << tmp_yx.min() << ", " << tmp_yy.min() << ", " << tmp_yz.min() << "\n"
+        << tmp_zx.min() << ", " << tmp_zy.min() << ", " << tmp_zz.min() << "\n";
+
+        os << "mean = \n" 
+        << tmp_xx.mean() << ", " << tmp_xy.mean() << ", " << tmp_xz.mean() << "\n"
+        << tmp_yx.mean() << ", " << tmp_yy.mean() << ", " << tmp_yz.mean() << "\n"
+        << tmp_zx.mean() << ", " << tmp_zy.mean() << ", " << tmp_zz.mean() << "\n";
+
+        os << "std  = \n"
+        << tmp_xx.std() << ", " << tmp_xy.std() << ", " << tmp_xz.std() << "\n"
+        << tmp_yx.std() << ", " << tmp_yy.std() << ", " << tmp_yz.std() << "\n"
+        << tmp_zx.std() << ", " << tmp_zy.std() << ", " << tmp_zz.std() << "\n";
+        
+        os << "rms  = \n"
+        << tmp_xx.rms() << ", " << tmp_xy.rms() << ", " << tmp_xz.rms() << "\n"
+        << tmp_yx.rms() << ", " << tmp_yy.rms() << ", " << tmp_yz.rms() << "\n"
+        << tmp_zx.rms() << ", " << tmp_zy.rms() << ", " << tmp_zz.rms() << "\n";
+    }
+    else // Scalar
+    {
+        for (size_t i = 0; i < datval_.size(); i++)
+        {
+            if (fabs(datval_[i] - nan_val_) > 1e-10) vn++;
+        }
+
+        array<double> tmp(vn);
+
+        vn = 0;
+        for (size_t i = 0; i < datval_.size(); i++)
+        {
+            if (fabs(datval_[i] - nan_val_) > 1e-10)
+            {
+                tmp[vn] = datval_[i];
+                vn++;
+            }
+        }
+        
+        os << "max = " << tmp.max() << ", min = " << tmp.min() << std::endl;
+        os << "mean = " << tmp.mean() << ", std = " << tmp.std() << ", rms = " << tmp.rms() << std::endl;
+    }
 	return;
 }
 
-bool gctl::meshdata::get_output()
+void gctl::meshdata::load_binary(std::ifstream &infile)
 {
-	return output_ok;
+    int name_size;
+    infile.read((char*) &name_size, sizeof(int));
+
+    name_.resize(name_size);
+    infile.read((char*) name_.c_str(), name_size);
+
+    infile.read((char*) &loctype_, sizeof(int));
+	infile.read((char*) &valtype_, sizeof(int));
+    infile.read((char*) &output_ok_, sizeof(bool));
+
+    int n;
+    infile.read((char*) &n, sizeof(int));
+    datval_.resize(n);
+
+	infile.read((char*) datval_.get(), sizeof(double)*datval_.size());
+    infile.read((char*) &nan_val_, sizeof(double));
+    return;
 }
 
-void gctl::meshdata::show_info(std::ostream &os)
+void gctl::meshdata::save_binary(std::ofstream &outfile)
 {
-	os << "Data: " << datname << " | ";
+    int name_size = name_.size();
-	if (dattype == NodeData) os << "Type: Node data | ";
+	outfile.write((char*) &name_size, sizeof(int));
-	if (dattype == ElemData) os << "Type: Element data | ";
+	outfile.write((char*) name_.c_str(), name_size);
-	if (dattype == ElemData2D) os << "Type: 2D element data | ";
+	outfile.write((char*) &loctype_, sizeof(int));
-	if (dattype == ElemData3D) os << "Type: 3D element data | ";
+	outfile.write((char*) &valtype_, sizeof(int));
-	if (valtype == Scalar) os << "Value: Scalar | ";
+	outfile.write((char*) &output_ok_, sizeof(bool));
-	if (valtype == Vector) os << "Value: Vector | ";
+
-	if (valtype == Tensor) os << "Value: Tensor | ";
+    int n = datval_.size();
-	if (output_ok) os << "Output: Yes" << std::endl;
+	outfile.write((char*) &n, sizeof(int));
-	else os<< "Output: No" << std::endl;
+	outfile.write((char*) datval_.get(), sizeof(double)*datval_.size());
-	return;
+    outfile.write((char*) &nan_val_, sizeof(double));
+    return;
 }

lib/mesh/meshdata.h

@@ -1,4 +1,4 @@
-/********************************************************
+/********************************************************
  *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
  * ██║  ███╗██║        ██║   ██║
@@ -30,11 +30,12 @@
 
 #include "gctl_mesh_config.h"
 #include "gctl/core.h"
-#include "gctl/geometry.h"
+#include "gctl/geometry/point3c.h"
+#include "gctl/geometry/tensor.h"
 
 namespace gctl
 {
-	/**
+    /**
 	 * @brief      网格数据的数值类型
 	 */
 	enum mesh_data_value_e
@@ -44,93 +45,90 @@ namespace gctl
 		Tensor, ///< 张量数据
 	};
 
-	/**
+    /**
-	 * @brief      网格文件中的数据对象（纯虚类）。此类为具体类型的数据对象的基类
+     * @brief 网格文件中的数据对象，可储存标量、矢量和张量数据。
-	 */
+     * 
-	class meshdata
+     */
-	{
+    struct meshdata
-	protected:
+    {
-		std::string datname; // 数据的名称
+        std::string name_; // 数据的名称
-		mesh_data_type_e dattype; // 数据的赋值属性 顶点或是元素（设置后不可更改）
+        mesh_data_type_e loctype_; // 数据的赋值位置属性 顶点或是元素
-		mesh_data_value_e valtype; // 数据的类型（设置后不可更改）
+        mesh_data_value_e valtype_; // 数据的类型
-		bool output_ok; // 是否可输出数据
+        bool output_ok_; // 是否可输出数据
-
+        array<double> datval_; // 数据值 注意目前我们只支持浮点数据存储
-	public:
+        double nan_val_; // 无效数据的标记值
-		// 构造函数
-		meshdata(std::string in_name, mesh_data_type_e in_type, bool if_output);
-
-		// 析构函数
-		virtual ~meshdata();
 
 		/**
-		 * @brief      设置数据名称
+		 * @brief 构造函数
-		 *
-		 * @param[in]  in_name  名称字符串
 		 */
-		void set_datname(std::string in_name);
+		meshdata();
 
 		/**
-		 * @brief      返回数据的名称
+		 * @brief 创建数据对象
 		 * 
-		 * @return     名称字符串
+		 * @param in_loctype 数据的赋值位置属性 顶点或是元素
+		 * @param in_valtype 数据的类型
+		 * @param size 数据的大小
+		 * @param name 数据的名称
+		 * @param if_output 是否可输出数据
+		 * @param nan_val 无效数据的标记值
 		 */
-		std::string get_datname();
+        meshdata(mesh_data_type_e in_loctype, mesh_data_value_e in_valtype, 
+            size_t size, std::string name = "untitled", 
+			bool if_output = true, double nan_val = GCTL_BDL_MAX);
 		
 		/**
-		 * @brief      返回数据的赋值类型
+		 * @brief 析构函数
+		 */
+        virtual ~meshdata();
+
+		/**
+		 * @brief 清空数据对象
+		 */
+		void clear();
+
+		/**
+		 * @brief 创建数据对象
 		 * 
-		 * @return     数据赋值类型的枚举
+		 * @param in_loctype 数据的赋值位置属性 顶点或是元素
+		 * @param in_valtype 数据的类型
+		 * @param size 数据的大小
+		 * @param name 数据的名称
+		 * @param if_output 是否可输出数据
+		 * @param nan_val 无效数据的标记值
 		 */
-		mesh_data_type_e get_dattype();
+		void create(mesh_data_type_e in_loctype, mesh_data_value_e in_valtype,
+            size_t size, std::string name = "untitled", 
+			bool if_output = true, double nan_val = GCTL_BDL_MAX);
 
 		/**
-		 * @brief      返回数据值的类型
+		 * @brief 返回适量数据。
-		 *
-		 * @return     数据值类型的枚举
 		 */
-		mesh_data_value_e get_valtype();
+		array<point3dc> export_vector() const;
 
 		/**
-		 * @brief      设置数据对象输出状态的设置情况
+		 * @brief 返回张量数据。
-		 *
-		 * @param[in]  if_output  是否可输出此对象
 		 */
-		void set_output(bool if_output);
+		array<tensor> export_tensor() const;
 
-		/**
+        /**
-		 * @brief      返回数据对象输出状态的设置情况
-		 *
-		 * @return     是否可输出此对象
-		 */
-		bool get_output();
-
-		/**
 		 * @brief      显示数据对象的头信息
 		 */
-		void show_info(std::ostream &os = std::clog);
+		void show_info(std::ostream &os = std::clog) const;
 
-		/**
+        /**
 		 * @brief     显示数据统计参数
 		 */
-		virtual void show_stats(std::ostream &os = std::clog) = 0;
+		void show_stats(std::ostream &os = std::clog) const;
 
-		/**
+        /**
-		 * @brief      返回数据块的指针 所有子类的数据块均为gctl::array类型，这是一个通用的接口。
-		 * 
-		 * @note       注意因为具体的返回类型未知 因此需在调用后按实际数据类型转换指针
-		 *
-		 * @return     数据块的指针
-		 */
-		virtual void *get_datval_ptr() = 0;
-
-		/**
 		 * @brief      载入二进制文件
 		 * 
 		 * @warning    不要直接调用此函数，应该在网格类型的相应函数中调用
 		 *
 		 * @param      infile  输入文件的流（以二进制方式打开）
 		 */
-		virtual void load_binary(std::ifstream &infile) = 0;
+		void load_binary(std::ifstream &infile);
 
 		/**
 		 * @brief      保存二进制文件
@@ -139,20 +137,8 @@ namespace gctl
 		 *
 		 * @param      outfile  输出的文件流（以二进制方式打开）
 		 */
-		virtual void save_binary(std::ofstream &outfile) = 0;
+		void save_binary(std::ofstream &outfile);
-
+    };
-		/**
-		 * @brief      销毁指针的对象 可销毁此基类或子类指针指向的对象
-		 *
-		 * @param      obj_ptr  网格数据指针
-		 */
-		static void destroy(meshdata *obj_ptr)
-		{
-			delete obj_ptr;
-			obj_ptr = nullptr;
-			return;
-		}
-	};
 }
 
-#endif //_GCTL_MESHDATA_H
+#endif // _GCTL_MESHDATA_H

lib/mesh/regular_grid.h

@@ -49,90 +49,284 @@ namespace gctl
 		 * 以下为mesh类型的虚函数实现
 		 */
 
-		void init(std::string in_name, std::string in_info, int xnum, int ynum, 
+		void load_gmsh(std::string filename, index_packed_e packed = Packed){}
-			double xmin, double ymin, double dx, double dy);
 		
-		void show_mesh_dimension(std::ostream &os) const;
+		/**
+		 * @brief （基类纯虚函数的具体实现）保存网格到Gmsh文件。
+		 * 
+		 * @param filename 文件名
+		 * @param packed 索引是否为紧凑（从0开始）
+		 */
+		void save_gmsh(std::string filename, index_packed_e packed = Packed);
 
+		/**
+		 * @brief （基类纯虚函数的具体实现）读入二进制网格文件
+		 *
+		 * @param filename 文件名
+		 */
 		void load_binary(std::string filename);
+
+		/**
+		 * @brief （基类纯虚函数的具体实现）写入二进制网格文件
+		 *
+		 * @param filename 文件名
+		 */
 		void save_binary(std::string filename);
 
+		/**
+		 * @brief （基类纯虚函数的具体实现）显示网格的维度信息
+		 */
+		void show_mesh_dimension(std::ostream &os) const;
+
 		/**
 		 * 以下为regular_grid的专有函数
 		 */
 
-		regular_grid();
+		regular_grid(); ///< 构造函数
+
+		/**
+		 * @brief 构造函数
+		 *
+		 * @param in_name 网格名称
+		 * @param in_info 网格说明
+		 * @param xnum 网格x方向单元数量
+		 * @param ynum 网格y方向单元数量
+		 * @param xmin 网格x方向最小值
+		 * @param ymin 网格y方向最小值
+		 * @param dx 网格x方向单元长度
+		 * @param dy 网格y方向单元长度
+		 */
 		regular_grid(std::string in_name, std::string in_info, int xnum, int ynum, 
 			double xmin, double ymin, double dx, double dy);
-		virtual ~regular_grid();
 
+		virtual ~regular_grid(); ///< 析构函数
+
+		/**
+		 * @brief 初始化规则网格
+		 *
+		 * @param in_name 网格名称
+		 * @param in_info 网格说明
+		 * @param xnum 网格x方向单元数量
+		 * @param ynum 网格y方向单元数量
+		 * @param xmin 网格x方向最小值
+		 * @param ymin 网格y方向最小值
+		 * @param dx 网格x方向单元长度
+		 * @param dy 网格y方向单元长度
+		 */
+		void init(std::string in_name, std::string in_info, int xnum, int ynum, 
+			double xmin, double ymin, double dx, double dy);
+		
+		/**
+		 * @brief 清除规则网格
+		 */
 		void clear_regular_grid();
+
+		/**
+		 * @brief 使用mathgl显示规则网格
+		 * 
+		 * @param datname 数据名称
+		 * @return 状态代码
+		 */
 		int view(std::string datname);
+
+		/**
+		 * @brief 使用gmt绘制规则网格
+		 *
+		 * @param datname 数据名称
+		 */
 		void plot(std::string datname);
 
-		int get_xdim() const;
+		int get_xdim() const; ///< 获取规则网格的x方向单元数量
-		int get_ydim() const;
+		int get_ydim() const; ///< 获取规则网格的y方向单元数量
-		double get_xmin() const;
+		double get_xmin() const; ///< 获取规则网格的x方向最小值
-		double get_ymin() const;
+		double get_ymin() const; ///< 获取规则网格的y方向最小值
-		double get_dx() const;
+		double get_dx() const; ///< 获取规则网格的x方向单元长度
-		double get_dy() const;
+		double get_dy() const; ///< 获取规则网格的y方向单元长度
 
 #ifdef GCTL_NETCDF
 
+		/**
+		 * @brief 从netcdf文件中读取规则网格数据
+		 *
+		 * @param filename 文件名
+		 * @param d_type 数据类型
+		 * @param xname x坐标变量名
+		 * @param yname y坐标变量名
+		 * @param zname z坐标变量名
+		 */
 		void load_netcdf_grid(std::string filename, mesh_data_type_e d_type, 
 			std::string xname = "x", std::string yname = "y", std::string zname = "null");
+
+		/**
+		 * @brief 保存规则网格数据到netcdf文件中
+		 *
+		 * @param filename 文件名
+		 * @param d_type 数据类型
+		 * @param xname x坐标变量名
+		 * @param yname y坐标变量名
+		 */
 		void save_netcdf_grid(std::string filename, mesh_data_type_e d_type, 
 			std::string xname = "x", std::string yname = "y");
+		
+		/**
+		 * @brief 保存规则网格数据到netcdf文件中
+		 *
+		 * @param filename 文件名
+		 * @param datname 数据名称
+		 * @param xname x坐标变量名
+		 * @param yname y坐标变量名
+		 */
 		void save_netcdf_grid(std::string filename, std::string datname, 
 			std::string xname = "x", std::string yname = "y");
 
 #endif // GCTL_NETCDF
 
+		/**
+		 * @brief 从surfer文件中读取规则网格数据
+		 *
+		 * @param filename 文件名
+		 * @param datname 数据名称
+		 * @param d_type 数据类型
+		 * @param grid_type 网格类型
+		 */
 		void load_surfer_grid(std::string filename, std::string datname, mesh_data_type_e d_type, surfer_file_type_e grid_type = Surfer7Grid);
+		
+		/**
+		 * @brief 保存规则网格数据到surfer文件中
+		 *
+		 * @param filename 文件名
+		 * @param datname 数据名称
+		 * @param grid_type 网格类型
+		 */
 		void save_surfer_grid(std::string filename, std::string datname, surfer_file_type_e grid_type = Surfer7Grid);
 
-		void save_gmsh(std::string filename, index_packed_e packed = Packed);
+		/**
-		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
+		 * @brief 保存规则网格数据到文本文件中
-		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
+		 *
-
+		 * @param filename 文件名
+		 * @param datname 数据名称
+		 */
 		void save_text(std::string filename, const array<std::string> &datname);
 
+		/**
+		 * @brief 从点云数据中提取规则网格数据
+		 *
+		 * @param in_posi 点云位置数组
+		 * @param in_val 点云数据数组
+		 * @param search_xlen 搜索范围x方向长度
+		 * @param search_ylen 搜索范围y方向长度
+		 * @param search_deg 搜索范围角度
+		 * @param datname 新建数据名称
+		 * @param d_type 新建数据类型
+		 */
 		void load_data_cloud(const array<point2dc> &in_posi, const array<double> &in_val, double search_xlen, 
 			double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type = NodeData);
 
+		/**
+		 * @brief 从规则网格数据中提取点云数据
+		 *
+		 * @param datname 数据名称
+		 * @param in_posi 点云位置数组
+		 * @param out_val 输出的点云数据数组
+		 */
 		void extract_points(std::string datname, const array<point2dc> &in_posi, array<double> &out_val);
+
+		/**
+		 * @brief 从规则网格数据中提取剖面数据
+		 *
+		 * @param datname 数据名称
+		 * @param start_p 剖面起点
+		 * @param end_p 剖面终点
+		 * @param size_p 剖面点数
+		 * @param out_posi 输出的剖面位置数组
+		 * @param out_val 输出的剖面数据数组
+		 */
 		void extract_profile(std::string datname, const point2dc &start_p, const point2dc &end_p, int size_p, 
 			array<point2dc> &out_posi, array<double> &out_val);
 
+		/**
+		 * @brief 计算规则网格数据的梯度数据
+		 *
+		 * @param datname 数据名称
+		 * @param gradname 梯度数据名称
+		 * @param d_type 梯度数据类型
+		 * @param order 梯度阶数
+		 */
 		void gradient(std::string datname, std::string gradname, gradient_type_e d_type, int order = 1);
+		
+		/**
+		 * @brief 计算两个规则网格数据的和
+		 *
+		 * @param newname 新数据名称
+		 * @param datname 数据名称
+		 * @param datname2 数据名称2
+		 */
 		void sum(std::string newname, std::string datname, std::string datname2);
+
+		/**
+		 * @brief 计算两个规则网格数据的差
+		 *
+		 * @param newname 新数据名称
+		 * @param datname 数据名称
+		 * @param datname2 数据名称2
+		 */
 		void diff(std::string newname, std::string datname, std::string datname2);
+
+		/**
+		 * @brief 计算规则网格数据的布尔数据
+		 *
+		 * @param newname 新数据名称
+		 * @param datname 数据名称
+		 * @param maskname 掩码数据名称
+		 * @param reverse 是否反转掩码
+		 */
 		void boolean(std::string newname, std::string datname, std::string maskname, bool reverse = false);
 
 #ifdef GCTL_MESH_EXPRTK
+
+		/**
+		 * @brief 计算规则网格数据的函数值，其中f为节点位置（x，y）的函数
+		 *
+		 * @param expression_str 表达式字符串
+		 * @param newname 新数据名称
+		 * @param x_str x变量名称
+		 * @param y_str y变量名称
+		 * @param f_str f变量名称
+		 */
 		void function(std::string expression_str, std::string newname, std::string x_str = "x", std::string y_str = "y", std::string f_str = "f");
+		
+		/**
+		 * @brief 计算规则网格数据的函数值，计算不同网格数据合成的函数值
+		 *
+		 * @param expression_str 表达式字符串
+		 * @param var_list 变量名称列表
+		 * @param data_list 数据名称列表
+		 */
 		void calculator(std::string expression_str, const array<std::string> &var_list, const array<std::string> &data_list);
+
 #endif // GCTL_MESH_EXPRTK
 
 #ifdef GCTL_MESH_WAVELIB
+
+		/**
+		 * @brief 计算规则网格数据的小波数据
+		 *
+		 * @param datname 数据名称
+		 * @param wavename 小波数据名称
+		 * @param order 小波阶数
+		 * @param summary 是否输出小波处理摘要
+		 */
 		void wavelet(std::string datname, std::string wavename, int order, bool summary = true);
+
 #endif // GCTL_MESH_WAVELIB
 
 	protected:
-		int rg_xnum, rg_ynum;
+		int rg_xnum, rg_ynum; ///< 规则网格的x方向单元数量和y方向单元数量
-		double rg_xmin, rg_ymin;
+		double rg_xmin, rg_ymin; ///< 规则网格的x方向最小值和y方向最小值
-		double rg_dx, rg_dy;
+		double rg_dx, rg_dy; ///< 规则网格的x方向单元长度和y方向单元长度
 
-#ifdef GCTL_GRAPHIC_MATHGL
+		array<vertex2dc> nodes; ///< 规则网格的节点数组
-		mathgl_dens plt;
+		array<rectangle2d> elements; ///< 规则网格的单元数组
-#endif // GCTL_GRAPHIC_MATHGL
-
-#ifdef GCTL_GRAPHIC_GMT
-		gmt_JX_single pic;
-#endif // GCTL_GRAPHIC_GMT
-
-		array<vertex2dc> nodes;
-		array<rectangle2d> elements;
 	};
 }
 

lib/mesh/regular_mesh_2d.cpp

@@ -30,22 +30,18 @@
 void gctl::regular_mesh_2d::init(std::string in_name, std::string in_info, int xbnum, int ybnum, 
 	double xmin, double ymin, double xsize, double ysize)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is already initialized.");
-	}
-
 	base_mesh::init(REGULAR_MESH, MESH_2D, in_name, in_info);
 
 	rm_xbnum = xbnum; rm_ybnum = ybnum;
 	rm_xmin = xmin; rm_ymin = ymin;
 	rm_xsize = xsize; rm_ysize = ysize;
 
-	node_num = (rm_xbnum+1) * (rm_ybnum+1);
+	node_num_ = (rm_xbnum+1) * (rm_ybnum+1);
-	ele_num  = rm_xbnum * rm_ybnum;
+	ele_num_  = rm_xbnum * rm_ybnum;
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	for (int j = 0; j < rm_ybnum+1; j++)
 	{
@@ -71,16 +67,24 @@ void gctl::regular_mesh_2d::init(std::string in_name, std::string in_info, int x
 		}
 	}
 
-	initialized = true;
+	initialized_ = true;
+	return;
+}
+
+void gctl::regular_mesh_2d::show_mesh_dimension(std::ostream &os) const
+{
+	os << "node num: " << node_num_ << std::endl;
+	os << "elem num: " << ele_num_ << std::endl;
+	os << "x-range: " << rm_xmin << " -> " << rm_xmin + (rm_xbnum - 1)*rm_xsize << "\n";
+	os << "y-range: " << rm_ymin << " -> " << rm_ymin + (rm_ybnum - 1)*rm_ysize << "\n";
+	os << "block-size: " << rm_xsize << ", " << rm_ysize << "\n";
+	os << "dimension: " << rm_xbnum << ", " << rm_ybnum << "\n";
 	return;
 }
 
 void gctl::regular_mesh_2d::load_binary(std::string filename)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化(如果已经初始化，那么就会抛出异常)
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is already initialized.");
-	}
 
 	std::ifstream infile;
 	gctl::open_infile(infile, filename, ".2m", 
@@ -97,11 +101,11 @@ void gctl::regular_mesh_2d::load_binary(std::string filename)
 	infile.read((char*)&rm_xsize, sizeof(double));
 	infile.read((char*)&rm_ysize, sizeof(double));
 
-	node_num = (rm_xbnum+1) * (rm_ybnum+1);
+	node_num_ = (rm_xbnum+1) * (rm_ybnum+1);
-	ele_num  = rm_xbnum * rm_ybnum;
+	ele_num_  = rm_xbnum * rm_ybnum;
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	for (int j = 0; j < rm_ybnum+1; j++)
 	{
@@ -127,21 +131,17 @@ void gctl::regular_mesh_2d::load_binary(std::string filename)
 		}
 	}
 
-	initialized = true;
-
 	// 读入模型数据单元
 	load_datablock(infile);
 
 	infile.close();
+	initialized_ = true;
 	return;
 }
 
 void gctl::regular_mesh_2d::save_binary(std::string filename)
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is not initialized.");
-	}
 
 	std::ofstream outfile;
 	gctl::open_outfile(outfile, filename, ".2m", 
@@ -177,61 +177,37 @@ gctl::regular_mesh_2d::~regular_mesh_2d(){}
 
 int gctl::regular_mesh_2d::get_xbnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is not initialized.");
-	}
-
 	return rm_xbnum;
 }
 
 int gctl::regular_mesh_2d::get_ybnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is not initialized.");
-	}
-
 	return rm_ybnum;
 }
 
 double gctl::regular_mesh_2d::get_xmin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is not initialized.");
-	}
-
 	return rm_xmin;
 }
 
 double gctl::regular_mesh_2d::get_ymin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is not initialized.");
-	}
-
 	return rm_ymin;
 }
 
 double gctl::regular_mesh_2d::get_xsize() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is not initialized.");
-	}
-
 	return rm_xsize;
 }
 
 double gctl::regular_mesh_2d::get_ysize() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_2d] The mesh is not initialized.");
-	}
-
 	return rm_ysize;
 }
 
@@ -243,15 +219,3 @@ void gctl::regular_mesh_2d::save_gmsh(std::string filename, index_packed_e packe
 	outfile.close();
 	return;
 }
-
-void gctl::regular_mesh_2d::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
-	return;
-}
-
-void gctl::regular_mesh_2d::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, datname, out_mode, packed);
-	return;
-}

lib/mesh/regular_mesh_2d.h

@@ -40,12 +40,14 @@ namespace gctl
 		 * 以下为mesh类型的虚函数实现
 		 */
 
-		void init(std::string in_name, std::string in_info, int xbnum, int ybnum, 
+		void load_gmsh(std::string filename, index_packed_e packed = Packed){}
-			double xmin, double ymin, double xsize, double ysize);
+		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;
+
 		/**
 		 * 以下为regular_mesh_2d的专有函数
 		 */
@@ -55,6 +57,9 @@ namespace gctl
 			double xmin, double ymin, double xsize, double ysize);
 		virtual ~regular_mesh_2d();
 
+		void init(std::string in_name, std::string in_info, int xbnum, int ybnum, 
+			double xmin, double ymin, double xsize, double ysize);
+
 		int get_xbnum() const;
 		int get_ybnum() const;
 		double get_xmin() const;
@@ -62,10 +67,6 @@ namespace gctl
 		double get_xsize() const;
 		double get_ysize() const;
 
-		void save_gmsh(std::string filename, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
-
 	protected:
 		int rm_xbnum, rm_ybnum;
 		double rm_xmin, rm_ymin;

lib/mesh/regular_mesh_3d.cpp

@@ -30,22 +30,18 @@
 void gctl::regular_mesh_3d::init(std::string in_name, std::string in_info, int xbnum, int ybnum, int zbnum, 
 	double xmin, double ymin, double zmin, double xsize, double ysize, double zsize)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is already initialized.");
-	}
-
 	base_mesh::init(REGULAR_MESH, MESH_3D, in_name, in_info);
 
 	rm_xbnum = xbnum; rm_ybnum = ybnum; rm_zbnum = zbnum;
 	rm_xmin = xmin; rm_ymin = ymin; rm_zmin = zmin;
 	rm_xsize = xsize; rm_ysize = ysize; rm_zsize = zsize;
 
-	node_num = (rm_xbnum+1) * (rm_ybnum+1) *(rm_zbnum+1);
+	node_num_ = (rm_xbnum+1) * (rm_ybnum+1) *(rm_zbnum+1);
-	ele_num  = rm_xbnum * rm_ybnum * rm_zbnum;
+	ele_num_  = rm_xbnum * rm_ybnum * rm_zbnum;
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	for (int k = 0; k < rm_zbnum+1; k++)
 	{
@@ -82,16 +78,25 @@ void gctl::regular_mesh_3d::init(std::string in_name, std::string in_info, int x
 		}
 	}
 
-	initialized = true;
+	initialized_ = true;
+	return;
+}
+
+void gctl::regular_mesh_3d::show_mesh_dimension(std::ostream &os) const
+{
+	os << "node num: " << node_num_ << std::endl;
+	os << "elem num: " << ele_num_ << std::endl;
+	os << "x-range: " << rm_xmin << " -> " << rm_xmin + (rm_xbnum - 1)*rm_xsize << "\n";
+	os << "y-range: " << rm_ymin << " -> " << rm_ymin + (rm_ybnum - 1)*rm_ysize << "\n";
+	os << "z-range: " << rm_zmin << " -> " << rm_zmin + (rm_zbnum - 1)*rm_zsize << "\n";
+	os << "block-size: " << rm_xsize << ", " << rm_ysize << ", " << rm_zsize << "\n";
+	os << "dimension: " << rm_xbnum << ", " << rm_ybnum << ", " << rm_zbnum << "\n";
 	return;
 }
 
 void gctl::regular_mesh_3d::load_binary(std::string filename)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is already initialized.");
-	}
 
 	std::ifstream infile;
 	gctl::open_infile(infile, filename, ".2m", 
@@ -111,11 +116,11 @@ void gctl::regular_mesh_3d::load_binary(std::string filename)
 	infile.read((char*)&rm_ysize, sizeof(double));
 	infile.read((char*)&rm_zsize, sizeof(double));
 
-	node_num = (rm_xbnum+1) * (rm_ybnum+1) *(rm_zbnum+1);
+	node_num_ = (rm_xbnum+1) * (rm_ybnum+1) *(rm_zbnum+1);
-	ele_num  = rm_xbnum * rm_ybnum * rm_zbnum;
+	ele_num_  = rm_xbnum * rm_ybnum * rm_zbnum;
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	for (int k = 0; k < rm_zbnum+1; k++)
 	{
@@ -152,21 +157,17 @@ void gctl::regular_mesh_3d::load_binary(std::string filename)
 		}
 	}
 
-	initialized = true;
-
 	// 读入模型数据单元
 	load_datablock(infile);
 
 	infile.close();
+	initialized_ = true;
 	return;
 }
 
 void gctl::regular_mesh_3d::save_binary(std::string filename)
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
 
 	std::ofstream outfile;
 	gctl::open_outfile(outfile, filename, ".2m", 
@@ -205,91 +206,55 @@ gctl::regular_mesh_3d::~regular_mesh_3d(){}
 
 int gctl::regular_mesh_3d::get_xbnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_xbnum;
 }
 
 int gctl::regular_mesh_3d::get_ybnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_ybnum;
 }
 
 int gctl::regular_mesh_3d::get_zbnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_zbnum;
 }
 
 double gctl::regular_mesh_3d::get_xmin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_xmin;
 }
 
 double gctl::regular_mesh_3d::get_ymin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_ymin;
 }
 
 double gctl::regular_mesh_3d::get_zmin() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_zmin;
 }
 
 double gctl::regular_mesh_3d::get_xsize() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_xsize;
 }
 
 double gctl::regular_mesh_3d::get_ysize() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_ysize;
 }
 
 double gctl::regular_mesh_3d::get_zsize() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_3d] The mesh is not initialized.");
-	}
-
 	return rm_zsize;
 }
 
@@ -301,43 +266,3 @@ void gctl::regular_mesh_3d::save_gmsh(std::string filename, index_packed_e packe
 	outfile.close();
 	return;
 }
-
-void gctl::regular_mesh_3d::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
-	return;
-}
-
-void gctl::regular_mesh_3d::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, datname, out_mode, packed);
-	return;
-}
-
-/*
-void gctl::regular_mesh_3d::edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, double in_val)
-{
-	std::string err_str;
-	meshdata *dat_ptr = get_data(datname);
-	if (dat_ptr->get_valtype() != Double)
-	{
-		err_str = "Incompatible value type of data: " << datname << ". From regular_mesh_3d::edit_data(...)";
-		throw runtime_error(err_str);
-	}
-
-	array<double>* val_ptr = (array<double>*) dat_ptr->get_datval();
-	// 按形体类型解析参数字符串
-	if (p_type == Block)
-	{
-		double xmin, xmax, ymin, ymax, zmin, zmax;
-		block_parameter_parser(para_str, xmin, xmax, ymin, ymax, zmin, zmax);
-
-
-	}
-	else
-	{
-		err_str = "Invalid physical entity type. From regular_mesh_3d::edit_data(...)";
-		throw runtime_error(err_str);
-	}
-}
-*/

lib/mesh/regular_mesh_3d.h

@@ -40,12 +40,14 @@ namespace gctl
 		 * 以下为mesh类型的虚函数实现
 		 */
 
-		void init(std::string in_name, std::string in_info, int xbnum, int ybnum, int zbnum, 
+		void load_gmsh(std::string filename, index_packed_e packed = Packed){}
-			double xmin, double ymin, double zmin, double xsize, double ysize, double zsize);
+		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;
+
 		/**
 		 * 以下为regular_mesh_3d的专有函数
 		 */
@@ -55,6 +57,9 @@ namespace gctl
 			double xmin, double ymin, double zmin, double xsize, double ysize, double zsize);
 		virtual ~regular_mesh_3d();
 
+		void init(std::string in_name, std::string in_info, int xbnum, int ybnum, int zbnum, 
+			double xmin, double ymin, double zmin, double xsize, double ysize, double zsize);
+
 		int get_xbnum() const;
 		int get_ybnum() const;
 		int get_zbnum() const;
@@ -65,13 +70,6 @@ namespace gctl
 		double get_ysize() const;
 		double get_zsize() const;
 
-		void save_gmsh(std::string filename, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
-
-		//void edit_data(std::string datname, physical_type_e p_type, value_operator_e v_type, std::string para_str, double in_val);
-		//void purge_data(std::string datname, double in_val);
-
 	protected:
 		int rm_xbnum, rm_ybnum, rm_zbnum;
 		double rm_xmin, rm_ymin, rm_zmin;

lib/mesh/regular_mesh_sph_3d.cpp

@@ -30,22 +30,18 @@
 void gctl::regular_mesh_sph_3d::init(std::string in_name, std::string in_info, double lon_min, double lat_min, 
 	double rad_min, double lon_size, double lat_size, double rad_size, int lon_bnum, int lat_bnum, int rad_bnum)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is already initialized.");
-	}
-
 	base_mesh::init(REGULAR_MESH_SPH, MESH_3D, in_name, in_info);
 
 	rm_lon_bnum = lon_bnum; rm_lat_bnum = lat_bnum; rm_rad_bnum = rad_bnum;
 	rm_lon_min = lon_min; rm_lat_min = lat_min; rm_rad_min = rad_min;
 	rm_lon_size = lon_size; rm_lat_size = lat_size; rm_rad_size = rad_size;
 
-	node_num = (rm_lon_bnum+1) * (rm_lat_bnum+1) *(rm_rad_bnum+1);
+	node_num_ = (rm_lon_bnum+1) * (rm_lat_bnum+1) *(rm_rad_bnum+1);
-	ele_num  = rm_lon_bnum * rm_lat_bnum * rm_rad_bnum;
+	ele_num_  = rm_lon_bnum * rm_lat_bnum * rm_rad_bnum;
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	int tmp_id;
 	for (int k = 0; k < rm_rad_bnum+1; k++)
@@ -85,16 +81,25 @@ void gctl::regular_mesh_sph_3d::init(std::string in_name, std::string in_info, d
 		}
 	}
 
-	initialized = true;
+	initialized_ = true;
+	return;
+}
+
+void gctl::regular_mesh_sph_3d::show_mesh_dimension(std::ostream &os) const
+{
+	os << "node num: " << node_num_ << std::endl;
+	os << "elem num: " << ele_num_ << std::endl;
+	os << "x-range: " << rm_lon_min << " -> " << rm_lon_min + (rm_lon_bnum - 1)*rm_lon_size << "\n";
+	os << "y-range: " << rm_lat_min << " -> " << rm_lat_min + (rm_lat_bnum - 1)*rm_lat_size << "\n";
+	os << "z-range: " << rm_rad_min << " -> " << rm_rad_min + (rm_rad_bnum - 1)*rm_rad_size << "\n";
+	os << "block-size: " << rm_lon_size << ", " << rm_lat_size << ", " << rm_rad_size << "\n";
+	os << "dimension: " << rm_lon_bnum << ", " << rm_lat_bnum << ", " << rm_rad_bnum << "\n";
 	return;
 }
 
 void gctl::regular_mesh_sph_3d::load_binary(std::string filename)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is already initialized.");
-	}
 
 	std::ifstream infile;
 	gctl::open_infile(infile, filename, ".2m", 
@@ -114,11 +119,11 @@ void gctl::regular_mesh_sph_3d::load_binary(std::string filename)
 	infile.read((char*)&rm_lat_size, sizeof(double));
 	infile.read((char*)&rm_rad_size, sizeof(double));
 
-	node_num = (rm_lon_bnum+1) * (rm_lat_bnum+1) *(rm_rad_bnum+1);
+	node_num_ = (rm_lon_bnum+1) * (rm_lat_bnum+1) *(rm_rad_bnum+1);
-	ele_num  = rm_lon_bnum * rm_lat_bnum * rm_rad_bnum;
+	ele_num_  = rm_lon_bnum * rm_lat_bnum * rm_rad_bnum;
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
 	int tmp_id;
 	for (int k = 0; k < rm_rad_bnum+1; k++)
@@ -158,7 +163,7 @@ void gctl::regular_mesh_sph_3d::load_binary(std::string filename)
 		}
 	}
 
-	initialized = true;
+	initialized_ = true;
 
 	// 读入模型数据单元
 	load_datablock(infile);
@@ -169,10 +174,7 @@ void gctl::regular_mesh_sph_3d::load_binary(std::string filename)
 
 void gctl::regular_mesh_sph_3d::save_binary(std::string filename)
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
 
 	std::ofstream outfile;
 	gctl::open_outfile(outfile, filename, ".2m", 
@@ -211,91 +213,55 @@ gctl::regular_mesh_sph_3d::~regular_mesh_sph_3d(){}
 
 int gctl::regular_mesh_sph_3d::get_lon_bnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_lon_bnum;
 }
 
 int gctl::regular_mesh_sph_3d::get_lat_bnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_lat_bnum;
 }
 
 int gctl::regular_mesh_sph_3d::get_rad_bnum() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_rad_bnum;
 }
 
 double gctl::regular_mesh_sph_3d::get_lon_min() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_lon_min;
 }
 
 double gctl::regular_mesh_sph_3d::get_lat_min() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_lat_min;
 }
 
 double gctl::regular_mesh_sph_3d::get_rad_min() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_rad_min;
 }
 
 double gctl::regular_mesh_sph_3d::get_lon_size() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_lon_size;
 }
 
 double gctl::regular_mesh_sph_3d::get_lat_size() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_lat_size;
 }
 
 double gctl::regular_mesh_sph_3d::get_rad_size() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::regular_mesh_sph_3d] The mesh is not initialized.");
-	}
-
 	return rm_rad_size;
 }
 
@@ -318,15 +284,3 @@ void gctl::regular_mesh_sph_3d::save_gmsh(std::string filename, index_packed_e p
 	tmp_nodes.clear();
 	return;
 }
-
-void gctl::regular_mesh_sph_3d::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
-	return;
-}
-
-void gctl::regular_mesh_sph_3d::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, datname, out_mode, packed);
-	return;
-}

lib/mesh/regular_mesh_sph_3d.h

@@ -39,11 +39,14 @@ namespace gctl
         /**
          * 以下为mesh类型的虚函数实现
          */
-        void init(std::string in_name, std::string in_info, double lon_min, double lat_min, double rad_min, 
-            double lon_size, double lat_size, double rad_size, int lon_bnum, int lat_bnum, int rad_bnum);
 
-        void load_binary(std::string filename);
+        void load_gmsh(std::string filename, index_packed_e packed = Packed){}
-        void save_binary(std::string filename);
+		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;
 
         /**
          * 以下为regular_mesh_sph_3d的专有函数
@@ -53,6 +56,9 @@ namespace gctl
             double rad_min, double lon_size, double lat_size, double rad_size, int lon_bnum, int lat_bnum, int rad_bnum);
         virtual ~regular_mesh_sph_3d();
 
+        void init(std::string in_name, std::string in_info, double lon_min, double lat_min, double rad_min, 
+            double lon_size, double lat_size, double rad_size, int lon_bnum, int lat_bnum, int rad_bnum);
+
         int get_lon_bnum() const;
         int get_lat_bnum() const;
         int get_rad_bnum() const;
@@ -63,10 +69,6 @@ namespace gctl
         double get_lat_size() const;
         double get_rad_size() const;
 
-        void save_gmsh(std::string filename, index_packed_e packed = Packed);
-        void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
-        void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
-
     protected:
         int rm_lon_bnum, rm_lat_bnum, rm_rad_bnum;
         double rm_lon_min, rm_lat_min, rm_rad_min;

lib/mesh/tet_mesh.cpp

@@ -30,24 +30,20 @@
 void gctl::tetrahedron_mesh::init(std::string in_name, std::string in_info, const array<vertex3dc> &in_nodes, 
 	const array<tetrahedron> &in_tets)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::tetrahedron_mesh] The mesh is already initialized.");
-	}
-
 	base_mesh::init(TRI_TET_MESH, MESH_3D, in_name, in_info);
 
-	node_num = in_nodes.size();
+	node_num_ = in_nodes.size();
-	ele_num  = in_tets.size();
+	ele_num_  = in_tets.size();
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	for (int i = 0; i < node_num; i++)
+	for (int i = 0; i < node_num_; i++)
 	{
 		nodes[i] = in_nodes[i];
 	}
 
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
-	for (int i = 0; i < ele_num; i++)
+	for (int i = 0; i < ele_num_; i++)
 	{
 		elements[i].id = i;
 		for (int j = 0; j < 4; j++)
@@ -57,16 +53,20 @@ void gctl::tetrahedron_mesh::init(std::string in_name, std::string in_info, cons
 		elements[i].deter_vert_order();
 	}
 
-	initialized = true;
+	initialized_ = true;
+	return;
+}
+
+void gctl::tetrahedron_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::tetrahedron_mesh::load_binary(std::string filename)
 {
-	if (initialized)
+	check_initiated(true); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::tetrahedron_mesh] The mesh is already initialized.");
-	}
 
 	std::ifstream infile;
 	gctl::open_infile(infile, filename, ".2m", 
@@ -76,19 +76,19 @@ void gctl::tetrahedron_mesh::load_binary(std::string filename)
 	load_headinfo(infile, TRI_TET_MESH, MESH_2D);
 
 	// 读入网格信息
-	infile.read((char*)&node_num, sizeof(int));
+	infile.read((char*)&node_num_, sizeof(int));
-	infile.read((char*)&ele_num, sizeof(int));
+	infile.read((char*)&ele_num_, sizeof(int));
 
-	nodes.resize(node_num);
+	nodes.resize(node_num_);
-	elements.resize(ele_num);
+	elements.resize(ele_num_);
 
-	for (int i = 0; i < node_num; i++)
+	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++)
+	for (int i = 0; i < ele_num_; i++)
 	{
 		elements[i].id = i;
 		for (int j = 0; j < 4; j++)
@@ -98,21 +98,18 @@ void gctl::tetrahedron_mesh::load_binary(std::string filename)
 		}
 		elements[i].deter_vert_order();
 	}
-	initialized = true;
 
 	// 读入模型数据单元
 	load_datablock(infile);
 
 	infile.close();
+	initialized_ = true;
 	return;
 }
 
 void gctl::tetrahedron_mesh::save_binary(std::string filename)
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::tetrahedron_mesh] The mesh is not initialized.");
-	}
 
 	std::ofstream outfile;
 	gctl::open_outfile(outfile, filename, ".2m", 
@@ -122,16 +119,16 @@ void gctl::tetrahedron_mesh::save_binary(std::string filename)
 	save_headinfo(outfile);
 
 	// 输出网格信息
-	outfile.write((char*)&node_num, sizeof(int));
+	outfile.write((char*)&node_num_, sizeof(int));
-	outfile.write((char*)&ele_num, sizeof(int));
+	outfile.write((char*)&ele_num_, sizeof(int));
 
-	for (int i = 0; i < node_num; i++)
+	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 i = 0; i < ele_num_; i++)
 	{
 		for (int j = 0; j < 4; j++)
 		{
@@ -159,21 +156,13 @@ gctl::tetrahedron_mesh::~tetrahedron_mesh(){}
 
 const gctl::array<gctl::vertex3dc> &gctl::tetrahedron_mesh::get_nodes() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::tetrahedron_mesh] The mesh is not initialized.");
-	}
-
 	return nodes;
 }
 
 const gctl::array<gctl::tetrahedron> &gctl::tetrahedron_mesh::get_elements() const
 {
-	if (!initialized)
+	check_initiated(); // 检查网格是否已经初始化
-	{
-		throw std::runtime_error("[gctl::tetrahedron_mesh] The mesh is not initialized.");
-	}
-
 	return elements;
 }
 
@@ -183,11 +172,10 @@ void gctl::tetrahedron_mesh::load_tetgen(std::string filename, index_packed_e pa
 	gctl::read_Tetgen_element(filename, elements, nodes, packed);
 
 	// 设置名称与信息等
-	node_num = nodes.size();
+	node_num_ = nodes.size();
-	ele_num  = elements.size();
+	ele_num_  = elements.size();
 	base_mesh::init(TRI_TET_MESH, MESH_3D, filename, "Imported from a .node and .ele file.");
-	initialized = true;
+	initialized_ = true;
-
 	return;
 }
 
@@ -200,11 +188,10 @@ void gctl::tetrahedron_mesh::load_gmsh(std::string filename, index_packed_e pack
 	infile.close();
 
 	// 设置名称与信息等
-	node_num = nodes.size();
+	node_num_ = nodes.size();
-	ele_num  = elements.size();
+	ele_num_  = elements.size();
 	base_mesh::init(TRI_TET_MESH, MESH_3D, filename, "Imported from a .msh file.");
-	initialized = true;
+	initialized_ = true;
-
 	return;
 }
 
@@ -216,15 +203,3 @@ void gctl::tetrahedron_mesh::save_gmsh(std::string filename, index_packed_e pack
 	outfile.close();
 	return;
 }
-
-void gctl::tetrahedron_mesh::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
-	return;
-}
-
-void gctl::tetrahedron_mesh::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, datname, out_mode, packed);
-	return;
-}

lib/mesh/tet_mesh.h

@@ -40,12 +40,14 @@ namespace gctl
 		 * 以下为mesh类型的虚函数实现
 		 */
 
-		void init(std::string in_name, std::string in_info, const array<vertex3dc> &in_nodes, 
+		void load_gmsh(std::string filename, index_packed_e packed = Packed);
-			const array<tetrahedron> &in_tets);
+		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的专有函数
 		 */
@@ -55,14 +57,13 @@ namespace gctl
 			const array<tetrahedron> &in_tets);
 		virtual ~tetrahedron_mesh();
 
+		void init(std::string in_name, std::string in_info, const array<vertex3dc> &in_nodes, 
+			const array<tetrahedron> &in_tets);
+
 		const array<vertex3dc> &get_nodes() const;
 		const array<tetrahedron> &get_elements() const;
 
 		void load_tetgen(std::string filename, index_packed_e packed = Packed);
-		void load_gmsh(std::string filename, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
-		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
 
 	protected:
 		array<vertex3dc> nodes;

lib/mesh/tri2d_mesh.cpp

@@ -0,0 +1,283 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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 "tri2d_mesh.h"
+
+void gctl::triangle2d_mesh::init(std::string in_name, std::string in_info, const array<vertex2dc> &in_nodes_, 
+	const array<triangle2d> &in_triangles)
+{
+	check_initiated(true); // 检查是否已经初始化
+	base_mesh::init(TRI_TET_MESH, MESH_2D, 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::triangle2d_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::triangle2d_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_2D);
+
+	// 读入网格信息
+	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::vertex2dc));
+	}
+
+	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::triangle2d_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::vertex2dc));
+	}
+
+	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::triangle2d_mesh::triangle2d_mesh() : base_mesh::base_mesh(){}
+
+gctl::triangle2d_mesh::triangle2d_mesh(std::string in_name, std::string in_info, const array<vertex2dc> &in_nodes_, 
+	const array<triangle2d> &in_triangles)
+{
+	init(in_name, in_info, in_nodes_, in_triangles);
+}
+
+gctl::triangle2d_mesh::~triangle2d_mesh(){}
+
+const gctl::array<gctl::vertex2dc> &gctl::triangle2d_mesh::get_nodes() const
+{
+	check_initiated(); // 检查是否已经初始化
+	return nodes_;
+}
+
+const gctl::array<gctl::triangle2d> &gctl::triangle2d_mesh::get_elements() const
+{
+	check_initiated(); // 检查是否已经初始化
+	return elems_;
+}
+
+void gctl::triangle2d_mesh::reorder2anticlockwise()
+{
+	check_initiated(); // 检查是否已经初始化
+	vertex2dc *v_ptr;
+    for (size_t i = 0; i < get_elenum(); i++)
+    {
+        if (cross(*elems_[i].vert[1] - *elems_[i].vert[0], *elems_[i].vert[2] - *elems_[i].vert[0]) < 0)
+        {
+            v_ptr = elems_[i].vert[1]; 
+            elems_[i].vert[1] = elems_[i].vert[2];
+            elems_[i].vert[2] = v_ptr;
+        }
+    }
+    return;
+}
+
+void gctl::triangle2d_mesh::load_triangle(std::string filename, index_packed_e packed)
+{
+	gctl::read_Triangle_node(filename, nodes_, packed);
+	gctl::read_Triangle_element(filename, elems_, nodes_, packed);
+
+	// 设置名称与信息等
+	node_num_ = nodes_.size();
+	ele_num_  = elems_.size();
+	base_mesh::init(TRI_TET_MESH, MESH_2D, filename, "Imported from a .node and .ele file.");
+	initialized_ = true;
+	return;
+}
+
+void gctl::triangle2d_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_2D, filename, "Imported from a .msh file.");
+	initialized_ = true;
+	return;
+}
+
+void gctl::triangle2d_mesh::load_gmsh_groups()
+{
+	check_initiated(); // 检查是否已经初始化
+	meshio_.read_physical_groups(groups_);
+	return;
+}
+
+void gctl::triangle2d_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::triangle2d_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/tri2d_mesh.h

@@ -32,7 +32,7 @@
 
 namespace gctl
 {
-	class triangle_mesh : public base_mesh
+	class triangle2d_mesh : public base_mesh
 	{
 	public:
 
@@ -40,33 +40,48 @@ namespace gctl
 		 * 以下为mesh类型的虚函数实现
 		 */
 
-		void init(std::string in_name, std::string in_info, const array<vertex2dc> &in_nodes, 
+		void load_gmsh(std::string filename, index_packed_e packed = Packed);
-			const array<triangle2d> &in_triangles);
+		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的专有函数
+		 * 以下为triangle2d_mesh的专有函数
 		 */
 
-		triangle_mesh();
+		triangle2d_mesh();
-		triangle_mesh(std::string in_name, std::string in_info, const array<vertex2dc> &in_nodes, 
+		triangle2d_mesh(std::string in_name, std::string in_info, const array<vertex2dc> &in_nodes, 
+			const array<triangle2d> &in_triangles);
+		virtual ~triangle2d_mesh();
+
+		void init(std::string in_name, std::string in_info, const array<vertex2dc> &in_nodes, 
 			const array<triangle2d> &in_triangles);
-		virtual ~triangle_mesh();
 
 		const array<vertex2dc> &get_nodes() const;
 		const array<triangle2d> &get_elements() const;
 
+		void reorder2anticlockwise();
+
 		void load_triangle(std::string filename, index_packed_e packed = Packed);
-		void load_gmsh(std::string filename, index_packed_e packed = Packed);
+		void load_gmsh_groups();
-		void save_gmsh(std::string filename, index_packed_e packed = Packed);
+
-		void save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed = Packed);
+		/**
-		void save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed = Packed);
+		 * @brief 将模型物理分组转换为模型数据组	
+		 * 
+		 * @param datname 赋值给模型数据组的名称（如果数组不存在则新建）
+		 * @param phynames 物理组名称
+		 * @param phyvals 物理组值
+		 */
+		void groups2data(std::string datname, _1s_vector phynames, _1d_vector phyvals);
 
 	protected:
-		array<vertex2dc> nodes;
+		array<vertex2dc> nodes_;
-		array<triangle2d> elements;
+		array<triangle2d> elems_;
+		array<gmsh_physical_group> groups_;
+		_2i_vector elems_tag_;
 	};
 }
 

lib/mesh/tri_mesh.cpp

@@ -1,227 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <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<vertex2dc> &in_nodes, 
-	const array<triangle2d> &in_triangles)
-{
-	if (initialized)
-	{
-		throw std::runtime_error("[gctl::triangle_mesh] The mesh is already initialized.");
-	}
-
-	base_mesh::init(TRI_TET_MESH, MESH_2D, 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];
-	}
-
-	elements.resize(ele_num);
-	for (int i = 0; i < ele_num; i++)
-	{
-		elements[i].id = i;
-		for (int j = 0; j < 3; j++)
-		{
-			elements[i].vert[j] = nodes.get(in_triangles[i].vert[j]->id);
-		}
-	}
-	initialized = true;
-	return;
-}
-
-void gctl::triangle_mesh::load_binary(std::string filename)
-{
-	if (initialized)
-	{
-		throw std::runtime_error("[gctl::triangle_mesh] The mesh is already initialized.");
-	}
-
-	std::ifstream infile;
-	gctl::open_infile(infile, filename, ".2m", 
-		std::ios::in|std::ios::binary);
-
-	// 读入网格头信息
-	load_headinfo(infile, TRI_TET_MESH, MESH_2D);
-
-	// 读入网格信息
-	infile.read((char*)&node_num, sizeof(int));
-	infile.read((char*)&ele_num, sizeof(int));
-
-	nodes.resize(node_num);
-	elements.resize(ele_num);
-
-	for (int i = 0; i < node_num; i++)
-	{
-		infile.read((char*)nodes.get(i), sizeof(gctl::vertex2dc));
-	}
-
-	int in_index;
-	for (int i = 0; i < ele_num; i++)
-	{
-		elements[i].id = i;
-		for (int j = 0; j < 3; j++)
-		{
-			infile.read((char*)&in_index, sizeof(int));
-			elements[i].vert[j] = nodes.get(in_index);
-		}
-	}
-	initialized = true;
-
-	// 读入模型数据单元
-	load_datablock(infile);
-
-	infile.close();
-	return;
-}
-
-void gctl::triangle_mesh::save_binary(std::string filename)
-{
-	if (!initialized)
-	{
-		throw std::runtime_error("[gctl::triangle_mesh] The mesh is not initialized.");
-	}
-
-	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::vertex2dc));
-	}
-
-	int in_index;
-	for (int i = 0; i < ele_num; i++)
-	{
-		for (int j = 0; j < 3; j++)
-		{
-			in_index = elements[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<vertex2dc> &in_nodes, 
-	const array<triangle2d> &in_triangles)
-{
-	init(in_name, in_info, in_nodes, in_triangles);
-}
-
-gctl::triangle_mesh::~triangle_mesh(){}
-
-const gctl::array<gctl::vertex2dc> &gctl::triangle_mesh::get_nodes() const
-{
-	if (!initialized)
-	{
-		throw std::runtime_error("[gctl::triangle_mesh] The mesh is not initialized.");
-	}
-
-	return nodes;
-}
-
-const gctl::array<gctl::triangle2d> &gctl::triangle_mesh::get_elements() const
-{
-	if (!initialized)
-	{
-		throw std::runtime_error("[gctl::triangle_mesh] The mesh is not initialized.");
-	}
-
-	return elements;
-}
-
-void gctl::triangle_mesh::load_triangle(std::string filename, index_packed_e packed)
-{
-	gctl::read_Triangle_node(filename, nodes, packed);
-	gctl::read_Triangle_element(filename, elements, nodes, packed);
-
-	// 设置名称与信息等
-	node_num = nodes.size();
-	ele_num  = elements.size();
-	base_mesh::init(TRI_TET_MESH, MESH_2D, filename, "Imported from a .node and .ele file.");
-	initialized = true;
-
-	return;
-}
-
-void gctl::triangle_mesh::load_gmsh(std::string filename, index_packed_e packed)
-{
-	std::ifstream infile;
-	gctl::open_infile(infile, filename, ".msh");
-	gctl::read_gmsh_node(infile, nodes, packed);
-	gctl::read_gmsh_element(infile, elements, nodes, packed);
-	infile.close();
-
-	// 设置名称与信息等
-	node_num = nodes.size();
-	ele_num  = elements.size();
-	base_mesh::init(TRI_TET_MESH, MESH_2D, filename, "Imported from a .msh file.");
-	initialized = true;
-
-	return;
-}
-
-void gctl::triangle_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();
-	return;
-}
-
-void gctl::triangle_mesh::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
-	return;
-}
-
-void gctl::triangle_mesh::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
-{
-	base_mesh::save_gmsh(filename, datname, out_mode, packed);
-	return;
-}

tool/gridmanager/gridmanager.cpp

@@ -377,8 +377,7 @@ void save_gmsh(const std::vector<std::string> &cmd_units)
     // save gmsh <file>
     if (cmd_units.size() < 3) throw std::runtime_error("save: insufficient parameters.");
 
-    rg.save_gmsh(cmd_units[2], NodeData, OverWrite, NotPacked);
+    rg.save_gmsh_withdata(cmd_units[2], OverWrite, NotPacked);
-    rg.save_gmsh(cmd_units[2], ElemData, Append, NotPacked);
     return;
 }
 
@@ -477,7 +476,7 @@ void data_cloud(const std::vector<std::string> &cmd_units)
     array<point2dc> posi_arr(posix_vec.size());
     array<double> posi_val;
 
-    posi_val.import_vector(data_vec);
+    posi_val.input(data_vec);
     for (size_t i = 0; i < posi_arr.size(); i++)
     {
         posi_arr[i].x = posix_vec[i];
@@ -622,21 +621,21 @@ void data_output(const std::vector<std::string> &cmd_units)
         copy_str[i - 1] = cmd_units[i];
     }
 
-    meshdata* data_ptr;
+    meshdata *data_ptr;
     if (copy_str[0] == "enable")
     {
         for (size_t i = 1; i < copy_str.size(); i++)
         {
-            data_ptr = rg.get_data(copy_str[i]);
+            data_ptr = rg.get_data_ptr(copy_str[i]);
-            data_ptr->set_output(true);
+            data_ptr->output_ok_ = true;
         }
     }
     else if (copy_str[0] == "disable")
     {
         for (size_t i = 1; i < copy_str.size(); i++)
         {
-            data_ptr = rg.get_data(copy_str[i]);
+            data_ptr = rg.get_data_ptr(copy_str[i]);
-            data_ptr->set_output(false);
+            data_ptr->output_ok_ = false;
         }
     }
     else throw std::runtime_error("data-output: invalid operation type.");
@@ -654,7 +653,8 @@ void data_rename(const std::vector<std::string> &cmd_units)
         copy_str[i - 1] = cmd_units[i];
     }
 
-    rg.rename_data(copy_str[0], copy_str[1]);
+    meshdata &curr_data = rg.get_data(copy_str[0]);
+    curr_data.name_ = copy_str[1];
     return;
 }
 
@@ -697,7 +697,7 @@ void get_stats(const std::vector<std::string> &cmd_units)
     std::vector<double> stats;
     for (size_t i = 0; i < copy_str.size(); i++)
     {
-        data_ptr = rg.get_data(copy_str[i]);
+        data_ptr = rg.get_data_ptr(copy_str[i]);
         data_ptr->show_info();
         data_ptr->show_stats();
     }