[Getfem-commits] (no subject)

[Liang Jin Lim]

branch: uniform_composite_fem
commit e2db190543a027ce641c31c94d23b1ec290735d2
Author: lj <address@hidden>
Date:   Mon Jun 4 16:52:27 2018 +0200

    Revert "Use single polynomial if mesh is uniform."
    
    This reverts commit f269855c03fab8439e2c36dc90a22bad69a0ed86.
---
 src/bgeot_poly_composite.cc       | 59 ++++++++++-----------------------------
 src/getfem/bgeot_poly_composite.h |  9 ++----
 2 files changed, 18 insertions(+), 50 deletions(-)

diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index 86c3d6e..a5d8418 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -70,18 +70,9 @@ namespace bgeot {
     for (size_type i = 0; i <= m.points().index().last_true(); ++i) {
       vertexes.add(m.points()[i]);
     }
-
-    bool is_uniform = true;
-    pgeometric_trans pgt_old = nullptr;
-
     for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
 
       pgeometric_trans pgt = m.trans_of_convex(cv);
-      if (is_uniform) {
-        if (!pgt_old) pgt_old = pgt;
-        else is_uniform = pgt_old == pgt;
-      }
-
       size_type N = pgt->structure()->dim();
       size_type P = m.dim();
       GMM_ASSERT1(pgt->is_linear() && N == P, "Bad geometric transformation");
@@ -100,7 +91,6 @@ namespace bgeot {
       gtrans[cv] = B0;
       orgs[cv] = m.points_of_convex(cv)[0];
     }
-    uniform = false;
   }
 
   scalar_type polynomial_composite::eval(const base_node &pt) const {
@@ -130,32 +120,30 @@ namespace bgeot {
             elt[ii] = false;
             p0 = pt; p0 -= mp->orgs[ii];
             gmm::mult(gmm::transposed(mp->gtrans[ii]), p0, p1);
-            auto &ref_poly = uniform ? poly : polytab[ii];
             if (mp->trans_of_convex(ii)->convex_ref()->is_in(p1) < 1E-10)
-              return local_coordinate ? to_scalar(ref_poly.eval(p1.begin()))
-              : to_scalar(ref_poly.eval(pt.begin()));
+              return local_coordinate ? to_scalar(polytab[ii].eval(p1.begin()))
+              : to_scalar(polytab[ii].eval(pt.begin()));
           }
         }
         ++it1; i1 = it1.index();
       }
 
-      if (i2 != size_type(-1))
+      if (i2 != size_type(-1)) 
       {
         const mesh_structure::ind_cv_ct &tc
           = mp->linked_mesh().convex_to_point(i2);
         itc = tc.begin(); itce = tc.end();
-        for (; itc != itce; ++itc)
+        for (; itc != itce; ++itc) 
         {
           size_type ii = *itc;
-          if (elt[ii])
+          if (elt[ii]) 
           {
             elt[ii] = false;
             p0 = pt; p0 -= mp->orgs[ii];
             gmm::mult(gmm::transposed(mp->gtrans[ii]), p0, p1);
-            auto &ref_poly = uniform ? poly : polytab[ii];
             if (mp->trans_of_convex(ii)->convex_ref()->is_in(p1) < 1E-10)
-              return  local_coordinate ? to_scalar(ref_poly.eval(p1.begin()))
-              : to_scalar(ref_poly.eval(pt.begin()));
+              return  local_coordinate ? 
to_scalar(polytab[ii].eval(p1.begin()))
+              : to_scalar(polytab[ii].eval(pt.begin()));
           }
         }
         --it2; i2 = it2.index();
@@ -167,9 +155,8 @@ namespace bgeot {
 
   polynomial_composite::polynomial_composite(const mesh_precomposite &m,
     bool lc)
-    : mp(&m), polytab(m.uniform ? 0 : m.nb_convex()),
-      local_coordinate(lc), poly(m.dim(), 0), uniform(m.uniform) {
-      if (!uniform) std::fill(polytab.begin(), polytab.end(), 
base_poly(m.dim(), 0));
+    : mp(&m), polytab(m.nb_convex()), local_coordinate(lc) {
+      std::fill(polytab.begin(), polytab.end(), base_poly(m.dim(), 0));
   }
 
   void polynomial_composite::derivative(short_type k) {
@@ -177,31 +164,15 @@ namespace bgeot {
       dim_type N = mp->dim();
       base_poly P(N, 0), Q;
       base_vector e(N), f(N);
-      if (uniform) {
+      for (size_type ic = 0; ic < mp->nb_convex(); ++ic) {
         gmm::clear(e); e[k] = 1.0;
-        gmm::mult(gmm::transposed(mp->gtrans[0]), e, f);
+        gmm::mult(gmm::transposed(mp->gtrans[ic]), e, f);
         P.clear();
         for (dim_type n = 0; n < N; ++n)
-        {
-          Q = poly;
-          Q.derivative(n);
-          P += Q * f[n];
-        }
-        poly = P;
-      }
-      else {
-        for (size_type ic = 0; ic < mp->nb_convex(); ++ic) {
-          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;
-        }
+        { Q = polytab[ic];
+        Q.derivative(n);
+        P += Q * f[n];  }
+        polytab[ic] = P;
       }
     }
     else
diff --git a/src/getfem/bgeot_poly_composite.h 
b/src/getfem/bgeot_poly_composite.h
index d6b9760..5c43f65 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -78,8 +78,7 @@ namespace bgeot {
     std::vector<base_matrix> gtrans;
     std::vector<scalar_type> det;
     std::vector<base_node> orgs;
-    bool uniform = false;
-
+    
     const basic_mesh &linked_mesh(void) const { return *msh; }
     size_type nb_convex(void) const { return gtrans.size(); }
     dim_type dim(void) const { return msh->dim(); }
@@ -97,18 +96,16 @@ namespace bgeot {
   protected :
     const mesh_precomposite *mp;
     std::vector<base_poly> polytab;
-    base_poly poly;
     bool local_coordinate; // are the polynomials described on the local
     // coordinates of each sub-element or on global coordinates.
-    bool uniform = true;
 
   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 &poly_of_subelt(size_type l) { return uniform ? poly : 
polytab[l]; }
-    const base_poly &poly_of_subelt(size_type l) const { return uniform ? poly 
: polytab[l]; }
+    base_poly &poly_of_subelt(size_type l) { return polytab[l]; }
+    const base_poly &poly_of_subelt(size_type l) const { return polytab[l]; }
     size_type nb_subelt() const { return polytab.size(); }
 
     polynomial_composite(bool lc = true) : local_coordinate(lc) {}