[Getfem-commits] (no subject)

[Konstantinos Poulios]

branch: merge-mesh-feature
commit 35a596ab9d82f81d35f31e0cb41d804109d41f97
Author: Konstantinos Poulios <address@hidden>
Date:   Fri Sep 7 20:23:40 2018 +0200

    Add merge_convexes_from_mesh method to the mesh object
---
 src/getfem/getfem_mesh.h |  7 +++++++
 src/getfem_mesh.cc       | 41 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 48 insertions(+)

diff --git a/src/getfem/getfem_mesh.h b/src/getfem/getfem_mesh.h
index 1748596..1c46be5 100644
--- a/src/getfem/getfem_mesh.h
+++ b/src/getfem/getfem_mesh.h
@@ -328,6 +328,13 @@ namespace getfem {
     void add_faces_of_convex(size_type)
     { GMM_ASSERT1(false, "Sorry, to be done"); }
 
+
+    /** Merge all convexes from a another mesh, possibly restricted to a
+        mesh region.
+    */
+    void merge_convexes_from_mesh(const mesh &m, size_type rg=size_type(-1),
+                                  scalar_type tol=scalar_type(0));
+
     /// Delete the convex of index ic from the mesh.
     void sup_convex(size_type ic, bool sup_points = false);
     /** Swap the indexes of the convex of indexes i and j
diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc
index 5a79eed..edf8d96 100644
--- a/src/getfem_mesh.cc
+++ b/src/getfem_mesh.cc
@@ -314,6 +314,47 @@ namespace getfem {
                            add_point(p3), add_point(p4));
   }
 
+  void mesh::merge_convexes_from_mesh(const mesh &msource, size_type rg,
+                                      scalar_type tol) {
+
+    size_type nbpt = points_index().last()+1;
+    GMM_ASSERT1(nbpt == nb_points(),
+                "Please call the optimize_structure() function before "
+                "merging elements from another mesh");
+    GMM_ASSERT1(rg == size_type(-1) || msource.region(rg).is_only_convexes(),
+                "The provided mesh region should only contain convexes");
+
+    const dal::bit_vector &convexes = (rg == size_type(-1))
+                                    ? msource.convex_index()
+                                    : msource.region(rg).index();
+    std::vector<size_type> old2new(msource.points_index().last()+1, 
size_type(-1));
+    for (dal::bv_visitor cv(convexes); !cv.finished(); ++cv) {
+
+      bgeot::pgeometric_trans pgt = msource.trans_of_convex(cv);
+      short_type nb = short_type(pgt->nb_points());
+      const ind_cv_ct &rct = msource.ind_points_of_convex(cv);
+      GMM_ASSERT1(nb == rct.size(), "Internal error");
+      std::vector<size_type> ind(nb);
+      for (short_type i = 0; i < nb; ++i) {
+        size_type old_pid = rct[i];
+        size_type new_pid = old2new[old_pid];
+        if (new_pid == size_type(-1)) {
+          size_type next_pid = points_index().last()+1;
+          base_node pt = msource.points()[old_pid];
+          new_pid = pts.add_node(pt, tol);
+          if (new_pid < next_pid && new_pid >= nbpt) {
+            // do not allow internal merging of nodes in the source mesh
+            new_pid = pts.add_node(pt, -1.);
+            GMM_ASSERT1(new_pid == next_pid, "Internal error");
+          }
+          old2new[old_pid] = new_pid;
+        }
+        ind[i] = new_pid;
+      }
+      add_convex(pgt, ind.begin());
+    }
+  }
+
   void mesh::sup_convex(size_type ic, bool sup_points) {
     static std::vector<size_type> ipt;
     if (sup_points) {