[Getfem-commits] r5350 - in /trunk/getfem/src: ./ getfem/

[logari81]

Author: logari81
Date: Fri Jul 22 16:43:34 2016
New Revision: 5350

URL: http://svn.gna.org/viewcvs/getfem?rev=5350&view=rev
Log:
allow fem precomputations for fem defined on the real element and implement 
precomputations for fem_global_function

Modified:
    trunk/getfem/src/getfem/getfem_fem.h
    trunk/getfem/src/getfem/getfem_fem_global_function.h
    trunk/getfem/src/getfem/getfem_mesh_fem_global_function.h
    trunk/getfem/src/getfem_fem.cc
    trunk/getfem/src/getfem_fem_global_function.cc
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/src/getfem_mesh_fem_global_function.cc

Modified: trunk/getfem/src/getfem/getfem_fem.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_fem.h?rev=5350&r1=5349&r2=5350&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_fem.h        (original)
+++ trunk/getfem/src/getfem/getfem_fem.h        Fri Jul 22 16:43:34 2016
@@ -226,7 +226,8 @@
   class fem_interpolation_context;
 
   /** @brief Base class for finite element description */
-  class virtual_fem : virtual public dal::static_stored_object {
+  class virtual_fem : virtual public dal::static_stored_object,
+                      public std::enable_shared_from_this<const virtual_fem> {
   public :
     enum vec_type { VECTORIAL_NOTRANSFORM_TYPE, VECTORIAL_PRIMAL_TYPE,
                     VECTORIAL_DUAL_TYPE };
@@ -458,7 +459,9 @@
       vtype = VECTORIAL_NOTRANSFORM_TYPE;
       cvs_node = bgeot::new_convex_structure();
     }
-    virtual_fem(const virtual_fem& f) : dal::static_stored_object()
+    virtual_fem(const virtual_fem& f)
+      : dal::static_stored_object(),
+        std::enable_shared_from_this<const virtual_fem>()
     { copy(f); DAL_STORED_OBJECT_DEBUG_CREATED(this, "Fem"); }
     virtual ~virtual_fem() { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Fem"); }
   private:

Modified: trunk/getfem/src/getfem/getfem_fem_global_function.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_fem_global_function.h?rev=5350&r1=5349&r2=5350&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_fem_global_function.h        (original)
+++ trunk/getfem/src/getfem/getfem_fem_global_function.h        Fri Jul 22 
16:43:34 2016
@@ -56,6 +56,9 @@
     mutable std::vector<std::vector<size_type> > index_of_global_dof_;
     mutable bgeot::pstored_point_tab pspt_override;
 
+    mutable std::vector< std::map<bgeot::pstored_point_tab, 
std::vector<base_tensor> > >
+      precompval, precompgrad, precomphess; // store values, gradients and 
hessians of base functions
+
     void init();
     virtual void update_from_context() const;
   public :

Modified: trunk/getfem/src/getfem/getfem_mesh_fem_global_function.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_mesh_fem_global_function.h?rev=5350&r1=5349&r2=5350&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_mesh_fem_global_function.h   (original)
+++ trunk/getfem/src/getfem/getfem_mesh_fem_global_function.h   Fri Jul 22 
16:43:34 2016
@@ -52,10 +52,6 @@
     getfem::pfem fem_;
   public :
 
-    void set_functions(pglobal_function f,
-                       const mesh_im &mim=dummy_mesh_im());
-    void set_functions(pglobal_function f1, pglobal_function f2,
-                       const mesh_im &mim=dummy_mesh_im());
     void set_functions(const std::vector<pglobal_function>& f,
                        const mesh_im &mim=dummy_mesh_im());
     // size_type memsize() const;

Modified: trunk/getfem/src/getfem_fem.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_fem.cc?rev=5350&r1=5349&r2=5350&view=diff
==============================================================================
--- trunk/getfem/src/getfem_fem.cc      (original)
+++ trunk/getfem/src/getfem_fem.cc      Fri Jul 22 16:43:34 2016
@@ -65,6 +65,18 @@
     return (face_num_ != short_type(-1));
   }
 
+  // Comment regarding fem defined on the real element:
+  // Precomputed values, gradients and hessians of fem defined on the real
+  // element are not dealt with within fem_interpolation_context.
+  // In that case, any precomputations can be performed within the fem
+  // itself, which has to store the corresponding values internally.
+  // The methods real_base_value, real_grad_base_value or real_hess_base_value
+  // of the fem have the possibility to check if the passed ctx has a pfp
+  // and extract the corresponding internally stored results based on
+  // ctx.convex_num(), ctx.pfp()->get_ppoint_tab() and ctx.ii().
+  // In that case, the storage available in ctx.pfp()->c, ctx.pfp()->pc
+  // and ctx.pfp()->hpc is not used.
+
   void fem_interpolation_context::base_value(base_tensor& t,
                                              bool withM) const {
     if (pf()->is_on_real_element())
@@ -77,7 +89,7 @@
           t.mat_transp_reduction(pfp_->val(ii()), K(), 1); break;
         case virtual_fem::VECTORIAL_DUAL_TYPE:
           t.mat_transp_reduction(pfp_->val(ii()), B(), 1); break;
-        default: t=pfp_->val(ii());
+        default: t = pfp_->val(ii());
         }
       }
       else {

Modified: trunk/getfem/src/getfem_fem_global_function.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_fem_global_function.cc?rev=5350&r1=5349&r2=5350&view=diff
==============================================================================
--- trunk/getfem/src/getfem_fem_global_function.cc      (original)
+++ trunk/getfem/src/getfem_fem_global_function.cc      Fri Jul 22 16:43:34 2016
@@ -67,6 +67,10 @@
   }
 
   void fem_global_function::update_from_context() const {
+
+    precompval.clear();
+    precompgrad.clear();
+    precomphess.clear();
 
     mib.resize(2);
     mib[0] = short_type(1);
@@ -183,12 +187,30 @@
     size_type nbdof = nb_dof(cv);
     mib[0] = short_type(nbdof);
     t.adjust_sizes(mib);
-    for (size_type i=0; i < nbdof; ++i) {
-      /*cerr << "fem_global_function: real_base_value(" << c.xreal() << ")\n";
-      if (c.have_G()) cerr << "G = " << c.G() << "\n";
-      else cerr << "no G\n";*/
-      t[i] = functions[index_of_global_dof_[cv][i]]->val(c);
-    }
+    if (c.have_pfp() && c.ii() != size_type(-1)) {
+      if (precompval.size() == 0)
+        precompval.resize(m.convex_index().last_true()+1);
+      const bgeot::pstored_point_tab ptab = c.pfp()->get_ppoint_tab();
+      auto it = precompval[cv].find(ptab);
+      if (it == precompval[cv].end()) {
+        it = precompval[cv]
+             .emplace(ptab, std::vector<base_tensor>(ptab->size())).first;
+        base_matrix G;
+        bgeot::vectors_to_base_matrix(G, m.points_of_convex(cv));
+        for (size_type k = 0; k < ptab->size(); ++k) {
+          const fem_interpolation_context
+            ctx(m.trans_of_convex(cv), shared_from_this(), (*ptab)[k], G, cv);
+          real_base_value(ctx, it->second[k]);
+        }
+      }
+      gmm::copy(it->second[c.ii()].as_vector(), t.as_vector());
+    } else
+      for (size_type i=0; i < nbdof; ++i) {
+        /*cerr << "fem_global_function: real_base_value(" << c.xreal() << 
")\n";
+        if (c.have_G()) cerr << "G = " << c.G() << "\n";
+        else cerr << "no G\n";*/
+        t[i] = functions[index_of_global_dof_[cv][i]]->val(c);
+      }
   }
 
   void fem_global_function::real_grad_base_value
