[Getfem-commits] (no subject)

[Liang Jin Lim]

branch: uniform_composite_fem
commit d21a14d99e2ffb394617b3e67452594c681b4669
Author: lj <address@hidden>
Date:   Mon Jun 11 08:43:05 2018 +0200

    Save memory usage by storing polynomial as static_stored_object and store 
its convex and key in a map.
---
 src/bgeot_poly_composite.cc       | 51 ++++++++++++++++++++++-----------------
 src/getfem/bgeot_poly.h           |  3 ++-
 src/getfem/bgeot_poly_composite.h | 15 ++++++------
 src/getfem_fem_composite.cc       |  8 +++---
 4 files changed, 42 insertions(+), 35 deletions(-)

diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index 7ce7569..5a2fb8c 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -120,11 +120,8 @@ namespace bgeot {
             elt[ii] = false;
             p0 = pt; p0 -= mp->orgs[ii];
             gmm::mult(gmm::transposed(mp->gtrans[ii]), p0, p1);
-            auto it = polytab.find(ii);
-            if (it == polytab.end()) return 0.0;
             if (mp->trans_of_convex(ii)->convex_ref()->is_in(p1) < 1E-10) {
-              return local_coordinate ? to_scalar(it->second.eval(p1.begin()))
-                                      : to_scalar(it->second.eval(pt.begin()));
+              return to_scalar(poly_of_subelt(ii).eval(local_coordinate ? 
p1.begin() : pt.begin()));
             }
           }
         }
@@ -144,11 +141,8 @@ namespace bgeot {
             elt[ii] = false;
             p0 = pt; p0 -= mp->orgs[ii];
             gmm::mult(gmm::transposed(mp->gtrans[ii]), p0, p1);
-            auto it = polytab.find(ii);
-            if (it == polytab.end()) return 0.0;
             if (mp->trans_of_convex(ii)->convex_ref()->is_in(p1) < 1E-10) {
-              return  local_coordinate ? to_scalar(it->second.eval(p1.begin()))
-                                       : 
to_scalar(it->second.eval(pt.begin()));
+              return to_scalar(poly_of_subelt(ii).eval(local_coordinate ? 
p1.begin() : pt.begin()));
             }
           }
         }
@@ -157,10 +151,12 @@ namespace bgeot {
     }
     GMM_ASSERT1(false, "Element not found in composite polynomial: " << pt);
   }
-
+  
+  DAL_TRIPLE_KEY(base_poly_key, short_type, short_type, 
std::vector<opt_long_scalar_type>);
 
   polynomial_composite::polynomial_composite(
-    const mesh_precomposite &m, bool lc) : mp(&m), local_coordinate(lc) {}
+    const mesh_precomposite &m, bool lc)
+    : mp(&m), local_coordinate(lc), default_poly(mp->dim(), 0) {}
 
   void polynomial_composite::derivative(short_type k) {
     if (local_coordinate) {
@@ -171,26 +167,37 @@ namespace bgeot {
         gmm::clear(e); e[k] = 1.0;
         gmm::mult(gmm::transposed(mp->gtrans[ic]), e, f);
         P.clear();
-        for (dim_type n = 0; n < N; ++n)
-        { Q = polytab[ic];
-        Q.derivative(n);
-        P += Q * f[n];  }
-        polytab[ic] = P;
+        auto &poly = poly_of_subelt(ic);
+        for (dim_type n = 0; n < N; ++n) {
+          Q = poly;
+          Q.derivative(n);
+          P += Q * f[n];
+        }
+        if (polytab.find(ic) != polytab.end()) set_poly_of_subelt(ic, P);
       }
     }
     else
-      for (size_type ic = 0; ic < mp->nb_convex(); ++ic)
-        polytab[ic].derivative(k);
+    for (size_type ic = 0; ic < mp->nb_convex(); ++ic) {
+      auto poly = poly_of_subelt(ic);
+      poly.derivative(k);
+      if (polytab.find(ic) != polytab.end()) set_poly_of_subelt(ic, poly);
+    }
   }
 
