tmp

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)

example/meshio_ex.cpp

@@ -33,23 +33,25 @@ using namespace gctl;
 int main(int argc, char const *argv[]) try
 {
     mesh_io mshio;
-
-    mshio.read_gmsh_v2_ascii("tmp/example2");
+    mshio.read_gmsh_v2_ascii("tmp/test");
     mshio.info();
-/*
-    array<vertex2dc> nodes;
-    array<triangle2d> tris;
-    //array<vertex3dc*> orinode_ptrs;
-    //array<meshio_element*> oriele_ptrs;
 
-    mshio.select_elements("Survey");
+    array<vertex3dc> nodes;
+    array<triangle> tris;
+
+    mshio.select_elements("Bz (pT)", NodeData);
     mshio.export_to(tris, nodes);
 
-    gmshio tmp_out;
-    tmp_out.init_file("tmp", Output);
-    tmp_out.set_packed(NotPacked, Output);
-    tmp_out.save_mesh(tris, nodes);
-*/
+    meshio_data &d = mshio.get_data("Bz (pT)", NodeData);
+    array<double> dat = d.val;
+
+    mesh_io mshio2;
+    mshio2.import_from(tris, nodes);
+    mshio2.select_elements();
+    mshio2.add_data("Bz (pT)", dat, NodeData, OverWrite);
+    mshio2.save_gmsh_v2_ascii("tmp");
+    mshio2.info();
+
 /*  
     mshio.read_tetgen_ascii("tmp/ex1.1");
     mshio.edit_group(Disable, GeometryTag, 5);

lib/io/mesh_io.cpp

@@ -567,6 +567,78 @@ void gctl::mesh_io::select_elements(std::string phys_name)
     return;
 }
 
+void gctl::mesh_io::select_elements(std::string dat_name, mesh_data_type_e dtype)
+{
+    if (!selected_nodes_.empty()) selected_nodes_.clear();
+    if (!selected_elems_.empty()) selected_elems_.clear();
+    if (!selected_groups_.empty()) selected_groups_.clear();
+
+    int d_id = if_saved_data(dat_name, dtype);
+    if (d_id == DEFAULT_INVALID_TAG) throw std::runtime_error("[gctl::mesh_io::select_elements] Data not found.");
+
+    if (dtype == ElemData)
+    {
+        int vnum;
+        meshio_element *e_ptr;
+        for (size_t i = 0; i < datas_[d_id].tar_ptrs.size(); i++)
+        {
+            e_ptr = reinterpret_cast<meshio_element*>(datas_[d_id].tar_ptrs[i]);
+            if (e_ptr->enabled)
+            {
+                selected_elems_.push_back(e_ptr);
+                vnum = elem_size(e_ptr->type);
+                for (size_t v = 0; v < vnum; v++)
+                {
+                    selected_nodes_.push_back(e_ptr->vert_ptrs[v]);
+                }
+            }
+        }
+
+        std::vector<vertex3dc*>::iterator pos;
+        std::sort(selected_nodes_.begin(), selected_nodes_.end()); //排序
+        pos = std::unique(selected_nodes_.begin(), selected_nodes_.end()); //获取重复序列开始的位置
+        selected_nodes_.erase(pos, selected_nodes_.end()); //删除重复点
+    }
+    else // NodeData
+    {
+        vertex3dc *v_ptr;
+        std::unordered_set<vertex3dc*> nodes_set;
+        for (size_t i = 0; i < datas_[d_id].tar_ptrs.size(); i++)
+        {
+            v_ptr = reinterpret_cast<vertex3dc*>(datas_[d_id].tar_ptrs[i]);
+            if (v_ptr->id != DEFAULT_INVALID_TAG)
+            {
+                selected_nodes_.push_back(v_ptr);
+                nodes_set.insert(v_ptr);
+            }
+        }
+
+        int vnum;
+        bool skip_ele;
+        for (size_t i = 0; i < elems_.size(); i++)
+        {
+            if (elems_[i].enabled)
+            {
+                skip_ele = false;
+                vnum = elem_size(elems_[i].type);
+                for (size_t v = 0; v < vnum; v++)
+                {
+                    if (nodes_set.count(elems_[i].vert_ptrs[v]) <= 0)
+                    {
+                        skip_ele = true;
+                        break;
+                    }
+                }
+
+                if (skip_ele) break;
+                else selected_elems_.push_back(elems_.get(i));
+            }
+        }
+        
+    }
+    return;
+}
+
 size_t gctl::mesh_io::node_size()
 {
     return selected_nodes_.size();
@@ -579,10 +651,14 @@ size_t gctl::mesh_io::element_size()
 
 const gctl::array<int> &gctl::mesh_io::node_tags()
 {
-    selected_node_tag_.resize(selected_nodes_.size());
-    for (size_t i = 0; i < selected_node_tag_.size(); i++)
+    selected_node_tag_.resize(selected_nodes_.size(), DEFAULT_INVALID_TAG);
+
+    if (!nodes_tag_.empty())
     {
-        selected_node_tag_[i] = nodes_tag_[selected_nodes_[i]->id];
+        for (size_t i = 0; i < selected_node_tag_.size(); i++)
+        {
+            selected_node_tag_[i] = nodes_tag_[selected_nodes_[i]->id];
+        }
     }
     return selected_node_tag_;
 }

lib/io/mesh_io.h

@@ -28,6 +28,8 @@
 #ifndef _GCTL_MESH_IO_H
 #define _GCTL_MESH_IO_H
 
+#include <unordered_set>
+
 // library's head files
 #include "../core.h"
 #include "../geometry.h"
@@ -283,6 +285,14 @@ namespace gctl
          */
         void select_elements(std::string phys_name);
 
+        /**
+         * @brief 选择所有符合条件的单元体（注意只会统计有效的单元体组）。
+         * 
+         * @param datname 数据名称
+         * @param dtype 数据类型
+         */
+        void select_elements(std::string dat_name, mesh_data_type_e dtype);
+
         /**
          * @brief 返回已选择的单元体所对应的顶点数量（使用select_elements函数选择）。
          * 
@@ -298,14 +308,14 @@ namespace gctl
         size_t element_size();
 
         /**
-         * @brief 返回已选择的单元体所对应的顶点标签
+         * @brief 返回已选择的单元体所对应的顶点标签（默认的无效标签为-9999）
          * 
          * @return 整型数组的引用。
          */
         const array<int> &node_tags();
 
         /**
-         * @brief 返回已选择的单元体所对应的元素标签
+         * @brief 返回已选择的单元体所对应的元素标签（默认的无效标签为-9999）
          * 
          * @return 整型数组的引用。
          */
@@ -319,7 +329,7 @@ namespace gctl
         void get_gmsh_physical_groups(std::vector<gmsh_physical_group> &g_groups);
 
         /**
-         * @brief 将对应标签类型转换为网格数据类型
+         * @brief 将对应标签类型转换为网格数据类型（注意只会转换有效的顶点与单元体的标签）
          * 
          * @param tag_type 标签类型。
          */
@@ -336,7 +346,7 @@ namespace gctl
         int if_saved_data(std::string name, mesh_data_type_e type);
 
         /**
-         * @brief 获取数据对象的引用
+         * @brief 获取数据对象的引用（注意此函数会直接返回数据的引用，注意可能会包含无效的顶点和单元体）
          * 
          * @param name 数据名称
          * @param type 数据类型
@@ -423,7 +433,7 @@ 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);
 
         /**
-         * @brief 提取对应类型和单元体到外部数组。
+         * @brief 提取所选择的单元体到外部数组。
          * 
          * @tparam T 单元体类型。
          * @param elems 新生成的单元体数组的引用。
@@ -433,7 +443,7 @@ namespace gctl
         void export_to(array<T> &elems, array<vertex3dc> &nodes);
 
         /**
-         * @brief 提取对应类型和单元体到外部数组。
+         * @brief 提取所选择的单元体到外部数组。
          * 
          * @tparam T 单元体类型。
          * @param elems 新生成的单元体数组的引用。
@@ -447,6 +457,31 @@ namespace gctl
         template <typename T>
         void export_to(array<T> &elems, array<vertex2dc> &nodes, int x_id = 0, int y_id = 1);
 
+        /**
+         * @brief 导入外部单元体数组。
+         * 
+         * @tparam T 单元体类型。
+         * @param elems 单元体数组。
+         * @param nodes 顶点数组。
+         */
+        template <typename T>
+        void import_from(const array<T> &elems, const array<vertex3dc> &nodes);
+        
+        /**
+         * @brief 导入外部单元体数组。
+         * 
+         * @tparam T 单元体类型。
+         * @param elems 单元体数组。
+         * @param nodes 顶点数组。
+         * @param x_id 二维x坐标到三维坐标系统的索引（缺省值为0）。
+         * @param y_id 三维y坐标到三维坐标系统的索引（缺省值为1）。
+         * 
+         * @note 二维(x,y)坐标映射到三维坐标(x,y,z)，{0,1}表示(x,y)映射至(x,y,0)、
+         * {1,0}表示(x,y)映射至(y,x,0)、{0,2}表示(x,y)映射至(x,0,z)、{1,2}表示(x,y)映射至(0,y,z)。
+         */
+        template <typename T>
+        void import_from(const array<T> &elems, const array<vertex2dc> &nodes, int x_id = 0, int y_id = 1);
+
     protected:
         bool initialized_; // 类型是否已经初始化完成
 
@@ -557,6 +592,120 @@ namespace gctl
         node_map.clear();
         return;
     }
