[Getfem-commits] (no subject)

[Yves Renard]

branch: master
commit f1fc36b58da6c236edb014d22b41795bd3c32da9
Author: Yves Renard <address@hidden>
Date:   Sat Jun 9 16:34:02 2018 +0200

    Making 'Normal' refering to the unit vector to a level-set when integrating 
on a level-set
---
 doc/sphinx/source/userdoc/gasm_high.rst            |  8 ++--
 .../getfem_generic_assembly_compile_and_exec.h     |  3 +-
 src/getfem/getfem_mesh_im_level_set.h              |  4 ++
 src/getfem_generic_assembly_compile_and_exec.cc    | 53 ++++++++++++++++++----
 src/getfem_global_function.cc                      | 17 ++++---
 src/getfem_mesh_im_level_set.cc                    | 31 +++++++++++++
 6 files changed, 96 insertions(+), 20 deletions(-)

diff --git a/doc/sphinx/source/userdoc/gasm_high.rst 
b/doc/sphinx/source/userdoc/gasm_high.rst
index 7a74a24..717d993 100644
--- a/doc/sphinx/source/userdoc/gasm_high.rst
+++ b/doc/sphinx/source/userdoc/gasm_high.rst
@@ -54,7 +54,7 @@ A specific language has been developed to describe the weak 
formulation of bound
 
   - ``X`` is the current coordinate on the real element, ``X(i)`` is its i-th 
component.
 
-  - ``Normal`` is the outward unit normal vector to a boundary (for 
integration on a domain boundary).
+  - ``Normal`` is the outward unit normal vector to a boundary, when 
integrating on a domain boundary, or the unit normal vector to a level-set when 
integrating on a level-set with a ``mesh_im_level_set`` method. In the latter 
case, the normal vector is in the direction of the level-set function gradient.
 
   - ``Reshape(t, i, j, ...)``: Reshape a vector/matrix/tensor. Note that all 
tensors in |gf| are stored in the Fortran order.
 