-  base_poly &polynomial_composite::set_poly_of_subelt(size_type l) {
-    return polytab[l];
+  void polynomial_composite::set_poly_of_subelt(size_type l, const base_poly 
&poly) {
+    auto poly_key = std::make_shared<base_poly_key>(poly.degree(), poly.dim(), 
poly);
+    auto o = dal::search_stored_object(poly_key);
+    if (!o) {
+      o = std::make_shared<base_poly>(poly);
+      dal::add_stored_object(poly_key, o);
+    }
+    polytab[l] = poly_key;
   }
 
-  const base_poly *polynomial_composite::poly_of_subelt(size_type l) const {
+  const base_poly &polynomial_composite::poly_of_subelt(size_type l) const {
     auto it = polytab.find(l);
-    if (it == polytab.end()) return nullptr;
-    else return &it->second;
+    if (it == polytab.end()) return default_poly;
+    return dynamic_cast<const base_poly 
&>(*dal::search_stored_object(it->second));
   }
 
 
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index 67e0502..edcd7dd 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -40,6 +40,7 @@
 */
 
 #include "bgeot_config.h"
+#include "dal_static_stored_objects.h";
 #include <vector>
 
 namespace bgeot
@@ -176,7 +177,7 @@ namespace bgeot
    *        The resulting polynomials have a smaller degree.
    *
    */
-  template<typename T> class polynomial : public std::vector<T> {
+  template<typename T> class polynomial : public std::vector<T>, public 
dal::static_stored_object {
   protected :
     
     short_type n, d;
diff --git a/src/getfem/bgeot_poly_composite.h 
b/src/getfem/bgeot_poly_composite.h
index 2a83a70..80df3f4 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -95,17 +95,18 @@ namespace bgeot {
 
   protected :
     const mesh_precomposite *mp;
-    std::map<size_type, base_poly> polytab;
-    bool local_coordinate; // are the polynomials described on the local
-    // coordinates of each sub-element or on global coordinates.
+    std::map<size_type, dal::pstatic_stored_object_key> polytab;
+    bool local_coordinate;  // are the polynomials described on the local
+                            // coordinates of each sub-element or on global 
coordinates.
+    base_poly default_poly;
 
   public :
     
     template <class ITER> scalar_type eval(const ITER &it) const;
     scalar_type eval(const base_node &pt) const;
     void derivative(short_type k);
-    base_poly &set_poly_of_subelt(size_type l);
-    const base_poly *poly_of_subelt(size_type l) const;
+    void set_poly_of_subelt(size_type l, const base_poly &poly);
+    const base_poly &poly_of_subelt(size_type l) const;
     size_type nb_subelt() const { return polytab.size(); }
 
     polynomial_composite(bool lc = true) : local_coordinate(lc) {}
@@ -117,9 +118,7 @@ namespace bgeot {
   (std::ostream &o, const polynomial_composite& P) {
     o << "poly_composite [";
     for (size_type i = 0; i < P.nb_subelt(); ++i) {
-      if (i == 0) continue;
-      auto poly = P.poly_of_subelt(i);
-      o << ", ";  o << (poly) ? *poly : base_poly();
+      if (i != 0) o << ", " << P.poly_of_subelt(i);
     }
     o << "]";
     return o;
diff --git a/src/getfem_fem_composite.cc b/src/getfem_fem_composite.cc
index 90d9ed4..2241b98 100644
--- a/src/getfem_fem_composite.cc
+++ b/src/getfem_fem_composite.cc
@@ -59,7 +59,7 @@ namespace getfem {
       for (size_type k = 0; k < pf->nb_dof(cv); ++k) {
        size_type igl = mf.ind_basic_dof_of_element(cv)[k];
        base_poly fu = pf->base()[k];
-       base[igl].set_poly_of_subelt(cv) = fu;
+       base[igl].set_poly_of_subelt(cv, fu);
        dofd[igl] = pf->dof_types()[k];
       }
     }
@@ -189,7 +189,7 @@ namespace getfem {
       (4, bgeot::polynomial_composite(mp, false));
     for (size_type k = 0; k < 4; ++k)
       for (size_type ic = 0; ic < 3; ++ic)
-        base()[k].set_poly_of_subelt(ic) = bgeot::read_base_poly(2, s);
+        base()[k].set_poly_of_subelt(ic, bgeot::read_base_poly(2, s));
 
     for (size_type i = 0; i < 3; ++i) {
       base_node pt(0.0, 0.0);
@@ -355,7 +355,7 @@ namespace getfem {
       (12, bgeot::polynomial_composite(mp, false));
     for (size_type k = 0; k < 12; ++k)
       for (size_type ic = 0; ic < 3; ++ic)
-        base()[k].set_poly_of_subelt(ic) = bgeot::read_base_poly(2, s);
+        base()[k].set_poly_of_subelt(ic, bgeot::read_base_poly(2, s));
 
     for (size_type i = 0; i < 3; ++i) {
       base_node pt(0.0, 0.0);
@@ -621,7 +621,7 @@ namespace getfem {
       (16, bgeot::polynomial_composite(mp, false));
     for (size_type k = 0; k < 16; ++k)
       for (size_type ic = 0; ic < 4; ++ic)
-        base()[k].set_poly_of_subelt(ic) = bgeot::read_base_poly(2, s);
+        base()[k].set_poly_of_subelt(ic, bgeot::read_base_poly(2, s));
 
     for (size_type i = 0; i < 4; ++i) {
       base_node pt(0.0, 0.0);