@@ -198,11 +220,30 @@
     size_type nbdof = nb_dof(cv);
     mig[0] = short_type(nbdof);
     t.adjust_sizes(mig);
-    base_small_vector G(dim());
-    for (size_type i=0; i < nbdof; ++i) {
-      functions[index_of_global_dof_[cv][i]]->grad(c,G);
-      for (size_type j=0; j < dim(); ++j)
-        t[j*nbdof + i] = G[j];
+    if (c.have_pfp() && c.ii() != size_type(-1)) {
+      if (precompgrad.size() == 0)
+        precompgrad.resize(m.convex_index().last_true()+1);
+      const bgeot::pstored_point_tab ptab = c.pfp()->get_ppoint_tab();
+      auto it = precompgrad[cv].find(ptab);
+      if (it == precompgrad[cv].end()) {
+        it = precompgrad[cv]
+             .emplace(ptab, std::vector<base_tensor>(ptab->size())).first;
+        base_matrix G;
+        bgeot::vectors_to_base_matrix(G, m.points_of_convex(cv));
+        for (size_type k = 0; k < ptab->size(); ++k) {
+          const fem_interpolation_context
+            ctx(m.trans_of_convex(cv), shared_from_this(), (*ptab)[k], G, cv);
+          real_grad_base_value(ctx, it->second[k]);
+        }
+      }
+      gmm::copy(it->second[c.ii()].as_vector(), t.as_vector());
+    } else {
+      base_small_vector G(dim());
+      for (size_type i=0; i < nbdof; ++i) {
+        functions[index_of_global_dof_[cv][i]]->grad(c,G);
+        for (size_type j=0; j < dim(); ++j)
+          t[j*nbdof + i] = G[j];
+      }
     }
   }
 