+
+    template <typename T>
+    void gctl::mesh_io::import_from(const array<T> &elems, const array<vertex3dc> &nodes)
+    {
+        reset(); // 重置网格数据
+
+        const std::type_info &tinfo = typeid(T);
+        element_type_enum oe_type = match_type(tinfo);
+        int vnum = elem_size(oe_type);
+
+        valid_node_size_ = nodes.size();
+        valid_elem_size_ = elems.size();
+
+        nodes_.resize(valid_node_size_);
+        for (size_t i = 0; i < valid_node_size_; i++)
+        {
+            nodes_[i].id = i;
+            nodes_[i].x = nodes[i].x;
+            nodes_[i].y = nodes[i].y;
+            nodes_[i].z = nodes[i].z;
+        }
+
+        elems_.resize(valid_elem_size_);
+        for (size_t i = 0; i < valid_elem_size_; i++)
+        {
+            elems_[i].id = i;
+            elems_[i].type = oe_type;
+
+            elems_[i].vert_ptrs.resize(vnum);
+            for (size_t v = 0; v < vnum; v++)
+            {
+                elems_[i].vert_ptrs[v] = nodes_.get(elems[i].vert[v]->id);
+            }
+        }
+
+        // 所有单元体都属于同一个组
+        valid_group_size_ = 1;
+        groups_.resize(1);
+        groups_[0].enabled = true;
+        groups_[0].type = oe_type;
+        groups_[0].phys_group = 0; // 默认组别为0
+        groups_[0].geom_group = 0; // 默认组别为0
+        groups_[0].part_group = 3;
+
+        for (size_t i = 0; i < elems_.size(); i++)
+        {
+            groups_[0].elem_ptrs.push_back(elems_.get(i));
+        }
+
+        groups_[0].enable_elements();
+        initialized_ = true;
+        return;
+    }
+
+    template <typename T>
+    void gctl::mesh_io::import_from(const array<T> &elems, const array<vertex2dc> &nodes, int x_id, int y_id)
+    {
+        reset(); // 重置网格数据
+
+        const std::type_info &tinfo = typeid(T);
+        element_type_enum oe_type = match_type(tinfo);
+        int vnum = elem_size(oe_type);
+
+        valid_node_size_ = nodes.size();
+        valid_elem_size_ = elems.size();
+
+        int xyz_ref[3];
+        nodes_.resize(valid_node_size_);
+        for (size_t i = 0; i < valid_node_size_; i++)
+        {
+            nodes_[i].id = i;
+
+            xyz_ref[0] = 0.0;
+            xyz_ref[1] = 0.0;
+            xyz_ref[2] = 0.0;
+            xyz_ref[x_id] = nodes[i].x;
+            xyz_ref[y_id] = nodes[i].y;
+
+            nodes_[i].x = xyz_ref[0];
+            nodes_[i].y = xyz_ref[1];
+            nodes_[i].z = xyz_ref[2];
+        }
+
+        elems_.resize(valid_elem_size_);
+        for (size_t i = 0; i < valid_elem_size_; i++)
+        {
+            elems_[i].id = i;
+            elems_[i].type = oe_type;
+
+            elems_[i].vert_ptrs.resize(vnum);
+            for (size_t v = 0; v < vnum; v++)
+            {
+                elems_[i].vert_ptrs[v] = nodes_.get(elems[i].vert[v]->id);
+            }
+        }
+
+        // 所有单元体都属于同一个组
+        valid_group_size_ = 1;
+        groups_.resize(1);
+        groups_[0].enabled = true;
+        groups_[0].type = oe_type;
+        groups_[0].phys_group = 0; // 默认组别为0
+        groups_[0].geom_group = 0; // 默认组别为0
+        groups_[0].part_group = 3;
+
+        for (size_t i = 0; i < elems_.size(); i++)
+        {
+            groups_[0].elem_ptrs.push_back(elems_.get(i));
+        }
+
+        groups_[0].enable_elements();
+        initialized_ = true;
+        return;
+    }
 };
 
 #endif // _GCTL_MESH_IO_H