diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 6b39b6d..fd48dbb 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -23,7 +23,7 @@ add_example(windowfunc_ex OFF) add_example(legendre_ex OFF) add_example(refellipsoid_ex OFF) add_example(kde_ex OFF) -add_example(meshio_ex OFF) +add_example(meshio_ex ON) add_example(autodiff_ex OFF) add_example(multinary_ex OFF) add_example(text_io_ex OFF) diff --git a/example/meshio_ex.cpp b/example/meshio_ex.cpp index 1b74f20..0bb3399 100644 --- a/example/meshio_ex.cpp +++ b/example/meshio_ex.cpp @@ -33,28 +33,35 @@ using namespace gctl; int main(int argc, char const *argv[]) try { mesh_io mshio; - - //mshio.read_tetgen_ascii("tmp/ex1.1"); - //mshio.edit_group(Disable, GeometryTag, 5); - //mshio.edit_group(GeometryTag, 1, PhysicalTag, 1); - //mshio.edit_group(GeometryTag, 2, PhysicalTag, 2); - //mshio.edit_group(GeometryTag, 3, PhysicalTag, 3); - //mshio.edit_group(GeometryTag, 4, PhysicalTag, 4); - //mshio.edit_group(GeometryTag, 1, "Boundary"); - //mshio.edit_group(GeometryTag, 2, "Body1"); - //mshio.edit_group(GeometryTag, 3, "Body2"); - //mshio.edit_group(GeometryTag, 4, "Body3"); - //mshio.save_gmsh_v2_ascii("tmp/ex1.1"); -/* + +/* + mshio.read_tetgen_ascii("tmp/ex1.1"); + mshio.edit_group(Disable, GeometryTag, 5); + mshio.edit_group(GeometryTag, 1, PhysicalTag, 1); + mshio.edit_group(GeometryTag, 2, PhysicalTag, 2); + mshio.edit_group(GeometryTag, 3, PhysicalTag, 3); + mshio.edit_group(GeometryTag, 4, PhysicalTag, 4); + mshio.edit_group(GeometryTag, 1, "Boundary"); + mshio.edit_group(GeometryTag, 2, "Body1"); + mshio.edit_group(GeometryTag, 3, "Body2"); + mshio.edit_group(GeometryTag, 4, "Body3"); + mshio.save_gmsh_v2_ascii("tmp/ex1.1"); +*/ + mshio.read_gmsh_v2_ascii("tmp/ex1.1"); mshio.convert_tags_to_data(GeometryTag); + array body_val(mshio.element_size("Body2"), 2.0); - mshio.create_data(body_val, "BodyValue", "Body2"); + mshio.add_element_data(body_val, "BodyValue", "Body2"); + + array body_val2(mshio.element_size("Body3"), 1.0); + mshio.add_element_data(body_val2, "BodyValue", "Body3"); + mshio.save_gmsh_v2_ascii("tmp/ex1.2"); //mshio.save_vtk_legacy_ascii("tmp/ex1.1"); mshio.info(); - array nodes = mshio.get_nodes(); + const array &nodes = mshio.get_nodes(); array body2_tets; mshio.export_elements_to(body2_tets, "All"); @@ -63,8 +70,7 @@ int main(int argc, char const *argv[]) try gio.init_file("tmp.msh", Output); gio.set_packed(NotPacked, Output); gio.save_mesh(body2_tets, nodes); -*/ - +/* mshio.read_gmsh_v2_ascii("tmp/wjb.1"); mshio.edit_group(Disable); mshio.edit_group(Enable, GeometryTag, 3); @@ -72,6 +78,7 @@ int main(int argc, char const *argv[]) try mshio.edit_group(Enable, GeometryTag, 9); mshio.save_gmsh_v2_ascii("tmp/wjb.2"); +*/ return 0; } catch(std::exception &e) diff --git a/lib/io/mesh_io.cpp b/lib/io/mesh_io.cpp index 42bdd7b..1fa8338 100644 --- a/lib/io/mesh_io.cpp +++ b/lib/io/mesh_io.cpp @@ -564,13 +564,12 @@ void gctl::mesh_io::convert_tags_to_data(element_tag_enum tag_type) } else { - int t = 0; tmp_data.enabled = true; tmp_data.d_type = ElemData; - if (tag_type == PhysicalTag) {tmp_data.str_tag.resize(1, "Physical Tag"); t = 0;} - else if (tag_type == GeometryTag) {tmp_data.str_tag.resize(1, "Geometry Tag"); t = 1;} - else if (tag_type == PartitionTag) {tmp_data.str_tag.resize(1, "Partition Tag"); t = 2;} + if (tag_type == PhysicalTag) tmp_data.str_tag.resize(1, "Physical Tag"); + else if (tag_type == GeometryTag) tmp_data.str_tag.resize(1, "Geometry Tag"); + else if (tag_type == PartitionTag) tmp_data.str_tag.resize(1, "Partition Tag"); tmp_data.real_tag.resize(1, 0.0); tmp_data.int_tag.resize(3, 0); @@ -686,7 +685,7 @@ void gctl::mesh_io::export_elements_to(array &tets, std::string phy tets.resize(s); s = 0; - for (size_t i = 0; i < elems_.size(); i++) + for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && (groups_[i].name == phys_name || phys_name == "All")) @@ -722,7 +721,7 @@ void gctl::mesh_io::export_elements_to(array &tets, element_tag_enu tets.resize(s); s = 0; - for (size_t i = 0; i < elems_.size(); i++) + for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].type == _4NodeTetrahedron && (tag_type == PhysicalTag && groups_[i].phys_group == tag) || @@ -743,35 +742,101 @@ void gctl::mesh_io::export_elements_to(array &tets, element_tag_enu return; } -void gctl::mesh_io::create_node_data(const array &data, std::string name) +int gctl::mesh_io::if_saved_data(std::string name) { + for (size_t i = 0; i < datas_.size(); i++) + { + if (datas_[i].str_tag.front() == name) return i; + } + return -1; +} + +void gctl::mesh_io::add_node_data(const array &data, std::string name) +{ + size_t s = nodes_.size(); + if (data.size()!= s) throw std::runtime_error("[gctl::mesh_io::create_node_data] Incompatible data size."); + + int d_id = if_saved_data(name); + if (d_id != -1) + { + for (size_t i = 0; i < s; i++) + { + datas_[d_id].val[i] = data[i]; + } + return; + } + mesh_data new_data; new_data.enabled = true; new_data.d_type = NodeData; new_data.str_tag.resize(1, name); + new_data.real_tag.resize(1, 0.0); new_data.int_tag.resize(3, 0); new_data.int_tag[1] = 1; + new_data.int_tag[2] = s; + new_data.vert_ptrs.resize(s, nullptr); + new_data.val.resize(s); - size_t c = 0; - for (size_t i = 0; i < nodes_.size(); i++) + for (size_t i = 0; i < s; i++) { - if (nodes_[i].id != DEFAULT_INVALID_TAG) c++; + new_data.vert_ptrs[i] = nodes_.get(i); + new_data.val[i] = data[i]; } - if (data.size() != c) throw std::runtime_error("[gctl::mesh_io::create_node_data] Incompatible data size."); - - new_data.int_tag[2] = c; - new_data.vert_ptrs.resize(c, nullptr); - new_data.val.resize(c); + datas_.push_back(new_data); + return; +} - c = 0; +void gctl::mesh_io::add_node_data(const array &data, const array &boolen, std::string name) +{ + size_t s = nodes_.size(); + if (data.size()!= s || boolen.size() != s) throw std::runtime_error("[gctl::mesh_io::create_node_data] Incompatible data size."); + + s = 0; for (size_t i = 0; i < nodes_.size(); i++) { - if (nodes_[i].id != DEFAULT_INVALID_TAG) + if (boolen[i]) s++; + } + + int d_id = if_saved_data(name); + if (d_id != -1) + { + datas_[d_id].val.resize(s); + datas_[d_id].vert_ptrs.resize(s, nullptr); + datas_[d_id].int_tag[2] = s; + + s = 0; + for (size_t i = 0; i < nodes_.size(); i++) { - new_data.vert_ptrs[c] = nodes_.get(i); - new_data.val[c] = data[c]; - c++; + if (boolen[i]) + { + datas_[d_id].vert_ptrs[s] = nodes_.get(i); + datas_[d_id].val[s] = data[i]; + s++; + } + } + return; + } + + mesh_data new_data; + new_data.enabled = true; + new_data.d_type = NodeData; + new_data.str_tag.resize(1, name); + new_data.real_tag.resize(1, 0.0); + new_data.int_tag.resize(3, 0); + new_data.int_tag[1] = 1; + new_data.int_tag[2] = s; + new_data.vert_ptrs.resize(s, nullptr); + new_data.val.resize(s); + + s = 0; + for (size_t i = 0; i < nodes_.size(); i++) + { + if (boolen[i]) + { + new_data.vert_ptrs[s] = nodes_.get(i); + new_data.val[s] = data[i]; + s++; } } @@ -779,12 +844,13 @@ void gctl::mesh_io::create_node_data(const array &data, std::string name return; } -void gctl::mesh_io::create_element_data(const array &data, std::string name, element_type_enum e_type) +void gctl::mesh_io::add_element_data(const array &data, std::string name, element_type_enum e_type) { mesh_data new_data; new_data.enabled = true; new_data.d_type = ElemData; new_data.str_tag.resize(1, name); + new_data.real_tag.resize(1, 0.0); new_data.int_tag.resize(3, 0); new_data.int_tag[1] = 1; @@ -818,12 +884,13 @@ void gctl::mesh_io::create_element_data(const array &data, std::string n return; } -void gctl::mesh_io::create_element_data(const array &data, std::string name, element_tag_enum tag_type, int tag) +void gctl::mesh_io::add_element_data(const array &data, std::string name, element_tag_enum tag_type, int tag) { mesh_data new_data; new_data.enabled = true; new_data.d_type = ElemData; new_data.str_tag.resize(1, name); + new_data.real_tag.resize(1, 0.0); new_data.int_tag.resize(3, 0); new_data.int_tag[1] = 1; @@ -866,15 +933,8 @@ void gctl::mesh_io::create_element_data(const array &data, std::string n return; } -void gctl::mesh_io::create_element_data(const array &data, std::string name, std::string phys_name) +void gctl::mesh_io::add_element_data(const array &data, std::string name, std::string phys_name) { - mesh_data new_data; - new_data.enabled = true; - new_data.d_type = ElemData; - new_data.str_tag.resize(1, name); - new_data.int_tag.resize(3, 0); - new_data.int_tag[1] = 1; - size_t e = 0; for (size_t i = 0; i < groups_.size(); i++) { @@ -886,11 +946,44 @@ void gctl::mesh_io::create_element_data(const array &data, std::string n if (data.size() != e) throw std::runtime_error("[gctl::mesh_io::create_element_data] Incompatible data size."); + int d_id = if_saved_data(name); + if (d_id != -1) + { + array more_elem_ptrs(e, nullptr); + array more_val(e, 0.0); + + e = 0; + for (size_t i = 0; i < groups_.size(); i++) + { + if (groups_[i].enabled && groups_[i].name == phys_name) + { + for (size_t j = 0; j < groups_[i].elem_ptrs.size(); j++) + { + more_elem_ptrs[e] = groups_[i].elem_ptrs[j]; + more_val[e] = data[e]; + e++; + } + } + } + + datas_[d_id].int_tag[2] += e; + datas_[d_id].elem_ptrs.concat(more_elem_ptrs); + datas_[d_id].val.concat(more_val); + return; + } + + mesh_data new_data; + new_data.enabled = true; + new_data.d_type = ElemData; + new_data.str_tag.resize(1, name); + new_data.real_tag.resize(1, 0.0); + new_data.int_tag.resize(3, 0); + new_data.int_tag[1] = 1; new_data.int_tag[2] = e; new_data.elem_ptrs.resize(e, nullptr); new_data.val.resize(e, 0.0); - e = 0; + e = 0; for (size_t i = 0; i < groups_.size(); i++) { if (groups_[i].enabled && groups_[i].name == phys_name) @@ -1205,7 +1298,7 @@ void gctl::mesh_io::read_gmsh_v2_ascii(std::string filename, index_packed_e is_p size_t p_size = 0; size_t n_size = 0; array phys; - + while(getline(infile,tmp_str)) { if (tmp_str == "$PhysicalNames") @@ -1268,28 +1361,12 @@ void gctl::mesh_io::read_gmsh_v2_ascii(std::string filename, index_packed_e is_p elems_[i].type = elem_gmsh_type(type_code); elems_[i].vert_ptrs.resize(elem_size(elems_[i].type)); - if (attri_num >= 3) // default tags will be assgined to DEFAULT_INVALID_TAG + // we get at least three tags. see below for tag types + // default tags will be assgined to DEFAULT_INVALID_TAG + int_tag[i].resize(GCTL_MAX(3, attri_num), DEFAULT_INVALID_TAG); + for (size_t a = 0; a < attri_num; a++) { - int_tag[i].resize(attri_num); - for (size_t a = 0; a < attri_num; a++) - { - tmp_ss >> int_tag[i][a]; - } - } - else if (attri_num == 2) - { - int_tag[i].resize(3, DEFAULT_INVALID_TAG); - tmp_ss >> int_tag[i][0]; - tmp_ss >> int_tag[i][1]; - } - else if (attri_num == 1) - { - int_tag[i].resize(3, DEFAULT_INVALID_TAG); - tmp_ss >> int_tag[i][0]; - } - else - { - int_tag[i].resize(3, DEFAULT_INVALID_TAG); + tmp_ss >> int_tag[i][a]; } for (size_t v = 0; v < elems_[i].vert_ptrs.size(); v++) @@ -1392,6 +1469,7 @@ void gctl::mesh_io::read_gmsh_v2_ascii(std::string filename, index_packed_e is_p tmp_group.elem_ptrs.push_back(elems_.get(0)); groups_.push_back(tmp_group); + // 元素类型与三个标记都一致的元素将被分到同一组 bool not_found; for (size_t i = 1; i < elems_.size(); i++) { @@ -1413,9 +1491,9 @@ void gctl::mesh_io::read_gmsh_v2_ascii(std::string filename, index_packed_e is_p { tmp_group.elem_ptrs.clear(); tmp_group.type = elems_[i].type; - tmp_group.phys_group = int_tag[i][0]; - tmp_group.geom_group = int_tag[i][1]; - tmp_group.part_group = int_tag[i][2]; + tmp_group.phys_group = int_tag[i][0]; // 物理组 + tmp_group.geom_group = int_tag[i][1]; // 几何组 + tmp_group.part_group = int_tag[i][2]; // 剖分组(一般以元素的维度区分) tmp_group.elem_ptrs.push_back(elems_.get(i)); groups_.push_back(tmp_group); } @@ -1423,6 +1501,7 @@ void gctl::mesh_io::read_gmsh_v2_ascii(std::string filename, index_packed_e is_p if (!phys.empty()) { + // 遍历所有元素组 按物理组为标准为元素组命名 并以将元素维度赋值为剖分组 for (size_t g = 0; g < groups_.size(); g++) { for (size_t p = 0; p < phys.size(); p++) diff --git a/lib/io/mesh_io.h b/lib/io/mesh_io.h index d7b63f1..7a3ca62 100644 --- a/lib/io/mesh_io.h +++ b/lib/io/mesh_io.h @@ -93,10 +93,10 @@ namespace gctl */ enum element_tag_enum { - PhysicalTag, - GeometryTag, - PartitionTag, - NodeTag, + PhysicalTag, // 元素的物理分组标签 + GeometryTag, // 元素的几何分组标签 + PartitionTag, // 元素的剖分分组标签 + NodeTag, // 顶点的标签(仅用于输出顶点标签数据) }; /** @@ -105,11 +105,11 @@ namespace gctl */ struct mesh_element { - bool enabled; - int id; - element_type_enum type; - array vert_ptrs; - array neigh_ptrs; + bool enabled; // 单元体是否有效 + int id; // 单元体编号 + element_type_enum type; // 单元体类型 + array vert_ptrs; // 顶点指针数组 + array neigh_ptrs; // 相邻单元体指针数组 mesh_element(); }; @@ -120,14 +120,14 @@ namespace gctl */ struct mesh_data { - bool enabled; - mesh_data_type_e d_type; - array str_tag; - array real_tag; - array int_tag; + bool enabled; // 数据体是否有效 + mesh_data_type_e d_type; // 数据类型 + array str_tag; // 字符串类型的标签(默认为一个,即为数据的名称) + array real_tag; // 实数类型的标签(默认为一个,等于0.0) + array int_tag; // 整数类型的标签(最少三个,最后一个为数据的长度) array vert_ptrs; // 两者只能存在一个 array elem_ptrs; // 两者只能存在一个 - array val; + array val; // 数据值 mesh_data(); @@ -150,13 +150,13 @@ namespace gctl */ struct mesh_element_group { - bool enabled; - element_type_enum type; - int phys_group; - int geom_group; - int part_group; - std::string name; - std::vector elem_ptrs; + bool enabled; // 组是否有效 + element_type_enum type; // 组内单元体类型 + int phys_group; // 物理分组标签 + int geom_group; // 几何分组标签 + int part_group; // 剖分分组标签 + std::string name; // 组名 + std::vector elem_ptrs; // 组内单元体指针数组 mesh_element_group(); @@ -256,14 +256,14 @@ namespace gctl const array &get_node_tag(); /** - * @brief 返回顶点数组的引用。 + * @brief 返回所有顶点数组的引用。 * * @return 顶点数组的引用。 */ const array &get_nodes(); /** - * @brief 返回单元体数组的引用。 + * @brief 返回所有单元体数组的引用。 * * @return 单元体数组的引用。 */ @@ -336,40 +336,68 @@ namespace gctl void export_elements_to(array &tets, element_tag_enum tag_type, int tag); /** - * @brief 创建一个顶点数据对象。 + * @brief 检查是否存在名为name的数据 * - * @param data 输入的数据数组。 - * @param name 新建的数据名称。 + * @param name 数据名称 + * + * @return 存在则返回数据索引,不存在则返回-1。 */ - void create_node_data(const array &data, std::string name); + int if_saved_data(std::string name); /** - * @brief 按单元体类型筛选创建一个单元体数据对象。 + * @brief 添加一个顶点数据对象。数据将依次添加到所有顶点位置。 + * + * @note 若对应名称的数据已经存在则会覆盖 + * + * @param data 输入的数据数组,长度与网格所有顶点数据相同。 + * @param name 新建的数据名称。 + */ + void add_node_data(const array &data, std::string name); + + /** + * @brief 添加一个顶点数据对象。数据将依次添加到布尔为真的顶点位置。 + * + * @note 若对应名称的数据已经存在则会覆盖 + * + * @param data 输入的数据数组,长度与网格所有顶点数据相同。 + * @param boolen 输入的布尔,只有为真元素位置的顶点数据将被保存。 + * @param name 新建的数据名称。 + */ + void add_node_data(const array &data, const array &boolen, std::string name); + + /** + * @brief 按单元体类型筛选创建一个单元体数据对象。数据将依次添加到所选元素位置。 + * + * @note 若对应名称的数据已经存在则会覆盖 * * @param data 输入的数据数组。 * @param name 新建数据名称。 * @param e_type 新建数据的单元体类型(缺省值为NotSet,表示选择所有有效的单元体)。 */ - void create_element_data(const array &data, std::string name, element_type_enum e_type = NotSet); + void add_element_data(const array &data, std::string name, element_type_enum e_type = NotSet); /** * @brief 按单元体标签值筛选创建一个单元体数据对象。 * + * @note 若对应名称的数据已经存在则会追加 + * * @param data 输入的数据数组。 * @param name 新建数据名称。 * @param tag_type 标签类型。 * @param tag 标签值。 */ - void create_element_data(const array &data, std::string name, element_tag_enum tag_type, int tag); + void add_element_data(const array &data, std::string name, element_tag_enum tag_type, int tag); /** * @brief 按单元体组的名称筛选创建一个单元体数据对象。 * + * @note 若对应名称的数据已经存在则会追加 + * * @param data 输入的数据数组。 * @param name 新建数据名称。 * @param phys_name 单元体组的名称。 */ - void create_element_data(const array &data, std::string name, std::string phys_name); + void add_element_data(const array &data, std::string name, std::string phys_name); /** * @brief 读入triangle软件输出的网格剖分文件。 @@ -422,30 +450,31 @@ namespace gctl void save_data_to_xyz(std::string filename, std::string dataname = "null", coordinate_system_e out_coor = Cartesian, double refr = GCTL_Earth_Radius, double refR = GCTL_Earth_Radius); private: - bool initialized_; - - size_t valid_node_size_, valid_elem_size_, valid_group_size_; - array nodes_; - array elems_; - std::vector datas_; - std::vector groups_; - array nodes_tag_; + bool initialized_; // 类型是否已经初始化完成 - element_type_enum elem_gmshcode2type_[94]; - element_type_enum elem_vtkcode2type_[14]; - std::map elem_type2gmshcode_; - std::map elem_type2vtkcode_; - std::map elem_type2size_; - std::map elem_type2name_; + // 有效的顶点、单元体和单元体组的数量 + size_t valid_node_size_, valid_elem_size_, valid_group_size_; + array nodes_; // 网格顶点 当顶点索引为无效值时将不会被输出 + array elems_; // 网格元素 + std::vector datas_; // 网格数据 + std::vector groups_; // 网格单元体组 + array nodes_tag_; // 顶点标签 + + element_type_enum elem_gmshcode2type_[94]; // gmsh的单元体类型数组 数组索引为gmsh的单元体类型码值 + element_type_enum elem_vtkcode2type_[14]; // vtk的单元体类型数组 数组索引为vtk的单元体类型码值 + std::map elem_type2gmshcode_; // 单元体类型到gmsh类型码值的映射 + std::map elem_type2vtkcode_; // 单元体类型到vtk类型码值的映射 + std::map elem_type2size_; // 单元体类型到单元体顶点数量的映射 + std::map elem_type2name_; // 单元体类型到单元体名称的映射 - std::string elem_name(element_type_enum e_type); - int elem_gmsh_code(element_type_enum e_type); - int elem_vtk_code(element_type_enum e_type); - int elem_size(element_type_enum e_type); - element_type_enum elem_gmsh_type(int code); - element_type_enum elem_vtk_type(int code); - void update_indexing(); - void sort_groups(); + std::string elem_name(element_type_enum e_type); // 获取单元体名称字符串 + int elem_gmsh_code(element_type_enum e_type); // 获取单元体gmsh类型码值 + int elem_vtk_code(element_type_enum e_type); // 获取单元体vtk类型码值 + int elem_size(element_type_enum e_type); // 获取单元体顶点数量 + element_type_enum elem_gmsh_type(int code); // 获取对应gmsh类型码的单元体类型 + element_type_enum elem_vtk_type(int code); // 获取对应vtk类型码的单元体类型 + void update_indexing(); // 更新索引(对网格的元素进行操作后需调用) + void sort_groups(); // 对单元体组进行梳理(对网格的元素进行操作后需调用) }; }