[Getfem-commits] r4627 - in /trunk/getfem: doc/sphinx/source/userdoc/ src/ src/getfem/

[Yves . Renard]

Author: renard
Date: Fri Apr 25 16:26:26 2014
New Revision: 4627

URL: http://svn.gna.org/viewcvs/getfem?rev=4627&view=rev
Log:
generic interpolation

Modified:
    trunk/getfem/doc/sphinx/source/userdoc/index.rst
    trunk/getfem/doc/sphinx/source/userdoc/interMM.rst
    trunk/getfem/src/getfem/getfem_generic_assembly.h
    trunk/getfem/src/getfem_generic_assembly.cc

Modified: trunk/getfem/doc/sphinx/source/userdoc/index.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/index.rst?rev=4627&r1=4626&r2=4627&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/index.rst    (original)
+++ trunk/getfem/doc/sphinx/source/userdoc/index.rst    Fri Apr 25 16:26:26 2014
@@ -23,6 +23,7 @@
    asm
    gasm_high
    gasm_low
+   interMM
    ifem
    iinteg
    xfem
@@ -30,7 +31,6 @@
    computeL2H1
    computeD
    export
-   interMM
    convect
    model
    examples

Modified: trunk/getfem/doc/sphinx/source/userdoc/interMM.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/interMM.rst?rev=4627&r1=4626&r2=4627&view=diff
==============================================================================
--- trunk/getfem/doc/sphinx/source/userdoc/interMM.rst  (original)
+++ trunk/getfem/doc/sphinx/source/userdoc/interMM.rst  Fri Apr 25 16:26:26 2014
@@ -6,8 +6,13 @@
 
 .. _ud-intermm:
 
-Interpolation on different meshes
-=================================
+Interpolation of arbitary quantities
+====================================
+
+Once a solution has been computed, it is quite easy to extract any quantity of 
interest on it with the interpolation functions for instance for post-treatment.
+
+Basic interpolation
+*******************
 
 The file :file:`getfem/getfem_interpolation.h` defines the function 
 ``getfem::interpolation(...)`` to interpolate a solution from a given 
mesh/finite 
@@ -44,3 +49,11 @@
 the interpolation is done with a simple matrix multiplication::
 
   gmm::mult(M, U, V);
+
+
+Interpolation based on the high-level generic assembly langage
+**************************************************************
+
+It is possible to extract some arbitraries expressions linking several fields 
thanks to the high-level generic assembly langage and the interpolation 
functions.
+
+This is specially dedicated to the model object. For instance if ``md`` is a 
valid object containing some defined variables ``u`` and ``p`` ...

Modified: trunk/getfem/src/getfem/getfem_generic_assembly.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_generic_assembly.h?rev=4627&r1=4626&r2=4627&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_generic_assembly.h   (original)
+++ trunk/getfem/src/getfem/getfem_generic_assembly.h   Fri Apr 25 16:26:26 2014
@@ -418,7 +418,8 @@
   struct ga_interpolation_context {
 
     virtual const bgeot::stored_point_tab &
-    points_for_element(size_type cv, short_type f) const = 0;
+    points_for_element(size_type cv, short_type f,
+                       std::vector<size_type> &ind) const = 0;
     virtual bool use_pgp(size_type cv) const = 0;
     virtual void store_result(size_type cv, size_type i, base_tensor &t) = 0;
     virtual void finalize(void) = 0;
@@ -442,12 +443,18 @@
   (const getfem::model &md, const std::string &expr, const mesh_fem &mf,
    base_vector &result, const mesh_region &rg=mesh_region::all_convexes());
 
+  // untested
   void ga_interpolation_mti
   (const getfem::model &md, const std::string &expr, mesh_trans_inv &mti,
    base_vector &result, int extrapolation = 0,
    const mesh_region &rg=mesh_region::all_convexes(),
    size_type nbdof_ = size_type(-1));
 
+  // untested
+  void ga_interpolation_im_data
+  (const getfem::model &md, const std::string &expr, im_data &imd,
+   base_vector &result, const mesh_region &rg=mesh_region::all_convexes());
+
 
 
 }  /* end of namespace getfem.                                             */

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4627&r1=4626&r2=4627&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Fri Apr 25 16:26:26 2014
@@ -6320,11 +6320,13 @@
       ga_instruction_list &gil = it->second.instructions;
 
       // iteration on elements (or faces of elements)