@@ -895,7 +895,7 @@ For |gf| 5.1. When using a fem cut by a level-set (using 
fem_level_set or mesh_f
   Xfem_minus(Test_Div_u)
   Xfem_minus(Test_Hess_u)
   
-which are only available when the evaluation (integration) is made on the 
curve/surface separating two zones of continuity, i.e. on the zero level-set of 
a considered level-set function (using a mesh_im_level_set object). For 
instance, a jump in the variable ``u`` will be given by::
+which are only available when the evaluation (integration) is made on the 
curve/surface separating two zones of continuity, i.e. on the zero level-set of 
a considered level-set function (using a ``mesh_im_level_set`` object). For 
instance, a jump in the variable ``u`` will be given by::
 
   Xfem_plus(u)-Xfem_minus(u)
 
@@ -903,7 +903,9 @@ and the average by::
 
   (Xfem_plus(u)+Xfem_minus(u))/2
 
-The value ``Xfem_plus(u)`` is the value of ``u`` on the side where the 
corresponding level-set function is positive and ``Xfem_minus(u)`` the value of 
``u`` on the side where the level-set function is negative.
+  The value ``Xfem_plus(u)`` is the value of ``u`` on the side where the 
corresponding level-set function is positive and ``Xfem_minus(u)`` the value of 
``u`` on the side where the level-set function is negative.
+
+  Additionally, note that, when integrating on a level-set with a 
``mesh_im_level_set`` object, ``Normal`` stands for the normal unit vector to 
the level-set in the direction of the gradient of the level-set function.
   
 Storage of sub-expressions in a getfem::im_data object during assembly
 ----------------------------------------------------------------------
diff --git a/src/getfem/getfem_generic_assembly_compile_and_exec.h 
b/src/getfem/getfem_generic_assembly_compile_and_exec.h
index 0f7766d..3802ca2 100644
--- a/src/getfem/getfem_generic_assembly_compile_and_exec.h
+++ b/src/getfem/getfem_generic_assembly_compile_and_exec.h
@@ -144,6 +144,7 @@ namespace getfem {
     struct region_mim_instructions {
 
       const mesh *m;
+      const mesh_im *im;
       ga_if_hierarchy current_hierarchy;
       std::map<std::string, base_vector> local_dofs;
       std::map<const mesh_fem *, pfem_precomp> pfps;
@@ -176,7 +177,7 @@ namespace getfem {
       ga_instruction_list instructions;
       std::map<scalar_type, std::list<pga_tree_node> > node_list;
 
-      region_mim_instructions(): m(0) {}
+    region_mim_instructions(): m(0), im(0) {}
     };
 
     std::list<ga_tree> trees; // The trees are stored mainly because they
diff --git a/src/getfem/getfem_mesh_im_level_set.h 
b/src/getfem/getfem_mesh_im_level_set.h
index 5ae02ef..623d267 100644
--- a/src/getfem/getfem_mesh_im_level_set.h
+++ b/src/getfem/getfem_mesh_im_level_set.h
@@ -111,6 +111,8 @@ namespace getfem {
       regular_simplex_pim = reg;
       base_singular_pim = sing;
     }
+
+    int location(void) const { return integrate_where; }
     
     size_type memsize() const {
       return mesh_im::memsize(); // + ... ;
@@ -171,6 +173,8 @@ namespace getfem {
     void set_level_set_boolean_operations(const std::string description) {
       ls_csg_description = description;
     }
+    void compute_normal_vector(const fem_interpolation_context &ctx,
+                              base_small_vector &vec) const;
   };
 
 
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc 
b/src/getfem_generic_assembly_compile_and_exec.cc
index 6afb52d..f54d6bd 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -19,6 +19,7 @@
 
 ===========================================================================*/
 
+#include "getfem/getfem_mesh_im_level_set.h"
 #include "getfem/getfem_generic_assembly_tree.h"
 #include "getfem/getfem_generic_assembly_semantic.h"
 #include "getfem/getfem_generic_assembly_compile_and_exec.h"
@@ -271,7 +272,7 @@ namespace getfem {
   struct ga_instruction_copy_Normal : public ga_instruction_copy_small_vect {
 
     virtual int exec() {
-      GA_DEBUG_INFO("Instruction: Normal");
+      GA_DEBUG_INFO("Instruction: unit normal vector");
       GMM_ASSERT1(t.size() == vec.size(), "Invalid outward unit normal "
                   "vector. Possible reasons: not on boundary or "
                   "transformation failed.");
@@ -283,6 +284,27 @@ namespace getfem {
       : ga_instruction_copy_small_vect(t_, Normal_)  {}
   };
 
+  struct ga_instruction_level_set_normal_vector : public ga_instruction {
+    base_tensor &t;
+    const mesh_im_level_set *mimls;
+    const fem_interpolation_context &ctx;
+    base_small_vector vec;
+    
+    virtual int exec() {
+      GA_DEBUG_INFO("Instruction: unit normal vector to a level-set");
+      mimls->compute_normal_vector(ctx, vec);
+      GMM_ASSERT1(t.size() == vec.size(), "Invalid outward unit normal "
+                  "vector. Possible reasons: not on boundary or "
+                  "transformation failed.");
+      gmm::copy(vec, t.as_vector());
+      return 0;
+    }
+    ga_instruction_level_set_normal_vector
+    (base_tensor &t_, const mesh_im_level_set *mimls_,
+     const fem_interpolation_context &ctx_)
+      : t(t_), mimls(mimls_), ctx(ctx_), vec(t.size())  {}
+  };
+
   struct ga_instruction_element_size : public ga_instruction {
     base_tensor &t;
     scalar_type &es;
@@ -4932,13 +4954,24 @@ namespace getfem {
       break;
 
     case GA_NODE_NORMAL:
-      GMM_ASSERT1(!function_case,
-                  "No use of Normal is allowed in functions");
-      if (pnode->tensor().size() != m.dim())
-        pnode->init_vector_tensor(m.dim());
-      pgai = std::make_shared<ga_instruction_copy_Normal>
-             (pnode->tensor(), gis.Normal);
-      rmi.instructions.push_back(std::move(pgai));
+      {
+       GMM_ASSERT1(!function_case,
+                   "No use of Normal is allowed in functions");
+       if (pnode->tensor().size() != m.dim())
+         pnode->init_vector_tensor(m.dim());
+       const mesh_im_level_set *mimls
+         = dynamic_cast<const mesh_im_level_set *>(rmi.im);
+       if (mimls && mimls->location()==mesh_im_level_set::INTEGRATE_BOUNDARY) {
+         // Appel avec ctx (pt de Gauss)
+         pgai = std::make_shared<ga_instruction_level_set_normal_vector>
+           (pnode->tensor(), mimls, gis.ctx);
+         rmi.instructions.push_back(std::move(pgai));
+       } else {
+         pgai = std::make_shared<ga_instruction_copy_Normal>
+           (pnode->tensor(), gis.Normal);
+         rmi.instructions.push_back(std::move(pgai));
+       }
+      }
       break;
 
     case GA_NODE_INTERPOLATE_X:
@@ -6349,7 +6382,7 @@ namespace getfem {
           ga_instruction_set::region_mim rm(td.mim, td.rg);
           ga_instruction_set::region_mim_instructions &rmi
             = gis.whole_instructions[rm];
-          rmi.m = td.m;
+          rmi.m = td.m; rmi.im = td.mim;
           // rmi.interpolate_infos.clear();
           ga_compile_interpolate_trans(root, workspace, gis, rmi, *(td.m));
           ga_compile_node(root, workspace, gis,rmi, *(td.m), false,
@@ -6399,7 +6432,7 @@ namespace getfem {
             ga_instruction_set::region_mim rm(td.mim, td.rg);
             ga_instruction_set::region_mim_instructions &rmi
               = gis.whole_instructions[rm];
-            rmi.m = td.m;
+            rmi.m = td.m; rmi.im = td.mim;
             // rmi.interpolate_infos.clear();
             ga_compile_interpolate_trans(root, workspace, gis, rmi, *(td.m));
             ga_compile_node(root, workspace, gis, rmi, *(td.m), false,
diff --git a/src/getfem_global_function.cc b/src/getfem_global_function.cc
index d443a25..fcdda89 100644
--- a/src/getfem_global_function.cc
+++ b/src/getfem_global_function.cc
@@ -31,24 +31,24 @@ namespace getfem {
   scalar_type global_function_simple::val
   (const fem_interpolation_context &c) const {
     base_node pt = c.xreal();
-    GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") 
" <<
-                                   "passed to a global function of dim = "<< 
dim_ <<".");
+    GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
+               << "passed to a global function of dim = "<< dim_ <<".");
     return this->val(pt);
   }
 
   void global_function_simple::grad
   (const fem_interpolation_context &c, base_small_vector &g) const {
     base_node pt = c.xreal();
-    GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") 
" <<
-                                   "passed to a global function of dim = "<< 
dim_ <<".");
+    GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
+               << "passed to a global function of dim = "<< dim_ <<".");
     this->grad(pt, g);
   }
 
   void global_function_simple::hess
   (const fem_interpolation_context &c, base_matrix &h) const {
     base_node pt = c.xreal();
-    GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") 
" <<
-                                   "passed to a global function of dim = "<< 
dim_ <<".");
+    GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
+               << "passed to a global function of dim = "<< dim_ <<".");
     this->hess(pt, h);
   }
 
@@ -744,6 +744,11 @@ namespace getfem {
       scalar_type y = (*mls_y)(c.xref());
       if (c.xfem_side() > 0 && y <= 0) y = 1E-13;
       if (c.xfem_side() < 0 && y >= 0) y = -1E-13;
+      // if (c.xfem_side()) {
+      //       cout << "point : " << c.xref() << " side : " << c.xfem_side() 
<< endl;
+      //       cout << "side -1 : " << fn->val(x,-1E-13) << " : " << 
fn->val(x,-1E-16) << endl;
+      //       cout << "side  1 : " << fn->val(x, 1E-13) << " : " << 
fn->val(x, 1E-16) << endl;
+      // }
 
       return fn->val(x,y);
     }
diff --git a/src/getfem_mesh_im_level_set.cc b/src/getfem_mesh_im_level_set.cc
index a030c0a..6fd7afd 100644
--- a/src/getfem_mesh_im_level_set.cc
+++ b/src/getfem_mesh_im_level_set.cc
@@ -441,6 +441,35 @@ namespace getfem {
     // cout << "Number of built methods : " << build_methods.size() << endl;
   }
 
+  void mesh_im_level_set::compute_normal_vector
+  (const fem_interpolation_context &ctx, base_small_vector &vec) const {
+    size_type nb_ls = mls->nb_level_sets(), j = 0;
+    std::vector<pmesher_signed_distance> mesherls0(nb_ls);
+    if (vec.size() != linked_mesh().dim()) vec.resize(linked_mesh().dim());
+    base_small_vector un(ctx.pgt()->dim());
+                        
+    if (nb_ls == 0) {
+      gmm::clear(vec); return; 
+    } else if (nb_ls == 1) {
+      mesherls0[0]
+       = mls->get_level_set(0)->mls_of_convex(ctx.convex_num(), 0, false);
+    } else {
+      scalar_type d_min(0);
+      for (unsigned i = 0; i < nb_ls; ++i) {
+       mesherls0[i]
+         = mls->get_level_set(i)->mls_of_convex(ctx.convex_num(), 0, false);
+       scalar_type d = gmm::abs((*(mesherls0[i]))(ctx.xref()));
+       if (i == 0 || d < d_min) { d_min = d; j = i; }
+      }
+    }
+    mesherls0[j]->grad(ctx.xref(), un);
+    gmm::mult(ctx.B(), un, vec);
+    scalar_type no = gmm::vect_norm2(vec);
+    if (no != scalar_type(0))
+      gmm::scale(vec, scalar_type(1) / gmm::vect_norm2(vec));
+  }
+
+
 
   void mesh_im_cross_level_set::update_from_context(void) const
   { is_adapted = false; }
@@ -713,6 +742,8 @@ namespace getfem {
     is_adapted = true; touch();
   }
 
+
+
   
 }  /* end of namespace getfem.                                             */