[Getfem-commits] (no subject)

[Tetsuo Koyama]

branch: devel-tetsuo-outer_faces_in_ball
commit 786709b28317abaa9c56ef281d7ecc9dfb5646dc
Author: Tetsuo Koyama <address@hidden>
AuthorDate: Sun Feb 23 16:14:20 2020 +0900

    :heavy_plus_sign: select_faces_in_ball
---
 src/getfem/getfem_mesh.h                       |  8 ++++++++
 src/getfem_mesh.cc                             | 24 ++++++++++++++++++++++++
 tests/thermo_elasticity_electrical_coupling.cc | 18 ++++++++++++++----
 3 files changed, 46 insertions(+), 4 deletions(-)

diff --git a/src/getfem/getfem_mesh.h b/src/getfem/getfem_mesh.h
index 49ce06e..d6fce2a 100644
--- a/src/getfem/getfem_mesh.h
+++ b/src/getfem/getfem_mesh.h
@@ -678,6 +678,14 @@ namespace getfem {
                       const base_node &pt1,
                       const base_node &pt2);
 
+  /** Select in the region mr the faces of the mesh m lying entirely in the
+      ball delimated by pt1 and radius.
+   */
+  mesh_region APIDECL
+  select_faces_in_ball(const mesh &m, const mesh_region &mr,
+                       const base_node &center,
+                       scalar_type radius);
+
   mesh_region APIDECL
   select_convexes_in_box(const mesh &m, const mesh_region &mr,
                          const base_node &pt1,
diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc
index 8ba79fa..5347b78 100644
--- a/src/getfem_mesh.cc
+++ b/src/getfem_mesh.cc
@@ -957,6 +957,30 @@ namespace getfem {
     return mrr;
   }
 
+  mesh_region select_faces_in_ball(const mesh &m, const mesh_region &mr,
+                                   const base_node &center, scalar_type 
radius) {
+    mesh_region mrr;
+    size_type N = m.dim();
+    GMM_ASSERT1(center.size() == N, "Wrong dimensions");
+    for (getfem::mr_visitor i(mr, m); !i.finished(); ++i)
+      if (i.f() != short_type(-1)) {
+        bgeot::mesh_structure::ind_pt_face_ct pt
+          = m.ind_points_of_face_of_convex(i.cv(), i.f());
+
+        bool is_in = true;
+        for (bgeot::mesh_structure::ind_pt_face_ct::iterator it = pt.begin();
+             it != pt.end(); ++it) {
+          scalar_type checked_radius = scalar_type(0.0);
+          for (size_type j = 0; j < N; ++j)
+            checked_radius += pow(m.points()[*it][j] - center[j], 2);
+          checked_radius = std::sqrt(checked_radius);
+          if (checked_radius > radius) { is_in = false; break; }
+        }
+        if (is_in) mrr.add(i.cv(), i.f());
+      }
+    return mrr;
+  }
+
   mesh_region select_convexes_in_box(const mesh &m, const mesh_region &mr,
                                      const base_node &pt1, const base_node 
&pt2) {
     mesh_region mrr;
diff --git a/tests/thermo_elasticity_electrical_coupling.cc 
b/tests/thermo_elasticity_electrical_coupling.cc
index a9968ce..df15913 100644
--- a/tests/thermo_elasticity_electrical_coupling.cc
+++ b/tests/thermo_elasticity_electrical_coupling.cc
@@ -171,10 +171,20 @@ int main(int argc, char *argv[]) {
   mesh.region( HOLE2_BOUND) = fb7;
   mesh.region( HOLE3_BOUND) = fb8;
 
-  mesh.region_merge(HOLE1_BOUND, HOLE2_BOUND);
-  mesh.region_merge(HOLE1_BOUND, HOLE3_BOUND);
-
-  GMM_ASSERT1(mesh.region(HOLE_BOUND) == mesh.region(HOLE1_BOUND), "Region is 
not equal.");
+  dal::bit_vector nn = mesh.convex_index();
+  bgeot::size_type i;
+  for (i << nn; i != bgeot::size_type(-1); i << nn) {
+    bgeot::pconvex_structure cvs = mesh.structure_of_convex(i);
+    for (bgeot::short_type f = 0; f < cvs->nb_faces(); ++f) {
+      if(mesh.region(HOLE_BOUND).is_in(i, f))
+        GMM_ASSERT1(
+          mesh.region(HOLE1_BOUND).is_in(i, f) ||
+          mesh.region(HOLE2_BOUND).is_in(i, f) ||
+          mesh.region(HOLE3_BOUND).is_in(i, f),
+          "Error in select region"
+        );
+    }
+  }
 
   if (export_mesh) {
     getfem::vtk_export exp("mesh.vtk", false);