+      std::vector<size_type> ind;
       for (getfem::mr_visitor v(region, m); !v.finished(); ++v) {
+        ind.resize(0);
         const bgeot::stored_point_tab &spt
-          = gic.points_for_element(v.cv(), v.f());
-
-        if (spt.size()) {
+          = gic.points_for_element(v.cv(), v.f(), ind);
+
+        if (&spt && ind.size() && spt.size()) {
           bgeot::vectors_to_base_matrix(G, m.points_of_convex(v.cv()));
           size_type N = G.nrows();
           bgeot::pgeometric_trans pgt = m.trans_of_convex(v.cv());
@@ -6333,7 +6335,7 @@
             gis.ctx = fem_interpolation_context(gis.ctx.pgp(), 0, 0, G,
                                                 v.cv(), v.f());
           } else {
-            if (gic.use_pgp(v.cv())) {
+            if (!(gic.use_pgp(v.cv()))) {
               gis.ctx = fem_interpolation_context(pgt, 0, spt[0], G,
                                                   v.cv(), v.f());
             } else {
@@ -6345,11 +6347,12 @@
           
           // iterations on interpolation points
           gis.nbpt = spt.size();
-          for (gis.ipt = 0; gis.ipt < gis.nbpt; ++(gis.ipt)) {
+          for (size_type ii = 0; ii < ind.size(); ++ii) {
+            gis.ipt = ind[ii];
             if (gis.ctx.have_pgp()) gis.ctx.set_ii(gis.ipt);
             else gis.ctx.set_xref(spt[gis.ipt]);
 
-            if (gis.ipt == 0 || !(pgt->is_linear())) {
+            if (ii == 0 || !(pgt->is_linear())) {
               // Computation of unit normal vector in case of a boundary
               if (v.f() != short_type(-1)) {
                 const base_matrix& B = gis.ctx.B();
@@ -6453,7 +6456,7 @@
   // Interpolation functions
   //=========================================================================
 
-
+  // general Interpolation 
   void ga_interpolation(ga_workspace &workspace,
                         ga_interpolation_context &gic) {
     ga_instruction_set gis;
@@ -6461,7 +6464,7 @@
     ga_interpolation_exec(gis, workspace, gic);
   }
 
-  
+  // Interpolation on a Lagrange fem on the same mesh
   struct ga_interpolation_context_fem_same_mesh
     : public ga_interpolation_context {
     base_vector &result;
@@ -6471,15 +6474,24 @@
     size_type s;
 
     virtual const bgeot::stored_point_tab &
-    points_for_element(size_type cv, short_type f) const {
+    points_for_element(size_type cv, short_type f,
+                       std::vector<size_type> &ind) const {
       pfem pf = mf.fem_of_element(cv);
       GMM_ASSERT1(pf->is_lagrange(),
                   "Only Lagrange fems can be used in interpolation");
+      
       if (f != short_type(-1)) {
-        return *(store_point_tab(pf->node_convex(cv).points_of_face(f)));
+        
+        for (size_type i = 0;
+             i < pf->node_convex(cv).structure()->nb_points_of_face(f); ++i)
+          ind.push_back
+            (pf->node_convex(cv).structure()->ind_points_of_face(f)[i]);
       } else {
-        return *(pf->node_tab(cv));
-      }
+        for (size_type i = 0; i < pf->node_convex(cv).nb_points(); ++i)
+          ind.push_back(i);
+      }
+      
+      return *(pf->node_tab(cv));
     }
 
     virtual bool use_pgp(size_type) const { return true; }
@@ -6540,7 +6552,7 @@
     ga_interpolation_Lagrange_fem(workspace, mf, result);
   }
 
-
+  // Interpolation on a cloud of points
   struct ga_interpolation_context_mti
     : public ga_interpolation_context {
     base_vector &result;
@@ -6550,12 +6562,15 @@
 
 
     virtual const bgeot::stored_point_tab &
-    points_for_element(size_type cv, short_type) const {
+    points_for_element(size_type cv, short_type,
+                       std::vector<size_type> &ind) const {
       std::vector<size_type> itab;
       mti.points_on_convex(cv, itab);
       std::vector<base_node> pt_tab(itab.size());
-      for (size_type i = 0; i < itab.size(); ++i)
+      for (size_type i = 0; i < itab.size(); ++i) {
         pt_tab[i] = mti.reference_coords()[itab[i]];
+        ind.push_back(i);
+      }
       return *(store_point_tab(pt_tab));
     }
 
@@ -6608,4 +6623,75 @@
   }
 
 
+  // Interpolation on a im_data
+  struct ga_interpolation_context_im_data
+    : public ga_interpolation_context {
+    base_vector &result;
+    im_data &imd;
+    bool initialized;
+    size_type s;
+
+    virtual const bgeot::stored_point_tab &
+    points_for_element(size_type cv, short_type /*f*/,
+                       std::vector<size_type> &ind) const {
+      pintegration_method pim =imd.get_mesh_im().int_method_of_element(cv);
+      if (pim->type() == IM_NONE) return *(bgeot::pstored_point_tab(0));
+      GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods cannot "
+                  "be used in high level generic assembly");
+      for (size_type i = 0;
+           i < pim->approx_method()->nb_points_on_convex(); ++i)
+        ind.push_back(i);
+      return pim->approx_method()->integration_points();
+    }
+
+    virtual bool use_pgp(size_type cv) const {
+      pintegration_method pim =imd.get_mesh_im().int_method_of_element(cv);
+      if (pim->type() == IM_NONE) return false;
+      GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods cannot "
+                  "be used in high level generic assembly");
+      return !(pim->approx_method()->is_built_on_the_fly());
+    }
+
+    virtual void store_result(size_type cv, size_type i, base_tensor &t) {
+      size_type si = t.size();
+      if (!initialized) {
+        s = si;
+        imd.set_tensor_size(t.sizes());
+        gmm::resize(result, s * imd.nb_filtered_index());
+        gmm::clear(result);
+        initialized = true;
+      }
+      GMM_ASSERT1(s == si, "Internal error");
+      size_type ipt = imd.filtered_index_of_point(cv, i);
+      gmm::add(t.as_vector(),
+               gmm::sub_vector(result, gmm::sub_interval(s*ipt, s)));
+    }
+    
+    virtual void finalize(void)
+    { if (!initialized) gmm::clear(result); }
+
+    virtual const mesh &linked_mesh(void)
+    { return imd.get_mesh_im().linked_mesh(); }
+
+    ga_interpolation_context_im_data(im_data &imd_, base_vector &r)
+      : result(r), imd(imd_), initialized(false) { }
+  };
+
+
+  // To be parallelized ...
+  void ga_interpolation_im_data
+  (const getfem::model &md, const std::string &expr, im_data &imd,
+   base_vector &result, const mesh_region &rg) {
+
+    ga_workspace workspace(md);
+    workspace.add_interpolation_expression
+      (expr, imd.get_mesh_im().linked_mesh(), rg);
+    
+    ga_interpolation_context_im_data gic(imd, result);
+    ga_interpolation(workspace, gic);
+  }
+  
+
+
+
 } /* end of namespace */