@@ -213,13 +254,33 @@
     size_type nbdof = nb_dof(cv);
     mih[0] = short_type(nbdof);
     t.adjust_sizes(mih);
-    base_matrix H(dim(),dim());
-    for (size_type i=0; i < nbdof; ++i) {
-      functions[index_of_global_dof_[cv][i]]->hess(c,H);
-      for (size_type jk=0; jk < size_type(dim()*dim()); ++jk)
-        t[jk*nbdof + i] = H[jk];
+    if (c.have_pfp() && c.ii() != size_type(-1)) {
+      if (precomphess.size() == 0)
+        precomphess.resize(m.convex_index().last_true()+1);
+      const bgeot::pstored_point_tab ptab = c.pfp()->get_ppoint_tab();
+      auto it = precomphess[cv].find(ptab);
+      if (it == precomphess[cv].end()) {
+        it = precomphess[cv]
+             .emplace(ptab, std::vector<base_tensor>(ptab->size())).first;
+        base_matrix G;
+        bgeot::vectors_to_base_matrix(G, m.points_of_convex(cv));
+        for (size_type k = 0; k < ptab->size(); ++k) {
+          const fem_interpolation_context
+            ctx(m.trans_of_convex(cv), shared_from_this(), (*ptab)[k], G, cv);
+          real_hess_base_value(ctx, it->second[k]);
+        }
+      }
+      gmm::copy(it->second[c.ii()].as_vector(), t.as_vector());
+    } else {
+      base_matrix H(dim(),dim());
+      for (size_type i=0; i < nbdof; ++i) {
+        functions[index_of_global_dof_[cv][i]]->hess(c,H);
+        for (size_type jk=0; jk < size_type(dim()*dim()); ++jk)
+          t[jk*nbdof + i] = H[jk];
+      }
     }
   }
+
 
   DAL_SIMPLE_KEY(special_fem_globf_key, pfem);
 

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5350&r1=5349&r2=5350&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Fri Jul 22 16:43:34 2016
@@ -2717,11 +2717,7 @@
         pfem pf = mf.fem_of_element(ctx.convex_num());
         if (!pfp || pf != pfp->get_pfem() ||
             ctx.pgp()->get_ppoint_tab() != pfp->get_ppoint_tab()) {
-          if (pf->is_on_real_element())
-            pfp = 0;
-          else {
-            pfp = fp_pool(pf, ctx.pgp()->get_ppoint_tab());
-          }
+          pfp = fp_pool(pf, ctx.pgp()->get_ppoint_tab());
         }
       } else {
         pfp = 0;

Modified: trunk/getfem/src/getfem_mesh_fem_global_function.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_mesh_fem_global_function.cc?rev=5350&r1=5349&r2=5350&view=diff
==============================================================================
--- trunk/getfem/src/getfem_mesh_fem_global_function.cc (original)
+++ trunk/getfem/src/getfem_mesh_fem_global_function.cc Fri Jul 22 16:43:34 2016
@@ -24,20 +24,6 @@
 
 namespace getfem {
 
-  void mesh_fem_global_function::set_functions
-  (pglobal_function f, const mesh_im &mim) { 
-    std::vector<pglobal_function> funcs(1);
-    funcs[0]=f;
-    set_functions(funcs, mim);
-  }
-
-  void mesh_fem_global_function::set_functions
-  (pglobal_function f1, pglobal_function f2, const mesh_im &mim) { 
-    std::vector<pglobal_function> funcs(2);
-    funcs[0]=f1;
-    funcs[1] = f2;
-    set_functions(funcs, mim);
-  }
 
   void mesh_fem_global_function::set_functions
   (const std::vector<pglobal_function>& funcs, const mesh_im &mim) {