[Getfem-commits] (no subject)

[Tetsuo Koyama]

branch: devel-tetsuo-outer_faces_in_ball
commit f391cabaf0ec9efe28b67b20a2f0cf91819bc235
Author: Tetsuo Koyama <address@hidden>
AuthorDate: Wed Feb 19 00:05:22 2020 +0900

    :heavy_plus_sign: outer_faces_in_ball
---
 interface/src/gf_mesh_get.cc                       | 46 +++++++++++++++++++++-
 .../demo_thermo_elasticity_electrical_coupling.py  |  2 +-
 2 files changed, 45 insertions(+), 3 deletions(-)

diff --git a/interface/src/gf_mesh_get.cc b/interface/src/gf_mesh_get.cc
index 037b182..78b0e84 100644
--- a/interface/src/gf_mesh_get.cc
+++ b/interface/src/gf_mesh_get.cc
@@ -161,10 +161,11 @@ outer_faces(const getfem::mesh &m, mexargs_in &in, 
mexargs_out &out, const std::
 
   bool with_normal(condition == "direction");
   bool with_box(condition == "box");
+  bool with_ball(condition == "ball");
 
   darray normal_vector;
-  bgeot::base_node un, pmin, pmax;
-  scalar_type threshold(0);
+  bgeot::base_node un, pmin, pmax, center;
+  scalar_type threshold(0), radius(0);
   if (with_normal) {
     normal_vector = in.pop().to_darray();
     scalar_type angle = in.pop().to_scalar();
@@ -178,6 +179,13 @@ outer_faces(const getfem::mesh &m, mexargs_in &in, 
mexargs_out &out, const std::
       pmin[k] = std::min(p1[k],p2[k]);
       pmax[k] = std::max(p1[k],p2[k]);
     }
+  } else if (with_ball) {
+    darray p1 = in.pop().to_darray(m.dim());
+    center.resize(m.dim());
+    for (size_type k=0; k < m.dim(); ++k) {
+      center[k] = p1[k];
+    }
+    radius = in.pop().to_scalar();
   }
 
   if (in.remaining()) cvlst = in.pop().to_bit_vector(&m.convex_index());
@@ -226,6 +234,27 @@ outer_faces(const getfem::mesh &m, mexargs_in &in, 
mexargs_out &out, const std::
             break;
           }
         }
+      } else if (with_ball) {
+        for (std::vector<size_type>::const_iterator pid=lst[i].ptid.begin();
+             pid != lst[i].ptid.end(); ++pid) {
+          if (!checked_pids.is_in(*pid)) {
+            checked_pids.add(*pid);
+            scalar_type checked_radius = scalar_type(0.0);
+            for (size_type k=0; k < m.dim(); ++k) {
+              checked_radius += pow(m.points()[*pid][k] - center[k], 2);
+            }
+            checked_radius = std::sqrt(checked_radius);
+            for (size_type k=0; k < m.dim(); ++k)
+              if (checked_radius > radius) {
+                rejected_pids.add(*pid);
+                break;
+              }
+          }
+          if (rejected_pids.is_in(*pid)) {
+            lst[i].cnt = -1;
+            break;
+          }
+        }
       }
       if (lst[i].cnt == 1) fcnt++;
     }
@@ -849,6 +878,19 @@ void gf_mesh_get(getfemint::mexargs_in& m_in,
        outer_faces(*pmesh, in, out, "box");
        );
 
+    /*@GET CVFIDs = ('outer faces in ball', @vec center, @scalar radius [, 
CVIDs])
+    Return the set of faces not shared by two convexes and lying within the 
ball of corresponding `center` and `radius`.
+
+    The output `CVFIDs` is a two-rows matrix, the first row lists convex
+    #ids, and the second one lists face numbers (local number in the
+    convex). If `CVIDs` is given, it returns portion of the boundary of
+    the convex set defined by the #ids listed in `CVIDs`.@*/
+    sub_command
+      ("outer faces in ball", 2, 3, 0, 1,
+       check_empty_mesh(pmesh);
+       outer_faces(*pmesh, in, out, "ball");
+       );
+
     /*@GET CVFIDs = ('adjacent face', @int cvid, @int fid)
     Return convex face of the neighbour element if it exists.
     If the convex have more than one neighbour
diff --git 
a/interface/tests/python/demo_thermo_elasticity_electrical_coupling.py 
b/interface/tests/python/demo_thermo_elasticity_electrical_coupling.py
index c26de45..8475c1d 100644
--- a/interface/tests/python/demo_thermo_elasticity_electrical_coupling.py
+++ b/interface/tests/python/demo_thermo_elasticity_electrical_coupling.py
@@ -129,7 +129,7 @@ mesh.region_subtract(BOTTOM_BOUND, HOLE_BOUND)
 mesh.region_merge(HOLE1_BOUND, HOLE2_BOUND)
 mesh.region_merge(HOLE1_BOUND, HOLE3_BOUND)
 
-assert mesh.region(HOLE_BOUND) == mesh.region(HOLE_BOUND1)
+np.testing.assert_array_equal(mesh.region(HOLE_BOUND), 
mesh.region(HOLE1_BOUND))
 
 if (export_mesh):
     mesh.export_to_vtk('mesh.vtk');