branch: devel-logari81
commit 267e7fa26603c3a0dd32fcd282c1a4105481c008
Author: Konstantinos Poulios <logar...@gmail.com>
Date:   Sun Aug 6 21:52:55 2017 +0200

    white space and minor coding style normalization
---
 src/bgeot_convex_ref.cc                   |  21 +-
 src/bgeot_convex_structure.cc             |  51 +-
 src/bgeot_geometric_trans.cc              |  12 +-
 src/bgeot_poly.cc                         |  56 +-
 src/getfem/bgeot_convex_ref.h             |  16 +-
 src/getfem/bgeot_convex_structure.h       |   2 +-
 src/getfem/bgeot_geometric_trans.h        |  12 +-
 src/getfem/bgeot_kdtree.h                 |  29 +-
 src/getfem_contact_and_friction_common.cc |   8 +-
 src/getfem_export.cc                      |  20 +-
 src/getfem_fem.cc                         |  80 +--
 src/getfem_import.cc                      |  32 +-
 src/getfem_integration.cc                 | 869 +++++++++++++++---------------
 src/getfem_mesh.cc                        |   4 +-
 src/getfem_mesh_fem.cc                    | 118 ++--
 15 files changed, 677 insertions(+), 653 deletions(-)

diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc
index 53c779f..2b3abe4 100644
--- a/src/bgeot_convex_ref.cc
+++ b/src/bgeot_convex_ref.cc
@@ -444,15 +444,18 @@ namespace bgeot {
     dal::pstatic_stored_object_key
       pk = std::make_shared<product_ref_key_>(a, b);
     dal::pstatic_stored_object o = dal::search_stored_object(pk);
-    if (o) return std::dynamic_pointer_cast<const convex_of_reference>(o);
-    pconvex_ref p = std::make_shared<product_ref_>(a, b);
-    dal::add_stored_object(pk, p, a, b,
-                           convex_product_structure(a->structure(),
-                                                    b->structure()),
-                           p->pspt(), dal::PERMANENT_STATIC_OBJECT);
-    pconvex_ref p1 = basic_convex_ref(p);
-    if (p != p1) add_dependency(p, p1);
-    return p;
+    if (o)
+      return std::dynamic_pointer_cast<const convex_of_reference>(o);
+    else {
+      pconvex_ref p = std::make_shared<product_ref_>(a, b);
+      dal::add_stored_object(pk, p, a, b,
+                             convex_product_structure(a->structure(),
+                                                      b->structure()),
+                             p->pspt(), dal::PERMANENT_STATIC_OBJECT);
+      pconvex_ref p1 = basic_convex_ref(p);
+      if (p != p1) add_dependency(p, p1);
+      return p;
+    }
   }
 
   pconvex_ref parallelepiped_of_reference(dim_type nc, dim_type k) {
diff --git a/src/bgeot_convex_structure.cc b/src/bgeot_convex_structure.cc
index bb936f3..fe14175 100644
--- a/src/bgeot_convex_structure.cc
+++ b/src/bgeot_convex_structure.cc
@@ -299,7 +299,7 @@ namespace bgeot {
             : convex_product_structure(cv1->faces_structure()[i], cv2);
 
         faces[i] = std::vector<short_type>(cv1->nb_points_of_face(i)
-                                              * cv2->nb_points());
+                                           * cv2->nb_points());
 
         for (short_type j = 0; j < cv1->nb_points_of_face(i); j++)
           for (short_type l = 0; l < cv2->nb_points(); l++) {
@@ -361,18 +361,23 @@ namespace bgeot {
   DAL_DOUBLE_KEY(parallelepiped_key_, dim_type, dim_type);
 
   pconvex_structure parallelepiped_structure(dim_type nc, dim_type k) {
-    if (nc <= 1) return simplex_structure(nc, k);
+    if (nc <= 1)
+      return simplex_structure(nc, k);
+
     dal::pstatic_stored_object_key
       pcsk = std::make_shared<parallelepiped_key_>(nc, k);
 
     dal::pstatic_stored_object o = dal::search_stored_object(pcsk);
-    if (o) return ((std::dynamic_pointer_cast<const parallelepiped_>(o))->p);
-    auto p = std::make_shared<parallelepiped_>();
-    p->p = convex_product_structure(parallelepiped_structure(dim_type(nc-1),k),
-                                    simplex_structure(1,k));
-    dal::add_stored_object(pcsk, dal::pstatic_stored_object(p), p->p,
-                           dal::PERMANENT_STATIC_OBJECT);
-    return p->p;
+    if (o)
+      return ((std::dynamic_pointer_cast<const parallelepiped_>(o))->p);
+    else {
+      auto p = std::make_shared<parallelepiped_>();
+      p->p = 
convex_product_structure(parallelepiped_structure(dim_type(nc-1),k),
+                                      simplex_structure(1,k));
+      dal::add_stored_object(pcsk, dal::pstatic_stored_object(p), p->p,
+                             dal::PERMANENT_STATIC_OBJECT);
+      return p->p;
+    }
   }
 
   /* ******************************************************************** */
@@ -398,7 +403,7 @@ namespace bgeot {
     p->Nc = nc;
     p->nbpt = (nc == 2) ? 8 : 20;
     p->nbf =  (nc == 2) ? 4 : 6;
-    p->basic_pcvs =  parallelepiped_structure(nc);
+    p->basic_pcvs = parallelepiped_structure(nc); // k=1
     p->faces_struct.resize(p->nbf);
     p->faces = std::vector< std::vector<short_type> >(p->nbf);
     p->dir_points_ = std::vector<short_type>(p->Nc + 1);
@@ -468,24 +473,24 @@ namespace bgeot {
 
   pconvex_structure pyramidal_structure(dim_type k) {
     GMM_ASSERT1(k == 1 || k == 2, "Sorry, pyramidal elements implemented "
-               "only for degree one or two.");
+                "only for degree one or two.");
     dal::pstatic_stored_object_key
       pcsk = std::make_shared<pyramidal_structure_key_>(k);
     dal::pstatic_stored_object o = dal::search_stored_object(pcsk);
-    if (o) return std::dynamic_pointer_cast<const convex_structure>(o);
+    if (o)
+      return std::dynamic_pointer_cast<const convex_structure>(o);
 
     auto p = std::make_shared<pyramidal_structure_>();
     pconvex_structure pcvs(p);
-    
+
     p->Nc = 3;
     p->dir_points_ = std::vector<short_type>(p->Nc + 1);
-    
 
     if (k == 1) {
       p->nbpt = 5;
       p->nbf = 5;
       p->auto_basic = true;
-      //    4 
+      //    4
       //   /|||
       //  / || |
       // 2-|--|-3
@@ -508,12 +513,12 @@ namespace bgeot {
 
       p->faces_struct[0] = parallelepiped_structure(2);
       for (int i = 1; i < p->nbf; i++)
-       p->faces_struct[i] = simplex_structure(2);
+        p->faces_struct[i] = simplex_structure(2);
 
       dal::add_stored_object(pcsk, pcvs, parallelepiped_structure(2),
-                            simplex_structure(2),
-                            dal::PERMANENT_STATIC_OBJECT);
-      
+                             simplex_structure(2),
+                             dal::PERMANENT_STATIC_OBJECT);
+
     } else {
       p->nbpt = 14;
       p->nbf = 5;
@@ -544,15 +549,15 @@ namespace bgeot {
 
       p->faces_struct[0] = parallelepiped_structure(2, 2);
       for (int i = 1; i < p->nbf; i++)
-       p->faces_struct[i] = simplex_structure(2, 2);
+        p->faces_struct[i] = simplex_structure(2, 2);
 
       dal::add_stored_object(pcsk, pcvs, parallelepiped_structure(2, 2),
-                            simplex_structure(2, 2),
-                            dal::PERMANENT_STATIC_OBJECT);
+                             simplex_structure(2, 2),
+                             dal::PERMANENT_STATIC_OBJECT);
     }
     return pcvs;
   }
-  
+
   /* ******************************************************************** */
   /*        Generic dummy convex with n global nodes.                     */
   /* ******************************************************************** */
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 1d611e5..d9aad12 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -472,7 +472,7 @@ namespace bgeot {
           { vertex = false; break; }
       if (vertex) vertices_.push_back(ip);
     }
-    auto dimension = dim();
+    dim_type dimension = dim();
     if (dynamic_cast<const torus_geom_trans *>(this)) dimension -= 1;
     GMM_ASSERT1(vertices_.size() > dimension, "Internal error");
   }
@@ -801,7 +801,7 @@ namespace bgeot {
             "2*x*x*y + 2*x*y*y - 3*x*y;");
 
         for (int i = 0; i < 8; ++i)
-          trans[i] = bgeot::read_base_poly(2, s);
+          trans[i] = read_base_poly(2, s);
       } else {
         std::stringstream s
           ("1 + 2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2"
@@ -831,7 +831,7 @@ namespace bgeot {
            "2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2 - 5*x*y*z;");
 
         for (int i = 0; i < 20; ++i)
-          trans[i] = bgeot::read_base_poly(3, s);
+          trans[i] = read_base_poly(3, s);
       }
       fill_standard_vertices();
     }
@@ -907,14 +907,14 @@ namespace bgeot {
   };
 
   static pgeometric_trans
-    pyramidal_gt(gt_param_list& params,
-                     std::vector<dal::pstatic_stored_object> &dependencies) {
+  pyramidal_gt(gt_param_list& params,
+               std::vector<dal::pstatic_stored_object> &deps) {
     GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
                 << params.size() << " should be 1.");
     GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
     int k = int(::floor(params[0].num() + 0.01));
 
-    dependencies.push_back(pyramidal_element_of_reference(dim_type(k)));
+    deps.push_back(pyramidal_element_of_reference(dim_type(k)));
     return std::make_shared<pyramidal_trans_>(dim_type(k));
   }
 
diff --git a/src/bgeot_poly.cc b/src/bgeot_poly.cc
index de5ba06..545125e 100644
--- a/src/bgeot_poly.cc
+++ b/src/bgeot_poly.cc
@@ -33,9 +33,9 @@ namespace bgeot {
     static bool init = false;
     if (!init) {
       for (short_type i = 0; i < STORED; ++i) {
-       alpha_M_(i, 0) = alpha_M_(0, i) = 1;
-       for (short_type j = 1; j <= i; ++j)
-         alpha_M_(i,j) = alpha_M_(j,i) = (alpha_M_(i, j-1) * (i+j)) / j;
+        alpha_M_(i, 0) = alpha_M_(0, i) = 1;
+        for (short_type j = 1; j <= i; ++j)
+          alpha_M_(i,j) = alpha_M_(j,i) = (alpha_M_(i, j-1) * (i+j)) / j;
       }
       init = true;
     }
@@ -46,7 +46,7 @@ namespace bgeot {
   size_type alpha(short_type n, short_type d) {
     alpha_init_();
     GMM_ASSERT1(n < STORED && d < STORED,
-               "alpha called with n = " << n << " and d = " << d);
+                "alpha called with n = " << n << " and d = " << d);
     return alpha_(n,d);
   }
 
@@ -56,7 +56,7 @@ namespace bgeot {
       size_type g_idx = global_index_; short_type deg = degree_;
       reverse_iterator it = rbegin() + 1;
       for (l = short_type(n-2); l != short_type(-1); --l, ++it)
-       if (*it != 0) break;
+        if (*it != 0) break;
       short_type a = (*this)[n-1]; (*this)[n-1] = 0;
       (*this)[short_type(l+1)] = short_type(a + 1);
       if (l != short_type(-1)) ((*this)[l])--;
@@ -73,12 +73,12 @@ namespace bgeot {
       size_type g_idx = global_index_; short_type deg = degree_;
       reverse_iterator it = rbegin();
       for (l = short_type(n-1); l != short_type(-1); --l, ++it) {
-       if (*it != 0) break;
+        if (*it != 0) break;
       }
       if (l != short_type(-1)) {
-       short_type a = (*this)[l];
-       (*this)[l] = 0; (*this)[n-1] = short_type(a - 1);
-       if (l > 0) ((*this)[l-1])++; 
+        short_type a = (*this)[l];
+        (*this)[l] = 0; (*this)[n-1] = short_type(a - 1);
+        if (l > 0) ((*this)[l-1])++; 
         else if (short_type(deg+1)) degree_ = short_type(deg-1);
       }
       if (g_idx+1) global_index_ = g_idx-1;
@@ -144,20 +144,20 @@ namespace bgeot {
       else if (s == "u" && n > 5) result = base_poly(n, 1, 5);
       else if (s == "t" && n > 6) result = base_poly(n, 1, 6);
       else if (s == "sqrt") {
-       base_poly p = read_expression(n, f);
-       if (p.degree() > 0) parse_error(1);
-       result.one();  result *= sqrt(p[0]);
+        base_poly p = read_expression(n, f);
+        if (p.degree() > 0) parse_error(1);
+        result.one();  result *= sqrt(p[0]);
       }
       else { parse_error(2); }
       break;
     case 5 :
       switch (s[0]) {
       case '(' :
-       result = read_base_poly(n, f);
-       j = get_next_token(s, f);
-       if (j != 5 || s[0] != ')') parse_error(3);
-       break;
-       default : parse_error(4);
+        result = read_base_poly(n, f);
+        j = get_next_token(s, f);
+        if (j != 5 || s[0] != ')') parse_error(3);
+        break;
+        default : parse_error(4);
       }
       break;
     default : parse_error(5);
@@ -178,8 +178,8 @@ namespace bgeot {
   }
 
   void do_bin_op(std::vector<base_poly> &value_list,
-                std::vector<int> &op_list,
-                std::vector<int> &prior_list) {
+                 std::vector<int> &op_list,
+                 std::vector<int> &prior_list) {
     base_poly &p2 = *(value_list.rbegin());
     if (op_list.back() != 6) {
       assert(value_list.size()>1);
@@ -190,14 +190,14 @@ namespace bgeot {
         case 3  : p1 += p2; break;
         case 4  : p1 -= p2; break;
         case 5  : 
-         {
-           if (p2.degree() > 0) parse_error(7);
-           int pow = int(to_scalar(p2[0]));
-           if (p2[0] !=  opt_long_scalar_type(pow) || pow < 0) parse_error(8);
-           base_poly p = p1; p1.one();
-           for (int i = 0; i < pow; ++i) p1 *= p;
-         }
-         break;
+          {
+            if (p2.degree() > 0) parse_error(7);
+            int pow = int(to_scalar(p2[0]));
+            if (p2[0] !=  opt_long_scalar_type(pow) || pow < 0) parse_error(8);
+            base_poly p = p1; p1.one();
+            for (int i = 0; i < pow; ++i) p1 *= p;
+          }
+          break;
         default: assert(0);
       }
       value_list.pop_back(); 
@@ -223,7 +223,7 @@ namespace bgeot {
     operator_priority_(i, i ? s[0] : '0', prior, op);
     while (op) {
       while (!prior_list.empty() && prior_list.back() <= prior)
-       do_bin_op(value_list, op_list, prior_list);
+        do_bin_op(value_list, op_list, prior_list);
 
       value_list.push_back(read_expression(n, f));
       op_list.push_back(op);
diff --git a/src/getfem/bgeot_convex_ref.h b/src/getfem/bgeot_convex_ref.h
index 8834afe..66db71e 100644
--- a/src/getfem/bgeot_convex_ref.h
+++ b/src/getfem/bgeot_convex_ref.h
@@ -111,12 +111,12 @@ namespace bgeot {
      */
     virtual scalar_type is_in_face(short_type, const base_node &) const =0;
     /// return the normal vector for each face.
-    const std::vector<base_small_vector> &normals(void) const
+    const std::vector<base_small_vector> &normals() const
     { return normals_; }
     /// return the vertices of the reference convex.
-    const stored_point_tab &points(void) const { return *ppoints; }
-    const stored_point_tab &points(void) { return *ppoints; }
-    pstored_point_tab pspt(void) const { return ppoints; }
+    const stored_point_tab &points() const { return *ppoints; }
+    const stored_point_tab &points() { return *ppoints; }
+    pstored_point_tab pspt() const { return ppoints; }
     virtual ~convex_of_reference()
     { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Convex of reference"); }
 
@@ -130,7 +130,7 @@ namespace bgeot {
 
   /// return the associated order 1 reference convex.
   inline pconvex_ref basic_convex_ref(pconvex_ref cvr)
-  { if (cvr->auto_basic) return cvr; else return cvr->basic_convex_ref_; }
+  { return cvr->auto_basic ? cvr : cvr->basic_convex_ref_; }
 
 
   //@name public functions for obtaining a convex of reference
@@ -147,13 +147,11 @@ namespace bgeot {
   */
   pconvex_ref Q2_incomplete_reference(dim_type d);
   /** tensorial product of two convex ref.
-      in order to ensure unicity, it is required the a->dim() >= b->dim()
-  */
+      in order to ensure unicity, it is required the a->dim() >= b->dim() */
   /** Pyramidal element of reference of degree k (k = 1 or 2 only) */
   pconvex_ref pyramidal_element_of_reference(dim_type k);
   pconvex_ref convex_ref_product(pconvex_ref a, pconvex_ref b);
-  /** equilateral simplex (degree 1). used only for mesh quality estimations
-   */
+  /** equilateral simplex (degree 1). used only for mesh quality estimations */
   pconvex_ref equilateral_simplex_of_reference(dim_type nc);
 
   /** generic convex with n global nodes      */
diff --git a/src/getfem/bgeot_convex_structure.h 
b/src/getfem/bgeot_convex_structure.h
index 0dcb509..9a8cbe3 100644
--- a/src/getfem/bgeot_convex_structure.h
+++ b/src/getfem/bgeot_convex_structure.h
@@ -118,7 +118,7 @@ namespace bgeot {
     /** Give a pointer array on the structures of the faces.
      *   faces_structure()[i] is a pointer on the structure of the face i.
      */
-    inline const convex_structure_faces_ct &faces_structure(void) const
+    inline const convex_structure_faces_ct &faces_structure() const
     { return faces_struct; }
     /** Return "direct" points indexes for a given face.
      *  @param i the face number.
diff --git a/src/getfem/bgeot_geometric_trans.h 
b/src/getfem/bgeot_geometric_trans.h
index 01d1deb..9851024 100644
--- a/src/getfem/bgeot_geometric_trans.h
+++ b/src/getfem/bgeot_geometric_trans.h
@@ -126,7 +126,7 @@ namespace bgeot {
     pconvex_structure structure() const { return cvr->structure(); }
     /// Basic structure of the reference element.
     pconvex_structure basic_structure() const
-    { return bgeot::basic_structure(cvr->structure()); }
+      { return bgeot::basic_structure(cvr->structure()); }
     /// Gives the value of the functions vector at a certain point.
     virtual void poly_vector_val(const base_node &pt, base_vector &val) const 
= 0;
     /// Gives the value of a subgroup of the functions vector at a certain 
point.
@@ -163,7 +163,7 @@ namespace bgeot {
         if the transformation is linear, x is not used at all */
     size_type complexity() const { return complexity_; }
     virtual ~geometric_trans()
-    { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Geometric transformation"); }
+      { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Geometric transformation"); }
     geometric_trans()
       { DAL_STORED_OBJECT_DEBUG_CREATED(this, "Geometric transformation"); }
   };
@@ -251,14 +251,14 @@ namespace bgeot {
    *
    *  pt is the position of the evaluation point on the reference element
    */
-  base_small_vector APIDECL compute_normal(const 
geotrans_interpolation_context& c,
-                                   size_type face);
+  base_small_vector APIDECL
+  compute_normal(const geotrans_interpolation_context& c, size_type face);
 
   /** return the local basis (i.e. the normal in the first column, and the
    *  tangent vectors in the other columns
    */
-  base_matrix
-  APIDECL compute_local_basis(const geotrans_interpolation_context& c,
+  base_matrix APIDECL
+  compute_local_basis(const geotrans_interpolation_context& c,
                       size_type face);
     //@}
 
diff --git a/src/getfem/bgeot_kdtree.h b/src/getfem/bgeot_kdtree.h
index 9d45555..81c00e9 100644
--- a/src/getfem/bgeot_kdtree.h
+++ b/src/getfem/bgeot_kdtree.h
@@ -52,7 +52,7 @@ namespace bgeot {
     virtual ~kdtree_elt_base() {}
   };
 
-  /// store a point and the associated index for the kdtree. 
+  /// store a point and the associated index for the kdtree.
   /* std::pair<size_type,base_node> is not ok since it does not
      have a suitable overloaded swap function ...
   */
@@ -69,15 +69,15 @@ namespace bgeot {
   typedef std::vector<index_node_pair> kdtree_tab_type;
 
   /** Balanced tree over a set of points.
-   
+
   Once the tree have been built, it is possible to query very
   quickly for the list of points lying in a given box. Note that
   this is not a dynamic structure: once you start to call
   kdtree::points_in_box, you should not use anymore kdtree::add_point.
-  
+
   Here is an example of use (which tries to find the mapping between
   the dof of the mesh_fem and the node numbers of its mesh):
-  
+
   @code
   void dof_to_nodes(const getfem::mesh_fem &mf) {
     const getfem::getfem_mesh &m = mf.linked_mesh();
@@ -94,9 +94,9 @@ namespace bgeot {
       tree.points_in_box(t, P, P2);
       if (t.size()) {
         assert(t.size() == 1);
-       cout << " dof " << d << " maps to mesh node " << t[0].i << "\n";
+        cout << " dof " << d << " maps to mesh node " << t[0].i << "\n";
       }
-    } 
+    }
   }
   @endcode
   */
@@ -110,30 +110,31 @@ namespace bgeot {
     void clear() { clear_tree(); pts = kdtree_tab_type(); N = 0; }
     void reserve(size_type n) { pts.reserve(n); }
     /// insert a new point
-    size_type add_point(const base_node& n) { 
+    size_type add_point(const base_node& n) {
       size_type i = pts.size(); add_point_with_id(n,i); return i;
     }
     /// insert a new point, with an associated number.
     void add_point_with_id(const base_node& n, size_type i) {
-      if (pts.size() == 0) N = n.size(); 
-      else if (N != n.size()) 
-       GMM_ASSERT2(false, "invalid dimension");
+      if (pts.size() == 0)
+        N = n.size();
+      else
+        GMM_ASSERT2(N == n.size(), "invalid dimension");
       if (tree) clear_tree();
       pts.push_back(index_node_pair(i, n));
     }
     size_type nb_points() const { return pts.size(); }
     const kdtree_tab_type &points() const { return pts; }
-    /* fills ipts with the indexes of points in the box 
+    /* fills ipts with the indexes of points in the box
        [min,max] */
     void points_in_box(kdtree_tab_type &ipts,
-                      const base_node &min, 
-                      const base_node &max);
+                       const base_node &min,
+                       const base_node &max);
     /* assigns at ipt the index of the nearest neighbor at location
        pos and returns the square of the distance to this point*/
     scalar_type nearest_neighbor(index_node_pair &ipt,
                                  const base_node &pos);
   private:
-    typedef std::vector<size_type>::const_iterator ITER;  
+    typedef std::vector<size_type>::const_iterator ITER;
     void clear_tree();
   };
 }
diff --git a/src/getfem_contact_and_friction_common.cc 
b/src/getfem_contact_and_friction_common.cc
index d65366d..385ce8c 100644
--- a/src/getfem_contact_and_friction_common.cc
+++ b/src/getfem_contact_and_friction_common.cc
@@ -1104,7 +1104,7 @@ namespace getfem {
 
             for (size_type subiter(0);;) {
               pps(a, hessa);
-             det = gmm::abs(bgeot::lu_inverse(&(*(hessa.begin())),N-1, false));
+              det = gmm::abs(bgeot::lu_inverse(&(*(hessa.begin())),N-1, 
false));
               if (det > 1E-15) break;
               for (size_type i = 0; i < N-1; ++i)
                 a[i] += gmm::random() * 1E-7;
@@ -1164,7 +1164,7 @@ namespace getfem {
 
             for (size_type subiter(0);;) {
               pps(a, hessa);
-             det = gmm::abs(bgeot::lu_inverse(&(*(hessa.begin())),N-1, false));
+              det = gmm::abs(bgeot::lu_inverse(&(*(hessa.begin())),N-1, 
false));
               if (det > 1E-15) break;
               for (size_type i = 0; i < N-1; ++i)
                 a[i] += gmm::random() * 1E-7;
@@ -2677,8 +2677,8 @@ namespace getfem {
       if (args.size() != 6) return false;
       size_type N =  args[0]->size();
       if (N == 0 || args[1]->size() != N || args[2]->size() != N
-         || args[3]->size() != 1 || args[4]->size() > 3
-         || args[4]->size() == 0 || args[5]->size() != 1 ) return false;
+          || args[3]->size() != 1 || args[4]->size() > 3
+          || args[4]->size() == 0 || args[5]->size() != 1 ) return false;
       ga_init_vector(sizes, N);
       return true;
     }
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index 14da099..7edd6cf 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -168,11 +168,11 @@ namespace getfem
         pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 
1)
                                            : classical_fem(pgt, 1);
 
-       short_type degree = 1;
-       if ((pf != classical_pf1 && pf->estimated_degree() > 1) ||
-               pgt->structure() != pgt->basic_structure())
-         degree = 2;
-         
+        short_type degree = 1;
+        if ((pf != classical_pf1 && pf->estimated_degree() > 1) ||
+            pgt->structure() != pgt->basic_structure())
+          degree = 2;
+
         pmf->set_finite_element(cv, discontinuous ?
                                 classical_discontinuous_fem(pgt, degree) :
                                 classical_fem(pgt, degree));
@@ -480,9 +480,9 @@ namespace getfem
     for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
       bgeot::pgeometric_trans pgt2 = mf.linked_mesh().trans_of_convex(cv);
       GMM_ASSERT1(basic_structure(pgt->structure()) ==
-                 basic_structure(pgt2->structure()),
-                 "Cannot export this mesh to opendx, it contains "
-                 "different convex types. Slice it first.");
+                  basic_structure(pgt2->structure()),
+                  "Cannot export this mesh to opendx, it contains "
+                  "different convex types. Slice it first.");
       pfem pf = mf.fem_of_element(cv);
       bool discontinuous = false;
       for (unsigned i=0; i < pf->nb_dof(cv); ++i) {
@@ -875,7 +875,7 @@ namespace getfem
         if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
       }
       pfem classical_pf1 = discontinuous ?
-       classical_discontinuous_fem(pgt, 1) : classical_fem(pgt, 1);
+        classical_discontinuous_fem(pgt, 1) : classical_fem(pgt, 1);
       pmf->set_finite_element(cv, classical_pf1);
     }
     psl = NULL;
@@ -890,7 +890,7 @@ namespace getfem
         case 4:
           if ( 2 == pmf->fem_of_element(cv)->dim() ) t = POS_QU;
           else if (3 == pmf->fem_of_element(cv)->dim()) t = POS_SI;
-         break;
+          break;
         case 6: t = POS_PR; break;
         case 8: t = POS_HE; break;
         case 5: t = POS_PY; break;
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index 96a1b8f..5d81180 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -36,6 +36,8 @@
 
 namespace getfem {
 
+  using bgeot::read_base_poly;
+  using bgeot::base_rational_fraction;
 
   const base_matrix& fem_interpolation_context::M() const {
     if (gmm::mat_nrows(M_) == 0) {
@@ -1157,7 +1159,7 @@ namespace getfem {
           "4*(x*y - x*x*y);"
           "2*x*x*y + 2*x*y*y - 3*x*y;");
 
-      for (int i = 0; i < 8; ++i) p->base()[i] = bgeot::read_base_poly(2, s);
+      for (int i = 0; i < 8; ++i) p->base()[i] = read_base_poly(2, s);
 
       p->add_node(lagrange_dof(2), base_small_vector(0.0, 0.0));
       p->add_node(lagrange_dof(2), base_small_vector(0.5, 0.0));
@@ -1195,7 +1197,7 @@ namespace getfem {
          "4*( - x^2*y*z + x*y*z);"
          "2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2 - 5*x*y*z;");
 
-      for (int i = 0; i < 20; ++i) p->base()[i] = bgeot::read_base_poly(3, s);
+      for (int i = 0; i < 20; ++i) p->base()[i] = read_base_poly(3, s);
 
       p->add_node(lagrange_dof(3), base_small_vector(0.0, 0.0, 0.0));
       p->add_node(lagrange_dof(3), base_small_vector(0.5, 0.0, 0.0));
@@ -1256,7 +1258,7 @@ namespace getfem {
   // 0---1---2
 
   static pfem build_pyramidal_pk_fem(short_type k, bool disc) {
-    auto p = std::make_shared<fem<bgeot::base_rational_fraction>>();
+    auto p = std::make_shared<fem<base_rational_fraction>>();
     p->mref_convex() = bgeot::pyramidal_element_of_reference(1);
     p->dim() = 3;
     p->is_standard() = p->is_equivalent() = true;
@@ -1268,17 +1270,17 @@ namespace getfem {
 
     if (k == 0) {
       p->base().resize(1);
-      p->base()[0] = bgeot::read_base_poly(3, "1");
+      p->base()[0] = read_base_poly(3, "1");
       p->add_node(lagrange_0_dof(3), base_small_vector(0.0, 0.0, 0.5));
     } else if (k == 1) {
       p->base().resize(5);
-      bgeot::base_rational_fraction // Q = xy/(1-z)
-        Q(bgeot::read_base_poly(3, "x*y"), bgeot::read_base_poly(3, "1-z"));
-      p->base()[0] = (bgeot::read_base_poly(3, "1-x-y-z") + Q)*0.25;
-      p->base()[1] = (bgeot::read_base_poly(3, "1+x-y-z") - Q)*0.25;
-      p->base()[2] = (bgeot::read_base_poly(3, "1-x+y-z") - Q)*0.25;
-      p->base()[3] = (bgeot::read_base_poly(3, "1+x+y-z") + Q)*0.25;
-      p->base()[4] = bgeot::read_base_poly(3, "z");
+      base_rational_fraction // Q = xy/(1-z)
+        Q(read_base_poly(3, "x*y"), read_base_poly(3, "1-z"));
+      p->base()[0] = (read_base_poly(3, "1-x-y-z") + Q)*0.25;
+      p->base()[1] = (read_base_poly(3, "1+x-y-z") - Q)*0.25;
+      p->base()[2] = (read_base_poly(3, "1-x+y-z") - Q)*0.25;
+      p->base()[3] = (read_base_poly(3, "1+x+y-z") + Q)*0.25;
+      p->base()[4] = read_base_poly(3, "z");
 
       p->add_node(lag_dof, base_small_vector(-1.0, -1.0, 0.0));
       p->add_node(lag_dof, base_small_vector( 1.0, -1.0, 0.0));
@@ -1289,16 +1291,16 @@ namespace getfem {
     } else if (k == 2) {
       p->base().resize(14);
 
-      base_poly xi0  = bgeot::read_base_poly(3, "(1-z-x)*0.5");
-      base_poly xi1  = bgeot::read_base_poly(3, "(1-z-y)*0.5");
-      base_poly xi2  = bgeot::read_base_poly(3, "(1-z+x)*0.5");
-      base_poly xi3  = bgeot::read_base_poly(3, "(1-z+y)*0.5");
-      base_poly x    = bgeot::read_base_poly(3, "x");
-      base_poly y    = bgeot::read_base_poly(3, "y");
-      base_poly z    = bgeot::read_base_poly(3, "z");
-      base_poly ones = bgeot::read_base_poly(3, "1");
-      base_poly un_z = bgeot::read_base_poly(3, "1-z");
-      bgeot::base_rational_fraction Q(bgeot::read_base_poly(3, "1"), un_z);
+      base_poly xi0  = read_base_poly(3, "(1-z-x)*0.5");
+      base_poly xi1  = read_base_poly(3, "(1-z-y)*0.5");
+      base_poly xi2  = read_base_poly(3, "(1-z+x)*0.5");
+      base_poly xi3  = read_base_poly(3, "(1-z+y)*0.5");
+      base_poly x    = read_base_poly(3, "x");
+      base_poly y    = read_base_poly(3, "y");
+      base_poly z    = read_base_poly(3, "z");
+      base_poly ones = read_base_poly(3, "1");
+      base_poly un_z = read_base_poly(3, "1-z");
+      base_rational_fraction Q(read_base_poly(3, "1"), un_z);
 
       p->base()[ 0] = Q*Q*xi0*xi1*(x*y-z*un_z);
       p->base()[ 1] = -Q*Q*xi0*xi1*xi2*y*4.;
@@ -1313,7 +1315,7 @@ namespace getfem {
       p->base()[10] = Q*z*xi1*xi2*4.;
       p->base()[11] = Q*z*xi3*xi0*4.;
       p->base()[12] = Q*z*xi2*xi3*4.;
-      p->base()[13] = bgeot::read_base_poly(3, "z*(2*z-1)");
+      p->base()[13] = read_base_poly(3, "z*(2*z-1)");
 
       p->add_node(lag_dof, base_small_vector(-1.0, -1.0, 0.0));
       p->add_node(lag_dof, base_small_vector( 0.0, -1.0, 0.0));
@@ -1338,8 +1340,7 @@ namespace getfem {
 
 
   static pfem pyramidal_pk_fem
-  (fem_param_list &params,
-   std::vector<dal::pstatic_stored_object> &dependencies) {
+  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps) {
     GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
     short_type k = 2;
     if (params.size() > 0) {
@@ -1347,14 +1348,13 @@ namespace getfem {
       k = dim_type(::floor(params[0].num() + 0.01));
     }
     pfem p = build_pyramidal_pk_fem(k, false);
-    dependencies.push_back(p->ref_convex(0));
-    dependencies.push_back(p->node_tab(0));
+    deps.push_back(p->ref_convex(0));
+    deps.push_back(p->node_tab(0));
     return p;
   }
 
   static pfem pyramidal_disc_pk_fem
-  (fem_param_list &params,
-   std::vector<dal::pstatic_stored_object> &dependencies) {
+  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps) {
     GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
     short_type k = 2;
     if (params.size() > 0) {
@@ -1362,18 +1362,18 @@ namespace getfem {
       k = dim_type(::floor(params[0].num() + 0.01));
     }
     pfem p = build_pyramidal_pk_fem(k, true);
-    dependencies.push_back(p->ref_convex(0));
-    dependencies.push_back(p->node_tab(0));
+    deps.push_back(p->ref_convex(0));
+    deps.push_back(p->node_tab(0));
     return p;
   }
 
-   /* ******************************************************************** */
-   /*        P1 element with a bubble base fonction on a face              */
-   /* ******************************************************************** */
+  /* ******************************************************************** */
+  /*        P1 element with a bubble base fonction on a face              */
+  /* ******************************************************************** */
 
-   struct P1_wabbfoaf_ : public PK_fem_ {
-     P1_wabbfoaf_(dim_type nc);
-   };
+  struct P1_wabbfoaf_ : public PK_fem_ {
+    P1_wabbfoaf_(dim_type nc);
+  };
 
   P1_wabbfoaf_::P1_wabbfoaf_(dim_type nc) : PK_fem_(nc, 1) {
     is_lag = false; es_degree = 2;
@@ -3136,7 +3136,7 @@ namespace getfem {
     base_node pt(3);
     for (unsigned k = 0; k < 5; ++k) {
       for (unsigned i = 0; i < 4; ++i) {
-        base_[k*4+i] = bgeot::read_base_poly(3, s);
+        base_[k*4+i] = read_base_poly(3, s);
         pt[0] = pt[1] = pt[2] = ((k == 4) ? 1.0/3.0 : 0.0);
         if (k == 4 && i) pt[i-1] = 0.0;
         if (k < 4 && k) pt[k-1] = 1.0;
@@ -3308,7 +3308,7 @@ namespace getfem {
     base_node pt(2);
     for (unsigned k = 0; k < 7; ++k) {
       for (unsigned i = 0; i < 3; ++i) {
-        base_[k*3+i] = bgeot::read_base_poly(2, s);
+        base_[k*3+i] = read_base_poly(2, s);
         if (k == 6) {
           pt[0] = pt[1] = 0.5; if (i) pt[i-1] = 0.0;
           add_node(normal_derivative_dof(2), pt);
@@ -3416,7 +3416,7 @@ namespace getfem {
                         "((x+y)^2 - x - y)*sqrt(2)/2;  x*(x-1);  y*(y-1);");
 
     for (unsigned k = 0; k < 6; ++k)
-      base_[k] = bgeot::read_base_poly(2, s);
+      base_[k] = read_base_poly(2, s);
 
     add_node(lagrange_dof(2), base_node(0.0, 0.0));
     add_node(lagrange_dof(2), base_node(1.0, 0.0));
@@ -3556,7 +3556,7 @@ namespace getfem {
     std::stringstream name;
 
     // Identifying if it is a torus structure
-    if(bgeot::is_torus_geom_trans(pgt) && n == 3) n = 2;
+    if (bgeot::is_torus_geom_trans(pgt) && n == 3) n = 2;
 
     /* Identifying P1-simplexes.                                          */
     if (nbp == n+1)
diff --git a/src/getfem_import.cc b/src/getfem_import.cc
index 6860d80..524f7ed 100644
--- a/src/getfem_import.cc
+++ b/src/getfem_import.cc
@@ -123,7 +123,7 @@ namespace getfem {
         nodes.resize(6);
       } break;
       case 7: { /* PYRAMID */
-       nodes.resize(5);
+        nodes.resize(5);
       } break;
       case 8: { /* 2ND ORDER LINE */
         nodes.resize(3);
@@ -329,10 +329,10 @@ namespace getfem {
         size_type j;
         f >> j;
         std::map<size_type, size_type>::iterator
-             it = msh_node_2_getfem_node.find(j);
+          it = msh_node_2_getfem_node.find(j);
         GMM_ASSERT1(it != msh_node_2_getfem_node.end(),
-                   "Invalid node ID " << j << " in gmsh element "
-                   << (ci.id + 1));
+                    "Invalid node ID " << j << " in gmsh element "
+                    << (ci.id + 1));
         ci.nodes[i] = it->second;
       }
       if(ci.type != 15) ci.set_pgt();
@@ -715,11 +715,11 @@ namespace getfem {
         for (size_type i=0; i < nnode; ++i) {
           size_type j;
           s >> j;
-         std::map<size_type, size_type>::iterator
-           it = msh_node_2_getfem_node.find(j);
+          std::map<size_type, size_type>::iterator
+            it = msh_node_2_getfem_node.find(j);
           GMM_ASSERT1(it != msh_node_2_getfem_node.end(),
-                     "Invalid node ID " << j << " in GiD element " << cv_id);
-         cv_nodes[i] = it->second;
+                      "Invalid node ID " << j << " in GiD element " << cv_id);
+          cv_nodes[i] = it->second;
         }
         getfem_cv_nodes.resize(nnode);
         for (size_type i=0; i < nnode; ++i) {
@@ -1105,7 +1105,7 @@ namespace getfem {
                                            getfem_cv_nodes.begin()));
             if (itype < elt_cnt.size())
               elt_cnt[itype] += 1;
-           }
+          }
         }
       }
     }
@@ -1157,12 +1157,14 @@ namespace getfem {
 
     /*NE: nombre d'elements (premier nombre du fichier .noboite)
       NP: nombre de points (deuxieme nombre du fichier .noboite)
-      ligne_debut_NP: la ligne commence la liste des coordonnees des points 
dans le            fichier .noboite*/
+      ligne_debut_NP: la ligne commence la liste des coordonnees des points 
dans le
+                      fichier .noboite*/
 
     f >> NE>>NP;
     ligne_debut_NP=NE*4/6+3;
 
-    //passer 3 premiers lignes du fichier .noboite (la liste des elements 
commence a la                quatrieme ligne)
+    //passer 3 premiers lignes du fichier .noboite (la liste des elements 
commence a la
+    //quatrieme ligne)
     string contenu;
     for (i=1; i<=ligne_debut_NP; i++){
       getline(f, contenu);
@@ -1211,16 +1213,16 @@ namespace getfem {
 
     if(fichier_GiD)  // si l'ouverture a r�ussi
       {
-       // instructions
-       fichier_GiD.close();  // on referme le fichier
+        // instructions
+        fichier_GiD.close();  // on referme le fichier
       }
     else  // sinon
       cerr << "Erreur � l'ouverture !" << endl;
 
     if(f)  // si l'ouverture a r�ussi
       {
-       // instructions
-       //f.close();  // on referme le fichier
+        // instructions
+        //f.close();  // on referme le fichier
       }
     else  // sinon
       cerr << "Erreur � l'ouverture !" << endl;
diff --git a/src/getfem_integration.cc b/src/getfem_integration.cc
index 047fc3f..85f5233 100644
--- a/src/getfem_integration.cc
+++ b/src/getfem_integration.cc
@@ -33,10 +33,10 @@
 namespace getfem {
 
   /*
-   * dummy integration method 
+   * dummy integration method
    */
   static pintegration_method im_none(im_param_list &params,
-                              std::vector<dal::pstatic_stored_object> &) {
+                               std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 0, "IM_NONE does not accept any parameter");
     return std::make_shared<integration_method>();
   }
@@ -49,7 +49,7 @@ namespace getfem {
       hum->resize(i);
       bgeot::power_index mi(P.dim()); mi[P.dim()-1] = P.degree();
       for (size_type k = i; k > j; --k, --mi)
-       (*hum)[k-1] = int_monomial(mi);
+        (*hum)[k-1] = int_monomial(mi);
     }
     base_poly::const_iterator it = P.begin(), ite = P.end();
     std::vector<long_scalar_type>::const_iterator itb = int_monomials.begin();
@@ -68,11 +68,11 @@ namespace getfem {
       hum->resize(i);
       bgeot::power_index mi(P.dim()); mi[P.dim()-1] = P.degree();
       for (size_type k = i; k > j; --k, --mi)
-       (*hum)[k-1] = int_monomial_on_face(mi, f);
+        (*hum)[k-1] = int_monomial_on_face(mi, f);
     }
     base_poly::const_iterator it = P.begin(), ite = P.end();
     std::vector<long_scalar_type>::const_iterator itb = hum->begin();
-    for ( ; it != ite; ++it, ++itb) 
+    for ( ; it != ite; ++it, ++itb)
       res += long_scalar_type(*it) * long_scalar_type(*itb);
     return res;
   }
@@ -85,8 +85,8 @@ namespace getfem {
 
     long_scalar_type int_monomial(const bgeot::power_index &power) const;
 
-    long_scalar_type int_monomial_on_face(const bgeot::power_index &power, 
-                                    short_type f) const;
+    long_scalar_type int_monomial_on_face(const bgeot::power_index &power,
+                                          short_type f) const;
 
     simplex_poly_integration_(bgeot::pconvex_structure c)
       { cvs = c;  int_face_monomials.resize(c->nb_faces()); }
@@ -101,38 +101,39 @@ namespace getfem {
       itme = power.end();
     for ( ; itm != itme; ++itm)
       for (int k = 1; k <= *itm; ++k, ++fa)
-       res *= long_scalar_type(k) / long_scalar_type(fa);
-    
+        res *= long_scalar_type(k) / long_scalar_type(fa);
+
     for (int k = 0; k < cvs->dim(); k++) { res /= long_scalar_type(fa); fa++; }
     return res;
   }
-  
+
   long_scalar_type simplex_poly_integration_::int_monomial_on_face
   (const bgeot::power_index &power, short_type f) const {
     long_scalar_type res = LONG_SCAL(0);
-    
+
     if (f == 0 || power[f-1] == 0.0) {
       res = (f == 0) ? sqrt(long_scalar_type(cvs->dim())) : LONG_SCAL(1);
       short_type fa = 1;
       bgeot::power_index::const_iterator itm = power.begin(),
-       itme = power.end();
+        itme = power.end();
       for ( ; itm != itme; ++itm)
-       for (int k = 1; k <= *itm; ++k, ++fa)
-         res *= long_scalar_type(k) / long_scalar_type(fa);
-      
+        for (int k = 1; k <= *itm; ++k, ++fa)
+          res *= long_scalar_type(k) / long_scalar_type(fa);
+
       for (int k = 1; k < cvs->dim(); k++) { res/=long_scalar_type(fa); fa++; }
     }
     return res;
   }
 
-  static pintegration_method exact_simplex(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  exact_simplex(im_param_list &params,
+                std::vector<dal::pstatic_stored_object> &dependencies) {
     GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
-               << params.size() << " should be 1.");
+                << params.size() << " should be 1.");
     GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
     int n = int(::floor(params[0].num() + 0.01));
     GMM_ASSERT1(n > 0 && n < 100 && double(n) == params[0].num(),
-               "Bad parameters");
+                "Bad parameters");
     dependencies.push_back(bgeot::simplex_structure(dim_type(n)));
     ppoly_integration ppi = std::make_shared<simplex_poly_integration_>
       (bgeot::simplex_structure(dim_type(n)));
@@ -148,8 +149,8 @@ namespace getfem {
 
     long_scalar_type int_monomial(const bgeot::power_index &power) const;
 
-    long_scalar_type int_monomial_on_face(const bgeot::power_index &power, 
-                                    short_type f) const;
+    long_scalar_type int_monomial_on_face(const bgeot::power_index &power,
+                                          short_type f) const;
 
     plyint_mul_structure_(ppoly_integration a, ppoly_integration b);
   };
@@ -161,7 +162,7 @@ namespace getfem {
     std::copy(power.begin() + cv1->dim(), power.end(), mi2.begin());
     return cv1->int_monomial(mi1) * cv2->int_monomial(mi2);
   }
-  
+
   long_scalar_type plyint_mul_structure_::int_monomial_on_face
   (const bgeot::power_index &power, short_type f) const {
     bgeot::power_index mi1(cv1->dim()), mi2(cv2->dim());
@@ -172,30 +173,31 @@ namespace getfem {
       return cv1->int_monomial_on_face(mi1,f) * cv2->int_monomial(mi2);
     else
       return cv1->int_monomial(mi1)
-       * cv2->int_monomial_on_face(mi2, short_type(f-nfx));
+             * cv2->int_monomial_on_face(mi2, short_type(f-nfx));
   }
-  
-  plyint_mul_structure_::plyint_mul_structure_(ppoly_integration a, 
-                                              ppoly_integration b) {
+
+  plyint_mul_structure_::plyint_mul_structure_(ppoly_integration a,
+                                               ppoly_integration b) {
     cv1 = a; cv2 = b;
     cvs = bgeot::convex_product_structure(cv1->structure(),
-                                         cv2->structure());
+                                          cv2->structure());
     int_face_monomials.resize(cvs->nb_faces());
   }
 
-  static pintegration_method product_exact(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  product_exact(im_param_list &params,
+                std::vector<dal::pstatic_stored_object> &dependencies) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
-               << params.size() << " should be 2.");
-    GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1, 
-               "Bad type of parameters");
+                << params.size() << " should be 2.");
+    GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
+                "Bad type of parameters");
     pintegration_method a = params[0].method();
     pintegration_method b = params[1].method();
     GMM_ASSERT1(a->type() == IM_EXACT && b->type() == IM_EXACT,
-               "Bad parameters");
+                "Bad parameters");
     dependencies.push_back(a); dependencies.push_back(b);
     dependencies.push_back(bgeot::convex_product_structure(a->structure(),
-                                                          b->structure()));
+                                                           b->structure()));
     ppoly_integration ppi = std::make_shared<plyint_mul_structure_>
       (a->exact_method(), b->exact_method());
     return  std::make_shared<integration_method>(ppi);
@@ -205,36 +207,38 @@ namespace getfem {
   /* integration on parallelepiped.                                       */
   /* ******************************************************************** */
 
-  static pintegration_method exact_parallelepiped(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &) {
+  static pintegration_method
+  exact_parallelepiped(im_param_list &params,
+                       std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
-               << params.size() << " should be 1.");
+                << params.size() << " should be 1.");
     GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
     int n = int(::floor(params[0].num() + 0.01));
     GMM_ASSERT1(n > 0 && n < 100 && double(n) == params[0].num(),
-               "Bad parameters");
+                "Bad parameters");
 
     std::stringstream name;
     if (n == 1)
       name << "IM_EXACT_SIMPLEX(1)";
-    else 
+    else
       name << "IM_PRODUCT(IM_EXACT_PARALLELEPIPED(" << n-1
-          << "),IM_EXACT_SIMPLEX(1)))";
+           << "),IM_EXACT_SIMPLEX(1)))";
     return int_method_descriptor(name.str());
   }
 
-  static pintegration_method exact_prism(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &) {
+  static pintegration_method
+  exact_prism(im_param_list &params,
+              std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
-               << params.size() << " should be 1.");
+                << params.size() << " should be 1.");
     GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
     int n = int(::floor(params[0].num() + 0.01));
     GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
-               "Bad parameters");
+                "Bad parameters");
 
     std::stringstream name;
     name << "IM_PRODUCT(IM_EXACT_SIMPLEX(" << n-1
-        << "),IM_EXACT_SIMPLEX(1))";
+         << "),IM_EXACT_SIMPLEX(1))";
     return int_method_descriptor(name.str());
   }
 
@@ -243,116 +247,116 @@ namespace getfem {
   /* ********************************************************************* */
 
   void approx_integration::add_point(const base_node &pt,
-                                    scalar_type w,short_type f) {
+                                     scalar_type w,short_type f) {
     GMM_ASSERT1(!valid, "Impossible to modify a valid integration method.");
     if (gmm::abs(w) > 1.0E-15) {
       ++f;
       GMM_ASSERT1(f <= cvr->structure()->nb_faces(), "Wrong argument.");
       size_type i = pt_to_store[f].search_node(pt);
       if (i == size_type(-1)) {
-       i = pt_to_store[f].add_node(pt);
-       int_coeffs.resize(int_coeffs.size() + 1); 
-       for (size_type j = f; j <= cvr->structure()->nb_faces(); ++j)
-         repartition[j] += 1;
-       for (size_type j = int_coeffs.size(); j > repartition[f]; --j)
-         int_coeffs[j-1] = int_coeffs[j-2];
-       int_coeffs[repartition[f]-1] = 0.0;
+        i = pt_to_store[f].add_node(pt);
+        int_coeffs.resize(int_coeffs.size() + 1);
+        for (size_type j = f; j <= cvr->structure()->nb_faces(); ++j)
+          repartition[j] += 1;
+        for (size_type j = int_coeffs.size(); j > repartition[f]; --j)
+          int_coeffs[j-1] = int_coeffs[j-2];
+        int_coeffs[repartition[f]-1] = 0.0;
       }
       int_coeffs[((f == 0) ? 0 : repartition[f-1]) + i] += w;
     }
   }
 
   void approx_integration::add_point_norepeat(const base_node &pt,
-                                             scalar_type w,
-                                             short_type f) {
+                                              scalar_type w,
+                                              short_type f) {
     short_type f2 = f; ++f2;
     if (pt_to_store[f2].search_node(pt) == size_type(-1)) add_point(pt,w,f);
   }
 
   void approx_integration::add_point_full_symmetric(base_node pt,
-                                                   scalar_type w) {
+                                                    scalar_type w) {
     dim_type n = cvr->structure()->dim();
     dim_type k;
     base_node pt2(n);
     if (n+1 == cvr->structure()->nb_faces()) {
       // simplices
-      // for a point at (x,y) you have to consider the other points 
+      // for a point at (x,y) you have to consider the other points
       // at (y,x) (x,1-x-y) (1-x-y,x) (y,1-x-y) (1-x-y,y)
       base_node pt3(n+1);
       std::copy(pt.begin(), pt.end(), pt3.begin());
       pt3[n] = 1;  for (k = 0; k < n; ++k) pt3[n] -= pt[k];
       std::vector<int> ind(n); std::fill(ind.begin(), ind.end(), 0);
-      std::vector<bool> ind2(n+1); 
+      std::vector<bool> ind2(n+1);
       for (;;) {
-       
-       std::fill(ind2.begin(), ind2.end(), false);
-       bool good = true;
-       for (k = 0; k < n; ++k) 
-         if (ind2[ind[k]]) { good = false; break; } else ind2[ind[k]] = true;
-       if (good) {
-         for (k = 0; k < n; ++k) pt2[k] = pt3[ind[k]];
-         add_point_norepeat(pt2, w);
-       }
-       ind[0]++; k = 0;
-       while(ind[k] == n+1) { ind[k++] = 0; if (k == n) return; ind[k]++; }
+
+        std::fill(ind2.begin(), ind2.end(), false);
+        bool good = true;
+        for (k = 0; k < n; ++k)
+          if (ind2[ind[k]]) { good = false; break; } else ind2[ind[k]] = true;
+        if (good) {
+          for (k = 0; k < n; ++k) pt2[k] = pt3[ind[k]];
+          add_point_norepeat(pt2, w);
+        }
+        ind[0]++; k = 0;
+        while(ind[k] == n+1) { ind[k++] = 0; if (k == n) return; ind[k]++; }
       }
     }
 
     else if (cvr->structure()->nb_points() == (size_type(1) << n)) {
       // parallelepidedic
       for (size_type i = 0; i < (size_type(1) << n); ++i) {
-       for (k = 0; k < n; ++k)
-         if (i & (size_type(1) << k)) pt2[k]=pt[k]; else pt2[k] = 1.0-pt[k];
-       add_point_norepeat(pt2, w);
+        for (k = 0; k < n; ++k)
+          if (i & (size_type(1) << k)) pt2[k]=pt[k]; else pt2[k] = 1.0-pt[k];
+        add_point_norepeat(pt2, w);
       }
     }
     else
       GMM_ASSERT1(false, "Fully symmetric option is only valid for"
-                 "simplices and parallelepipedic elements");
+                  "simplices and parallelepipedic elements");
   }
 
   void approx_integration::add_method_on_face(pintegration_method ppi,
-                                             short_type f) {
+                                              short_type f) {
     papprox_integration pai = ppi->approx_method();
     GMM_ASSERT1(!valid, "Impossible to modify a valid integration method.");
     GMM_ASSERT1(pai->structure() == structure()->faces_structure()[f],
-               "structures missmatch");
+                "structures missmatch");
     GMM_ASSERT1(ppi->type() == IM_APPROX, "Impossible with an exact method.");
 
     dim_type N = pai->structure()->dim();
     scalar_type det = 1.0;
     base_node pt(N+1);
-    std::vector<base_node> pts(N);  
+    std::vector<base_node> pts(N);
     for (size_type i = 0; i < N; ++i)
       pts[i] = (cvr->dir_points_of_face(f))[i+1]
-       - (cvr->dir_points_of_face(f))[0];
+             - (cvr->dir_points_of_face(f))[0];
     if (N) {
       base_matrix a(N+1, N), b(N, N), tmp(N, N);
       for (dim_type i = 0; i < N+1; ++i)
-       for (dim_type j = 0; j < N; ++j)
-         a(i, j) = pts[j][i];
-      
+        for (dim_type j = 0; j < N; ++j)
+          a(i, j) = pts[j][i];
+
       gmm::mult(gmm::transposed(a), a, b);
       det = ::sqrt(gmm::abs(gmm::lu_det(b)));
     }
     for (size_type i = 0; i < pai->nb_points_on_convex(); ++i) {
       pt = (cvr->dir_points_of_face(f))[0];
       for (dim_type j = 0; j < N; ++j)
-       pt += pts[j] * ((*(pai->pintegration_points()))[i])[j];
+        pt += pts[j] * ((*(pai->pintegration_points()))[i])[j];
       add_point(pt, pai->coeff(i) * det, f);
     }
   }
 
-  void approx_integration::valid_method(void) {
+  void approx_integration::valid_method() {
     std::vector<base_node> ptab(int_coeffs.size());
     // std::vector<scalar_type> int_coeffs2(int_coeffs);
     size_type i = 0;
     for (short_type f = 0; f <= cvr->structure()->nb_faces(); ++f) {
       // size_type j = i;
       for (PT_TAB::const_iterator it = pt_to_store[f].begin();
-          it != pt_to_store[f].end(); ++it /* , ++j */) {
-       // int_coeffs[i] = int_coeffs2[j];
-       ptab[i++] = *it;
+           it != pt_to_store[f].end(); ++it /* , ++j */) {
+        // int_coeffs[i] = int_coeffs2[j];
+        ptab[i++] = *it;
       }
     }
     GMM_ASSERT1(i == int_coeffs.size(), "internal error.");
@@ -365,60 +369,61 @@ namespace getfem {
   /* ********************************************************************* */
   /* methods stored in getfem_im_list.h                                    */
   /* ********************************************************************* */
-  
+
   /// search a method in getfem_im_list.h
-  static pintegration_method im_list_integration(std::string name,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  im_list_integration(std::string name,
+                      std::vector<dal::pstatic_stored_object> &dependencies) {
     // cerr << "searching " << name << endl;
     for (int i = 0; i < NB_IM; ++i)
       if (!name.compare(im_desc_tab[i].method_name)) {
-       bgeot::pgeometric_trans pgt
-         = bgeot::geometric_trans_descriptor(im_desc_tab[i].geotrans_name);
-       dim_type N = pgt->structure()->dim();
-       base_node pt(N);
-       auto pai = std::make_shared<approx_integration>(pgt->convex_ref());
-       size_type fr = im_desc_tab[i].firstreal;
-       for (size_type j = 0; j < im_desc_tab[i].nb_points; ++j) {
-         for (dim_type k = 0; k < N; ++k)
-           pt[k] = atof(im_desc_real[fr+j*(N+1)+k]);
-         // pt[k] = LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+k]);
-         
-         switch (im_desc_node_type[im_desc_tab[i].firsttype + j]) {
-         case 2: {
+        bgeot::pgeometric_trans pgt
+          = bgeot::geometric_trans_descriptor(im_desc_tab[i].geotrans_name);
+        dim_type N = pgt->structure()->dim();
+        base_node pt(N);
+        auto pai = std::make_shared<approx_integration>(pgt->convex_ref());
+        size_type fr = im_desc_tab[i].firstreal;
+        for (size_type j = 0; j < im_desc_tab[i].nb_points; ++j) {
+          for (dim_type k = 0; k < N; ++k)
+            pt[k] = atof(im_desc_real[fr+j*(N+1)+k]);
+          // pt[k] = LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+k]);
+
+          switch (im_desc_node_type[im_desc_tab[i].firsttype + j]) {
+          case 2: {
             base_node pt2(pt.size());
             for (bgeot::permutation p(pt.size()); !p.finished(); ++p) {
               p.apply_to(pt,pt2);
-              pai->add_point_full_symmetric(pt2, 
-                                           atof(im_desc_real[fr+j*(N+1)+N]));
-             // pai->add_point_full_symmetric(pt2,
-             //                LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
+              pai->add_point_full_symmetric(pt2,
+                                            atof(im_desc_real[fr+j*(N+1)+N]));
+              // pai->add_point_full_symmetric(pt2,
+              //                LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
             }
-         } break;
-         case 1: {
-           pai->add_point_full_symmetric(pt,atof(im_desc_real[fr+j*(N+1)+N]));
-           // pai->add_point_full_symmetric(pt,
-           //                   LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
-         } break;
-         case 0: {
-           pai->add_point(pt, atof(im_desc_real[fr+j*(N+1)+N]));
-           // pai->add_point(pt,LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
-         } break;
-         default: GMM_ASSERT1(false, "");
-         }
-       }
-
-       for (short_type f = 0; N > 0 && f < pgt->structure()->nb_faces(); ++f)
-         pai->add_method_on_face
-           (int_method_descriptor
-            (im_desc_face_meth[im_desc_tab[i].firstface + f]), f);
-
-       pai->valid_method();
+          } break;
+          case 1: {
+            pai->add_point_full_symmetric(pt,atof(im_desc_real[fr+j*(N+1)+N]));
+            // pai->add_point_full_symmetric(pt,
+            //                   LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
+          } break;
+          case 0: {
+            pai->add_point(pt, atof(im_desc_real[fr+j*(N+1)+N]));
+            // pai->add_point(pt,LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
+          } break;
+          default: GMM_ASSERT1(false, "");
+          }
+        }
+
+        for (short_type f = 0; N > 0 && f < pgt->structure()->nb_faces(); ++f)
+          pai->add_method_on_face
+            (int_method_descriptor
+             (im_desc_face_meth[im_desc_tab[i].firstface + f]), f);
+
+        pai->valid_method();
         // cerr << "finding " << name << endl;
 
-       pintegration_method p(std::make_shared<integration_method>(pai));
-       dependencies.push_back(p->approx_method()->ref_convex());
-       dependencies.push_back(p->approx_method()->pintegration_points());
-       return p;
+        pintegration_method p(std::make_shared<integration_method>(pai));
+        dependencies.push_back(p->approx_method()->ref_convex());
+        dependencies.push_back(p->approx_method()->pintegration_points());
+        return p;
       }
     return pintegration_method();
   }
@@ -450,7 +455,7 @@ namespace getfem {
           (base_poly(1,1,0) * polynomials[nb_lp-1]
            * ((2.0 * long_scalar_type(nb_lp) - 1.0) / long_scalar_type(nb_lp)))
           + (polynomials[nb_lp-2]
-          * ((1.0 - long_scalar_type(nb_lp)) / long_scalar_type(nb_lp)));
+           * ((1.0 - long_scalar_type(nb_lp)) / long_scalar_type(nb_lp)));
         roots[nb_lp].resize(nb_lp);
         roots[nb_lp][nb_lp/2] = 0.0;
         long_scalar_type a = -1.0, b, c, d, e, cv, ev, ecart, ecart2;
@@ -474,7 +479,7 @@ namespace getfem {
           roots[nb_lp][nb_lp-k-1] = -c;
           a = b;
         }
-      } 
+      }
     }
   };
 
@@ -484,30 +489,30 @@ namespace getfem {
 
   gauss_approx_integration_::gauss_approx_integration_(short_type nbpt) {
     GMM_ASSERT1(nbpt <= 32000, "too much points");
-    
+
     cvr = bgeot::simplex_of_reference(1);
     std::vector<base_node> int_points(nbpt+2);
     int_coeffs.resize(nbpt+2);
     repartition.resize(3);
-    repartition[0] = nbpt; 
+    repartition[0] = nbpt;
     repartition[1] = nbpt + 1;
-    repartition[2] = nbpt + 2; 
-    
+    repartition[2] = nbpt + 2;
+
     Legendre_polynomials Lp;
     Lp.init(nbpt);
-    
+
     for (short_type i = 0; i < nbpt; ++i) {
       int_points[i].resize(1);
       long_scalar_type lr = Lp.roots[nbpt][i];
       int_points[i][0] = 0.5 + 0.5 * bgeot::to_scalar(lr);
       int_coeffs[i] = bgeot::to_scalar((1.0 - gmm::sqr(lr))
-       / gmm::sqr( long_scalar_type(nbpt)
-                   * (Lp.polynomials[nbpt-1].eval(&lr))));
+                    / gmm::sqr( long_scalar_type(nbpt)
+                               * (Lp.polynomials[nbpt-1].eval(&lr))));
     }
-    
+
     int_points[nbpt].resize(1);
     int_points[nbpt][0] = 1.0; int_coeffs[nbpt] = 1.0;
-    
+
     int_points[nbpt+1].resize(1);
     int_points[nbpt+1][0] = 0.0; int_coeffs[nbpt+1] = 1.0;
     pint_points = bgeot::store_point_tab(int_points);
@@ -515,14 +520,15 @@ namespace getfem {
   }
 
 
-  static pintegration_method gauss1d(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  gauss1d(im_param_list &params,
+          std::vector<dal::pstatic_stored_object> &dependencies) {
     GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
-               << params.size() << " should be 1.");
+                << params.size() << " should be 1.");
     GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
     int n = int(::floor(params[0].num() + 0.01));
     GMM_ASSERT1(n >= 0 && n < 32000 && double(n) == params[0].num(),
-               "Bad parameters");
+                "Bad parameters");
     if (n & 1) {
       std::stringstream name;
       name << "IM_GAUSS1D(" << n-1 << ")";
@@ -530,7 +536,7 @@ namespace getfem {
     }
     else {
       papprox_integration
-       pai = std::make_shared<gauss_approx_integration_>(short_type(n/2 + 1));
+        pai = std::make_shared<gauss_approx_integration_>(short_type(n/2 + 1));
       pintegration_method p = std::make_shared<integration_method>(pai);
       dependencies.push_back(p->approx_method()->ref_convex());
       dependencies.push_back(p->approx_method()->pintegration_points());
@@ -553,97 +559,98 @@ namespace getfem {
     : approx_integration(bgeot::simplex_of_reference(nc)) {
     size_type R = bgeot::alpha(nc,k);
 
-    base_node c(nc); 
+    base_node c(nc);
     if (nc == 0) {
       add_point(c, scalar_type(1));
     }
     else {
-      
+
       std::stringstream name;
       name << "IM_EXACT_SIMPLEX(" << int(nc) << ")";
       ppoly_integration ppi = 
int_method_descriptor(name.str())->exact_method();
-      
+
       size_type sum = 0, l;
       c.fill(scalar_type(0.0));
       if (k == 0) c.fill(1.0 / scalar_type(nc+1));
-      
+
       gmm::dense_matrix<long_scalar_type> M(R, R);
       std::vector<long_scalar_type> F(R), U(R);
       std::vector<bgeot::power_index> base(R);
       std::vector<base_node> nodes(R);
-      
+
       bgeot::power_index pi(nc);
-      
+
       for (size_type r = 0; r < R; ++r, ++pi) {
-       base[r] = pi; nodes[r] = c;
-       if (k != 0 && nc > 0) {
-         l = 0; c[l] += 1.0 / scalar_type(k); sum++;
-         while (sum > k) {
-           sum -= int(floor(0.5+(c[l] * k)));
-           c[l] = 0.0; l++; if (l == nc) break;
-           c[l] += 1.0 / scalar_type(k); sum++;
-         }
-       }
+        base[r] = pi; nodes[r] = c;
+        if (k != 0 && nc > 0) {
+          l = 0; c[l] += 1.0 / scalar_type(k); sum++;
+          while (sum > k) {
+            sum -= int(floor(0.5+(c[l] * k)));
+            c[l] = 0.0; l++; if (l == nc) break;
+            c[l] += 1.0 / scalar_type(k); sum++;
+          }
+        }
       }
 
 //       if (nc == 1) {
-//     M = bgeot::vsmatrix<long_scalar_type>((R+1)/2, (R+1)/2);
-//     U = F = bgeot::vsvector<long_scalar_type>((R+1)/2);
-//     gmm::clear(M);
+//         M = bgeot::vsmatrix<long_scalar_type>((R+1)/2, (R+1)/2);
+//         U = F = bgeot::vsvector<long_scalar_type>((R+1)/2);
+//         gmm::clear(M);
 //       }
 
       for (size_type r = 0; r < R; ++r) {
-//     if (nc == 1) {
-//       if (r < (R+1)/2) {
-//         F[r] = ppi->int_monomial(base[R-1-r]);
-//         cout << "F[" << r << "] = " << F[r] << endl; 
-//       }
-//     }
-//     else {
-         F[r] = ppi->int_monomial(base[r]);
-         //cout << "F[" << r << "] = " << F[r] << endl;
-//     }
-       for (size_type q = 0; q < R; ++q) {
-//       if (nc == 1) {
-//         if (r < (R+1)/2) {
-//           if (q < (R+1)/2) 
-//             M(r, q) += bgeot::eval_monomial(base[R-1-r], nodes[q].begin());
-//           else
-//             M(r, R-1-q) += bgeot::eval_monomial(base[R-1-r],
-//       nodes[q].begin());
-//         }
-//       }
-//       else
-           M(r, q) = bgeot::eval_monomial(base[r], nodes[q].begin());
-       }
+//         if (nc == 1) {
+//           if (r < (R+1)/2) {
+//             F[r] = ppi->int_monomial(base[R-1-r]);
+//             cout << "F[" << r << "] = " << F[r] << endl;
+//           }
+//         }
+//         else {
+          F[r] = ppi->int_monomial(base[r]);
+          //cout << "F[" << r << "] = " << F[r] << endl;
+//         }
+        for (size_type q = 0; q < R; ++q) {
+//           if (nc == 1) {
+//             if (r < (R+1)/2) {
+//               if (q < (R+1)/2)
+//                 M(r, q) += bgeot::eval_monomial(base[R-1-r], 
nodes[q].begin());
+//               else
+//                 M(r, R-1-q) += bgeot::eval_monomial(base[R-1-r],
+//          nodes[q].begin());
+//             }
+//           }
+//           else
+            M(r, q) = bgeot::eval_monomial(base[r], nodes[q].begin());
+        }
       }
-      
+
       gmm::lu_solve(M, U, F);
       for (size_type r = 0; r < R; ++r)
-       add_point(nodes[r], bgeot::to_scalar(U[r]));
-      
+        add_point(nodes[r], bgeot::to_scalar(U[r]));
+
       std::stringstream name2;
       name2 << "IM_NC(" << int(nc-1) << "," << int(k) << ")";
       for (short_type f = 0; f < structure()->nb_faces(); ++f)
-       add_method_on_face(int_method_descriptor(name2.str()), f);
+        add_method_on_face(int_method_descriptor(name2.str()), f);
     }
     valid_method();
   }
 
-  static pintegration_method Newton_Cotes(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  Newton_Cotes(im_param_list &params,
+               std::vector<dal::pstatic_stored_object> &dependencies) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
-               << params.size() << " should be 2.");
+                << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
-               "Bad type of parameters");
+                "Bad type of parameters");
     int n = int(::floor(params[0].num() + 0.01));
     int k = int(::floor(params[1].num() + 0.01));
     GMM_ASSERT1(n >= 0 && n < 100 && k >= 0 && k <= 150 &&
-               double(n) == params[0].num() && double(k) == params[1].num(),
-               "Bad parameters");
+                double(n) == params[0].num() && double(k) == params[1].num(),
+                "Bad parameters");
     papprox_integration
       pai = std::make_shared<Newton_Cotes_approx_integration_>(dim_type(n),
-                                                              short_type(k));
+                                                               short_type(k));
     pintegration_method p = std::make_shared<integration_method>(pai);
     dependencies.push_back(p->approx_method()->ref_convex());
     dependencies.push_back(p->approx_method()->pintegration_points());
@@ -661,7 +668,7 @@ namespace getfem {
 
 
   a_int_pro_integration::a_int_pro_integration(papprox_integration a,
-                                              papprox_integration b) {
+                                               papprox_integration b) {
     cvr = bgeot::convex_ref_product(a->ref_convex(), b->ref_convex());
     size_type n1 = a->nb_points_on_convex();
     size_type n2 = b->nb_points_on_convex();
@@ -672,12 +679,12 @@ namespace getfem {
     repartition[0] = n1 * n2;
     for (size_type i1 = 0; i1 < n1; ++i1)
       for (size_type i2 = 0; i2 < n2; ++i2) {
-       int_coeffs[i1 + i2 * n1] = a->coeff(i1) * b->coeff(i2);
-       int_points[i1 + i2 * n1].resize(dim());
-       std::copy(a->point(i1).begin(), a->point(i1).end(),
-                 int_points[i1 + i2 * n1].begin());
-       std::copy(b->point(i2).begin(), b->point(i2).end(),
-                 int_points[i1 + i2 * n1].begin() + a->dim());
+        int_coeffs[i1 + i2 * n1] = a->coeff(i1) * b->coeff(i2);
+        int_points[i1 + i2 * n1].resize(dim());
+        std::copy(a->point(i1).begin(), a->point(i1).end(),
+                  int_points[i1 + i2 * n1].begin());
+        std::copy(b->point(i2).begin(), b->point(i2).end(),
+                  int_points[i1 + i2 * n1].begin() + a->dim());
       }
     short_type f = 0;
     for (short_type f1 = 0; f1 < a->structure()->nb_faces(); ++f1, ++f) {
@@ -688,16 +695,16 @@ namespace getfem {
       int_points.resize(repartition[f+1]);
       int_coeffs.resize(repartition[f+1]);
       for (size_type i1 = 0; i1 < n1; ++i1)
-       for (size_type i2 = 0; i2 < n2; ++i2) {
-         int_coeffs[w + i1 + i2 * n1] = a->coeff_on_face(f1, i1)
-           * b->coeff(i2);
-         int_points[w + i1 + i2 * n1].resize(dim());
-         std::copy(a->point_on_face(f1, i1).begin(),
-                   a->point_on_face(f1, i1).end(),
-                   int_points[w + i1 + i2 * n1].begin());
-         std::copy(b->point(i2).begin(), b->point(i2).end(),
-                   int_points[w + i1 + i2 * n1].begin() + a->dim());
-       }
+        for (size_type i2 = 0; i2 < n2; ++i2) {
+          int_coeffs[w + i1 + i2 * n1] = a->coeff_on_face(f1, i1)
+            * b->coeff(i2);
+          int_points[w + i1 + i2 * n1].resize(dim());
+          std::copy(a->point_on_face(f1, i1).begin(),
+                    a->point_on_face(f1, i1).end(),
+                    int_points[w + i1 + i2 * n1].begin());
+          std::copy(b->point(i2).begin(), b->point(i2).end(),
+                    int_points[w + i1 + i2 * n1].begin() + a->dim());
+        }
     }
     for (short_type f2 = 0; f2 < b->structure()->nb_faces(); ++f2, ++f) {
       n1 = a->nb_points_on_convex();
@@ -707,46 +714,48 @@ namespace getfem {
       int_points.resize(repartition[f+1]);
       int_coeffs.resize(repartition[f+1]);
       for (size_type i1 = 0; i1 < n1; ++i1)
-       for (size_type i2 = 0; i2 < n2; ++i2) {
-         int_coeffs[w + i1 + i2 * n1] = a->coeff(i1)
-           * b->coeff_on_face(f2, i2);
-         int_points[w + i1 + i2 * n1].resize(dim());
-         std::copy(a->point(i1).begin(), a->point(i1).end(),
-                   int_points[w + i1 + i2 * n1].begin());
-         std::copy(b->point_on_face(f2, i2).begin(),
-                   b->point_on_face(f2, i2).end(),
-                   int_points[w + i1 + i2 * n1].begin() + a->dim());
-       }
+        for (size_type i2 = 0; i2 < n2; ++i2) {
+          int_coeffs[w + i1 + i2 * n1] = a->coeff(i1)
+            * b->coeff_on_face(f2, i2);
+          int_points[w + i1 + i2 * n1].resize(dim());
+          std::copy(a->point(i1).begin(), a->point(i1).end(),
+                    int_points[w + i1 + i2 * n1].begin());
+          std::copy(b->point_on_face(f2, i2).begin(),
+                    b->point_on_face(f2, i2).end(),
+                    int_points[w + i1 + i2 * n1].begin() + a->dim());
+        }
     }
     pint_points = bgeot::store_point_tab(int_points);
     valid = true;
   }
 
-  static pintegration_method product_approx(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  product_approx(im_param_list &params,
+                 std::vector<dal::pstatic_stored_object> &dependencies) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
-               << params.size() << " should be 2.");
+                << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
-               "Bad type of parameters");
+                "Bad type of parameters");
     pintegration_method a = params[0].method();
     pintegration_method b = params[1].method();
     GMM_ASSERT1(a->type() == IM_APPROX && b->type() == IM_APPROX,
-               "Bad parameters");
+                "Bad parameters");
     papprox_integration
       pai = std::make_shared<a_int_pro_integration>(a->approx_method(),
-                                                   b->approx_method());
+                                                    b->approx_method());
     pintegration_method p = std::make_shared<integration_method>(pai);
     dependencies.push_back(p->approx_method()->ref_convex());
     dependencies.push_back(p->approx_method()->pintegration_points());
     return p;
   }
 
-  static pintegration_method product_which(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  product_which(im_param_list &params,
+                std::vector<dal::pstatic_stored_object> &dependencies) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
-               << params.size() << " should be 2.");
+                << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
-               "Bad type of parameters");
+                "Bad type of parameters");
     pintegration_method a = params[0].method();
     pintegration_method b = params[1].method();
     if (a->type() == IM_EXACT || b->type() == IM_EXACT)
@@ -759,42 +768,44 @@ namespace getfem {
   /* integration on parallelepiped with Newton Cotes formulae              */
   /* ********************************************************************* */
 
-  static pintegration_method Newton_Cotes_para(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &) {
+  static pintegration_method
+  Newton_Cotes_para(im_param_list &params,
+                    std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
-               << params.size() << " should be 2.");
+                << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
-               "Bad type of parameters");
+                "Bad type of parameters");
     int n = int(::floor(params[0].num() + 0.01));
     int k = int(::floor(params[1].num() + 0.01));
     GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
-               double(n) == params[0].num() && double(k) == params[1].num(),
-               "Bad parameters");
+                double(n) == params[0].num() && double(k) == params[1].num(),
+                "Bad parameters");
 
     std::stringstream name;
     if (n == 1)
       name << "IM_NC(1," << k << ")";
-    else 
+    else
       name << "IM_PRODUCT(IM_NC_PARALLELEPIPED(" << n-1 << "," << k
-          << "),IM_NC(1," << k << "))";
+           << "),IM_NC(1," << k << "))";
     return int_method_descriptor(name.str());
   }
 
-  static pintegration_method Newton_Cotes_prism(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &) {
+  static pintegration_method
+  Newton_Cotes_prism(im_param_list &params,
+                     std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
-               << params.size() << " should be 2.");
+                << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
-               "Bad type of parameters");
+                "Bad type of parameters");
     int n = int(::floor(params[0].num() + 0.01));
     int k = int(::floor(params[1].num() + 0.01));
     GMM_ASSERT1(n > 1 && n < 100 && k >= 0 && k <= 150 &&
-               double(n) == params[0].num() && double(k) == params[1].num(),
-               "Bad parameters");
+                double(n) == params[0].num() && double(k) == params[1].num(),
+                "Bad parameters");
 
     std::stringstream name;
     name << "IM_PRODUCT(IM_NC(" << n-1 << "," << k << "),IM_NC(1,"
-        << k << "))";
+         << k << "))";
     return int_method_descriptor(name.str());
   }
 
@@ -802,24 +813,25 @@ namespace getfem {
   /* integration on parallelepiped with Gauss formulae                     */
   /* ********************************************************************* */
 
-  static pintegration_method Gauss_paramul(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &) {
+  static pintegration_method
+  Gauss_paramul(im_param_list &params,
+                std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
-               << params.size() << " should be 2.");
+                << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
-               "Bad type of parameters");
+                "Bad type of parameters");
     int n = int(::floor(params[0].num() + 0.01));
     int k = int(::floor(params[1].num() + 0.01));
     GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
-               double(n) == params[0].num() && double(k) == params[1].num(),
-               "Bad parameters");
+                double(n) == params[0].num() && double(k) == params[1].num(),
+                "Bad parameters");
 
     std::stringstream name;
     if (n == 1)
       name << "IM_GAUSS1D(" << k << ")";
-    else 
+    else
       name << "IM_PRODUCT(IM_GAUSS_PARALLELEPIPED(" << n-1 << "," << k
-          << "),IM_GAUSS1D(" << k << "))";
+           << "),IM_GAUSS1D(" << k << "))";
     return int_method_descriptor(name.str());
   }
 
@@ -828,8 +840,8 @@ namespace getfem {
   /* ******************************************************************** */
 
   struct quasi_polar_integration : public approx_integration {
-    quasi_polar_integration(papprox_integration base_im, 
-                           size_type ip1, size_type ip2=size_type(-1)) : 
+    quasi_polar_integration(papprox_integration base_im,
+                            size_type ip1, size_type ip2=size_type(-1)) :
       approx_integration
       ((base_im->structure() == bgeot::parallelepiped_structure(3)) ?
        bgeot::pyramidal_element_of_reference(1)
@@ -839,11 +851,11 @@ namespace getfem {
       enum { SQUARE, PRISM, TETRA_CYL, PRISM2, PYRAMID } what;
       if (N == 2) what = SQUARE;
       else if (base_im->structure() == bgeot::prism_structure(3))
-       what = (ip2 == size_type(-1) || ip1 == ip2) ? PRISM2 : PRISM;
+        what = (ip2 == size_type(-1) || ip1 == ip2) ? PRISM2 : PRISM;
       else if (base_im->structure() == bgeot::simplex_structure(3))
-       what = TETRA_CYL;
+        what = TETRA_CYL;
       else if (base_im->structure() == bgeot::parallelepiped_structure(3))
-       what = PYRAMID;
+        what = PYRAMID;
       else GMM_ASSERT1(false, "Incoherent integration method");
 
       // The first geometric transformation collapse a face of
@@ -855,8 +867,8 @@ namespace getfem {
       bgeot::pgeometric_trans pgt2 = bgeot::simplex_geotrans(N, 1);
       std::vector<base_node> nodes2(N+1);
       if (what == PYRAMID) {
-       pgt2 = bgeot::pyramidal_geotrans(1);
-       nodes2.resize(5);
+        pgt2 = bgeot::pyramidal_geotrans(1);
+        nodes2.resize(5);
       }
       std::vector<size_type> other_nodes; // for the construction of node2
       bgeot::pgeometric_trans pgt3 = bgeot::simplex_geotrans(N, 1);
@@ -864,162 +876,163 @@ namespace getfem {
 
       switch (what) {
       case SQUARE :
-       nodes1[3] = nodes1[1];
-       nodes2[ip1] = nodes1[1]; ip2 = ip1;
-       other_nodes.push_back(0);
-       other_nodes.push_back(2);
-       break;
+        nodes1[3] = nodes1[1];
+        nodes2[ip1] = nodes1[1]; ip2 = ip1;
+        other_nodes.push_back(0);
+        other_nodes.push_back(2);
+        break;
       case PRISM :
-       nodes1[4] = nodes1[0]; nodes1[5] = nodes1[1];
-       nodes2[ip1] = nodes1[0];
-       nodes2[ip2] = nodes1[1];
-       other_nodes.push_back(2);
-       other_nodes.push_back(6);
-       break;
+        nodes1[4] = nodes1[0]; nodes1[5] = nodes1[1];
+        nodes2[ip1] = nodes1[0];
+        nodes2[ip2] = nodes1[1];
+        other_nodes.push_back(2);
+        other_nodes.push_back(6);
+        break;
       case TETRA_CYL :
-       nodes3[0] = nodes1[0]; nodes3[1] = nodes1[1];
-       nodes3[2] = nodes1[2]; nodes3[3] = nodes1[4];
-       // nodes1[4] = nodes1[0]; nodes1[7] = base_node(1.0, 1.0, 2.0);
-       nodes2[ip1] = nodes1[1]; ip2 = ip1;
-       other_nodes.push_back(0);
-       other_nodes.push_back(2);
-       other_nodes.push_back(4);
-       break;
+        nodes3[0] = nodes1[0]; nodes3[1] = nodes1[1];
+        nodes3[2] = nodes1[2]; nodes3[3] = nodes1[4];
+        // nodes1[4] = nodes1[0]; nodes1[7] = base_node(1.0, 1.0, 2.0);
+        nodes2[ip1] = nodes1[1]; ip2 = ip1;
+        other_nodes.push_back(0);
+        other_nodes.push_back(2);
+        other_nodes.push_back(4);
+        break;
       case PRISM2 :
-       nodes2[ip1] = nodes1[4];
-       other_nodes.push_back(0);
-       other_nodes.push_back(1);
-       other_nodes.push_back(2);
-       break;
+        nodes2[ip1] = nodes1[4];
+        other_nodes.push_back(0);
+        other_nodes.push_back(1);
+        other_nodes.push_back(2);
+        break;
       case PYRAMID :
-       ip2 = ip1 = 0;
-       nodes1[0] =  base_node(-1.,-1., 0.);
-       nodes1[1] =  base_node( 1.,-1., 0.);
-       nodes1[2] =  base_node(-1., 1., 0.);
-       nodes1[3] =  base_node( 1., 1., 0.);
-       nodes1[4] =  base_node( 0., 0., 1.);
-       nodes1[5] =  nodes1[6] = nodes1[7] =  nodes1[4];
-       nodes2[ip1] = nodes1[0];
-       other_nodes.push_back(4);
-       other_nodes.push_back(3);
-       other_nodes.push_back(2);
-       other_nodes.push_back(1);
+        ip2 = ip1 = 0;
+        nodes1[0] =  base_node(-1.,-1., 0.);
+        nodes1[1] =  base_node( 1.,-1., 0.);
+        nodes1[2] =  base_node(-1., 1., 0.);
+        nodes1[3] =  base_node( 1., 1., 0.);
+        nodes1[4] =  base_node( 0., 0., 1.);
+        nodes1[5] =  nodes1[6] = nodes1[7] =  nodes1[4];
+        nodes2[ip1] = nodes1[0];
+        other_nodes.push_back(4);
+        other_nodes.push_back(3);
+        other_nodes.push_back(2);
+        other_nodes.push_back(1);
       }
 
       for (size_type i = 0; i <= nodes2.size()-1; ++i)
-       if (i != ip1 && i != ip2) {
-         GMM_ASSERT3(!other_nodes.empty(), "");
-         nodes2[i] = nodes1[other_nodes.back()];
-         other_nodes.pop_back();
-       }
+        if (i != ip1 && i != ip2) {
+          GMM_ASSERT3(!other_nodes.empty(), "");
+          nodes2[i] = nodes1[other_nodes.back()];
+          other_nodes.pop_back();
+        }
 
-      base_matrix G1, G2, G3; 
+      base_matrix G1, G2, G3;
       base_matrix K(N, N), grad(N, N), K3(N, N), K4(N, N);
       base_node normal1(N), normal2(N);
       bgeot::geotrans_inv_convex gic(nodes2, pgt2);
       scalar_type J1, J2, J3(1), J4(1);
-      
+
       bgeot::vectors_to_base_matrix(G1, nodes1);
       bgeot::vectors_to_base_matrix(G2, nodes2);
 
       for (size_type nc = 0; nc < 2; ++nc) {
-       
-       if (what == TETRA_CYL) {
-         if (nc == 1) nodes3[0] = nodes1[3];
-         bgeot::vectors_to_base_matrix(G3, nodes3);
-         pgt3->poly_vector_grad(nodes1[0], grad);
-         gmm::mult(gmm::transposed(grad), gmm::transposed(G3), K3);
-         J3 = gmm::abs(gmm::lu_inverse(K3)); /* = 1 */
-       }
-
-       for (size_type i=0; i <  base_im->nb_points(); ++i) {
-
-         gmm::copy(gmm::identity_matrix(), K4); J4 = J1 = scalar_type(1);
-
-         size_type fp = size_type(-1); // Search the face number in the
-         if (i >= base_im->nb_points_on_convex()) { // original element
-           size_type ii = i - base_im->nb_points_on_convex();
-           for (short_type ff=0; ff < base_im->structure()->nb_faces(); ++ff) {
-             if (ii < base_im->nb_points_on_face(ff)) { fp = ff; break; }
-             else ii -= base_im->nb_points_on_face(ff);
-           }
-           normal1 =  base_im->ref_convex()->normals()[fp];
-         }
-
-         base_node P = base_im->point(i);
-         if (what == TETRA_CYL) { 
-           P = pgt3->transform(P, nodes3);
-           scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
-           K4(0, 1) = - y / one_minus_z;
-           K4(1, 1) = 1.0 - x / one_minus_z;
-           K4(2, 1) = - x * y / gmm::sqr(one_minus_z);
-           J4 = gmm::abs(gmm::lu_det(K4));
-           P[1] *= (1.0 - x / one_minus_z);
-         }
-         if (what == PRISM2) {
-            scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
-            K4(0,0) = one_minus_z; K4(2,0) = -x;
-            K4(1,1) = one_minus_z; K4(2,1) = -y;
-            J4 = gmm::abs(gmm::lu_det(K4));
-            P[0] *= one_minus_z;
-            P[1] *= one_minus_z;
-         }
-
-         base_node P1 = pgt1->transform(P, nodes1), P2(N);
-         pgt1->poly_vector_grad(P, grad);
-         gmm::mult(gmm::transposed(grad), gmm::transposed(G1), K);
-         J1 = gmm::abs(gmm::lu_det(K)) * J3 * J4;
-
-         if (fp != size_type(-1) && J1 > 1E-10 && J4 > 1E-10) {
-           if (what == TETRA_CYL) {
-             gmm::mult(K3, normal1, normal2);
-             normal1 = normal2;
-           }
-           gmm::lu_inverse(K4);
-           gmm::lu_inverse(K);
-           gmm::mult(K4, normal1, normal2);
-           gmm::mult(K, normal2, normal1);
-           normal2 = normal1;
-           J1 *= gmm::vect_norm2(normal2);
-           normal2 /= gmm::vect_norm2(normal2);
-         }
-         gic.invert(P1, P2);
-         GMM_ASSERT1(pgt2->convex_ref()->is_in(P2) < 1E-8,
-                     "Point not in the convex ref : " << P2);
-         
-         pgt2->poly_vector_grad(P2, grad);
-         gmm::mult(gmm::transposed(grad), gmm::transposed(G2), K);
-         J2 = gmm::abs(gmm::lu_det(K)); /* = 1 */
-         
-         if (i <  base_im->nb_points_on_convex())
-           add_point(P2, base_im->coeff(i)*J1/J2, short_type(-1));
-         else if (J1 > 1E-10) {
-           short_type f = short_type(-1);
-           for (short_type ff = 0; ff <= N; ++ff)
-             if (gmm::abs(pgt2->convex_ref()->is_in_face(ff, P2)) < 1E-8) {
-               GMM_ASSERT1(f == short_type(-1),
-                           "An integration point is common to two faces");
-               f = ff;
-             }
-           if (f != short_type(-1)) {
-             gmm::mult(K, normal2, normal1);
-             add_point(P2,base_im->coeff(i)*J1*gmm::vect_norm2(normal1)/J2,f);
-           }
-           // else { cout << "Point " << P2 << " eliminated" << endl; }
-         }  
-       }
-       if (what != TETRA_CYL) break;
+
+        if (what == TETRA_CYL) {
+          if (nc == 1) nodes3[0] = nodes1[3];
+          bgeot::vectors_to_base_matrix(G3, nodes3);
+          pgt3->poly_vector_grad(nodes1[0], grad);
+          gmm::mult(gmm::transposed(grad), gmm::transposed(G3), K3);
+          J3 = gmm::abs(gmm::lu_inverse(K3)); /* = 1 */
+        }
+
+        for (size_type i=0; i <  base_im->nb_points(); ++i) {
+
+          gmm::copy(gmm::identity_matrix(), K4); J4 = J1 = scalar_type(1);
+
+          size_type fp = size_type(-1); // Search the face number in the
+          if (i >= base_im->nb_points_on_convex()) { // original element
+            size_type ii = i - base_im->nb_points_on_convex();
+            for (short_type ff=0; ff < base_im->structure()->nb_faces(); ++ff) 
{
+              if (ii < base_im->nb_points_on_face(ff)) { fp = ff; break; }
+              else ii -= base_im->nb_points_on_face(ff);
+            }
+            normal1 =  base_im->ref_convex()->normals()[fp];
+          }
+
+          base_node P = base_im->point(i);
+          if (what == TETRA_CYL) {
+            P = pgt3->transform(P, nodes3);
+            scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
+            K4(0, 1) = - y / one_minus_z;
+            K4(1, 1) = 1.0 - x / one_minus_z;
+            K4(2, 1) = - x * y / gmm::sqr(one_minus_z);
+            J4 = gmm::abs(gmm::lu_det(K4));
+            P[1] *= (1.0 - x / one_minus_z);
+          }
+          if (what == PRISM2) {
+             scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
+             K4(0,0) = one_minus_z; K4(2,0) = -x;
+             K4(1,1) = one_minus_z; K4(2,1) = -y;
+             J4 = gmm::abs(gmm::lu_det(K4));
+             P[0] *= one_minus_z;
+             P[1] *= one_minus_z;
+          }
+
+          base_node P1 = pgt1->transform(P, nodes1), P2(N);
+          pgt1->poly_vector_grad(P, grad);
+          gmm::mult(gmm::transposed(grad), gmm::transposed(G1), K);
+          J1 = gmm::abs(gmm::lu_det(K)) * J3 * J4;
+
+          if (fp != size_type(-1) && J1 > 1E-10 && J4 > 1E-10) {
+            if (what == TETRA_CYL) {
+              gmm::mult(K3, normal1, normal2);
+              normal1 = normal2;
+            }
+            gmm::lu_inverse(K4);
+            gmm::lu_inverse(K);
+            gmm::mult(K4, normal1, normal2);
+            gmm::mult(K, normal2, normal1);
+            normal2 = normal1;
+            J1 *= gmm::vect_norm2(normal2);
+            normal2 /= gmm::vect_norm2(normal2);
+          }
+          gic.invert(P1, P2);
+          GMM_ASSERT1(pgt2->convex_ref()->is_in(P2) < 1E-8,
+                      "Point not in the convex ref : " << P2);
+
+          pgt2->poly_vector_grad(P2, grad);
+          gmm::mult(gmm::transposed(grad), gmm::transposed(G2), K);
+          J2 = gmm::abs(gmm::lu_det(K)); /* = 1 */
+
+          if (i <  base_im->nb_points_on_convex())
+            add_point(P2, base_im->coeff(i)*J1/J2, short_type(-1));
+          else if (J1 > 1E-10) {
+            short_type f = short_type(-1);
+            for (short_type ff = 0; ff <= N; ++ff)
+              if (gmm::abs(pgt2->convex_ref()->is_in_face(ff, P2)) < 1E-8) {
+                GMM_ASSERT1(f == short_type(-1),
+                            "An integration point is common to two faces");
+                f = ff;
+              }
+            if (f != short_type(-1)) {
+              gmm::mult(K, normal2, normal1);
+              add_point(P2,base_im->coeff(i)*J1*gmm::vect_norm2(normal1)/J2,f);
+            }
+            // else { cout << "Point " << P2 << " eliminated" << endl; }
+          }
+        }
+        if (what != TETRA_CYL) break;
       }
       valid_method();
     }
   };
 
 
-  static pintegration_method quasi_polar(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  quasi_polar(im_param_list &params,
+              std::vector<dal::pstatic_stored_object> &dependencies) {
     GMM_ASSERT1(params.size() >= 1 && params[0].type() == 1,
-               "Bad parameters for quasi polar integration: the first "
-               "parameter should be an integration method");
+                "Bad parameters for quasi polar integration: the first "
+                "parameter should be an integration method");
     pintegration_method a = params[0].method();
     GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
     int ip1 = 0, ip2 = 0;
@@ -1027,31 +1040,32 @@ namespace getfem {
       GMM_ASSERT1(params.size() == 1, "Bad number of parameters");
     } else {
       GMM_ASSERT1(params.size() == 2 || params.size() == 3,
-                 "Bad number of parameters : " << params.size()
-                 << " should be 2 or 3.");
+                  "Bad number of parameters : " << params.size()
+                  << " should be 2 or 3.");
       GMM_ASSERT1(params[1].type() == 0
-                 && params.back().type() == 0, "Bad type of parameters");
+                  && params.back().type() == 0, "Bad type of parameters");
       ip1 = int(::floor(params[1].num() + 0.01));
       ip2 = int(::floor(params.back().num() + 0.01));
     }
     int N = a->approx_method()->dim();
     GMM_ASSERT1(N >= 2 && N <= 3 && ip1 >= 0 && ip2 >= 0 && ip1 <= N
-               && ip2 <= N, "Bad parameters");
+                && ip2 <= N, "Bad parameters");
 
     papprox_integration
       pai = std::make_shared<quasi_polar_integration>(a->approx_method(),
-                                                     ip1, ip2);
+                                                      ip1, ip2);
     pintegration_method p = std::make_shared<integration_method>(pai);
     dependencies.push_back(p->approx_method()->ref_convex());
     dependencies.push_back(p->approx_method()->pintegration_points());
     return p;
   }
 
-  static pintegration_method pyramid(im_param_list &params,
-       std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pintegration_method
+  pyramid(im_param_list &params,
+          std::vector<dal::pstatic_stored_object> &dependencies) {
     GMM_ASSERT1(params.size() == 1 && params[0].type() == 1,
-               "Bad parameters for pyramid integration: the first "
-               "parameter should be an integration method");
+                "Bad parameters for pyramid integration: the first "
+                "parameter should be an integration method");
     pintegration_method a = params[0].method();
     GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
     int N = a->approx_method()->dim();
@@ -1071,8 +1085,9 @@ namespace getfem {
   /*    Naming system                                                     */
   /* ******************************************************************** */
 
-  pintegration_method structured_composite_int_method(im_param_list &,
-                     std::vector<dal::pstatic_stored_object> &);
+  pintegration_method
+  structured_composite_int_method(im_param_list &,
+                                  std::vector<dal::pstatic_stored_object> &);
   pintegration_method HCT_composite_int_method(im_param_list &params,
    std::vector<dal::pstatic_stored_object> &dependencies);
 
@@ -1109,7 +1124,7 @@ namespace getfem {
   };
 
   pintegration_method int_method_descriptor(std::string name,
-                                           bool throw_if_not_found) {
+                                            bool throw_if_not_found) {
     size_type i = 0;
     return dal::singleton<im_naming_system>::instance().method
       (name, i, throw_if_not_found);
@@ -1122,13 +1137,13 @@ namespace getfem {
 
   // allows the add of an integration method.
   void add_integration_name(std::string name,
-                       dal::naming_system<integration_method>::pfunction f) {
+                        dal::naming_system<integration_method>::pfunction f) {
     dal::singleton<im_naming_system>::instance().add_suffix(name, f);
   }
 
 
   /* Fonctions pour la ref. directe.                                     */
-  
+
   pintegration_method exact_simplex_im(size_type n) {
     DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pintegration_method, pim, 0);
     DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(size_type, d, -2);
@@ -1186,20 +1201,20 @@ namespace getfem {
 
     if (nbp == n+1)
       if (cvs == bgeot::simplex_structure(dim_type(n)))
-       { name << "IM_EXACT_SIMPLEX("; found = true; }
-    
+        { name << "IM_EXACT_SIMPLEX("; found = true; }
+
     /* Identifying Q1-parallelepiped.                                     */
 
     if (!found && nbp == (size_type(1) << n))
       if (cvs == bgeot::parallelepiped_structure(dim_type(n)))
-       { name << "IM_EXACT_PARALLELEPIPED("; found = true; }
+        { name << "IM_EXACT_PARALLELEPIPED("; found = true; }
 
     /* Identifying Q1-prisms.                                             */
- 
+
     if (!found && nbp == 2 * n)
       if (cvs == bgeot::prism_structure(dim_type(n)))
-       { name << "IM_EXACT_PRISM("; found = true; }
-     
+        { name << "IM_EXACT_PRISM("; found = true; }
+
     // To be completed
 
     if (found) {
@@ -1208,7 +1223,7 @@ namespace getfem {
       cvs_last = cvs;
       return im_last;
     }
- 
+
     GMM_ASSERT1(false, "This element is not taken into account. Contact us");
   }
 
@@ -1247,7 +1262,7 @@ namespace getfem {
                   << " of degree >= " << int(degree));
     } else if (cvs->is_product(&a,&b) ||
                (bgeot::basic_structure(cvs).get() && 
bgeot::basic_structure(cvs)->is_product(&a,&b))) {
-      name << "IM_PRODUCT(" 
+      name << "IM_PRODUCT("
            << name_of_int_method(classical_approx_im_(a,degree)) << ","
            << name_of_int_method(classical_approx_im_(b,degree)) << ")";
     } else GMM_ASSERT1(false, "unknown convex structure!");
@@ -1255,7 +1270,7 @@ namespace getfem {
   }
 
   pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt,
-                                         dim_type degree) {
+                                          dim_type degree) {
     DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(bgeot::pgeometric_trans, pgt_last, 
0);
     DEFINE_STATIC_THREAD_LOCAL(dim_type, degree_last);
     DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pintegration_method, im_last, 0);
@@ -1267,13 +1282,13 @@ namespace getfem {
     return im_last;
   }
 
-  pintegration_method im_none(void) { 
+  pintegration_method im_none() {
      DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pintegration_method,im_last,0);
      if (!im_last) im_last = int_method_descriptor("IM_NONE");
      return im_last;
   }
 
-  /* try to integrate all monomials up to order 'order' and return the 
+  /* try to integrate all monomials up to order 'order' and return the
      max. error */
   scalar_type test_integration_error(papprox_integration pim, dim_type order) {
     short_type dim = pim->dim();
@@ -1282,10 +1297,10 @@ namespace getfem {
     for (bgeot::power_index idx(dim); idx.degree() <= order; ++idx) {
       opt_long_scalar_type sum(0), realsum;
       for (size_type i=0; i < pim->nb_points_on_convex(); ++i) {
-       opt_long_scalar_type prod = pim->coeff(i);
-       for (size_type d=0; d < dim; ++d) 
-         prod *= pow(opt_long_scalar_type(pim->point(i)[d]), idx[d]);
-       sum += prod;
+        opt_long_scalar_type prod = pim->coeff(i);
+        for (size_type d=0; d < dim; ++d)
+          prod *= pow(opt_long_scalar_type(pim->point(i)[d]), idx[d]);
+        sum += prod;
       }
       realsum = exact->exact_method()->int_monomial(idx);
       error = std::max(error, gmm::abs(realsum-sum));
@@ -1295,7 +1310,7 @@ namespace getfem {
 
   papprox_integration get_approx_im_or_fail(pintegration_method pim) {
     GMM_ASSERT1(pim->type() == IM_APPROX,  "error estimate work only with "
-               "approximate integration methods");
+                "approximate integration methods");
     return pim->approx_method();
   }
 
diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc
index caeb31d..84a21d5 100644
--- a/src/getfem_mesh.cc
+++ b/src/getfem_mesh.cc
@@ -299,8 +299,8 @@ namespace getfem {
   }
 
   size_type mesh::add_pyramid(size_type a, size_type b,
-                                  size_type c, size_type d, size_type e) {
-    size_type ipt[5]; ipt[0] = a; ipt[1] = b; ipt[2] = c; ipt[3] = d; ipt[4] = 
e;
+                              size_type c, size_type d, size_type e) {
+    size_type ipt[5] = {a, b, c, d, e};
     return add_convex(bgeot::pyramidal_geotrans(1), &(ipt[0]));
   }
 
diff --git a/src/getfem_mesh_fem.cc b/src/getfem_mesh_fem.cc
index c3b4927..0816453 100644
--- a/src/getfem_mesh_fem.cc
+++ b/src/getfem_mesh_fem.cc
@@ -143,20 +143,20 @@ namespace getfem {
                   " and target_dim " << int(pf->target_dim()) << " of " <<
                   name_of_fem(pf));
 
-      
+
       if (cv == f_elems.size()) {
-       f_elems.push_back(pf);
-       fe_convex.add(cv);
-       dof_enumeration_made = false;
+        f_elems.push_back(pf);
+        fe_convex.add(cv);
+        dof_enumeration_made = false;
         touch(); v_num = act_counter();
       } else {
-       if (cv > f_elems.size()) f_elems.resize(cv+1);
-       if (!fe_convex.is_in(cv) || f_elems[cv] != pf) {
-         fe_convex.add(cv);
-         f_elems[cv] = pf;
-         dof_enumeration_made = false;
-         touch(); v_num = act_counter();
-       }
+        if (cv > f_elems.size()) f_elems.resize(cv+1);
+        if (!fe_convex.is_in(cv) || f_elems[cv] != pf) {
+          fe_convex.add(cv);
+          f_elems[cv] = pf;
+          dof_enumeration_made = false;
+          touch(); v_num = act_counter();
+        }
       }
     }
   }
@@ -335,7 +335,7 @@ namespace getfem {
     // Information for global dof
     dal::bit_vector encountered_global_dof, processed_elt;
     dal::dynamic_array<size_type> ind_global_dof;
-    
+
     // Auxilliary variables
     std::vector<size_type> itab;
     base_node P(linked_mesh().dim());
@@ -349,25 +349,25 @@ namespace getfem {
     bgeot::pgeotrans_precomp pgp = 0;
 
     for (dal::bv_visitor cv(linked_mesh().convex_index());
-        !cv.finished(); ++cv) {
+         !cv.finished(); ++cv) {
       if (fe_convex.is_in(cv)) {
-       gmm::copy(linked_mesh().points_of_convex(cv)[0], bmin);
-       gmm::copy(bmin, bmax);
-       for (size_type i = 0; i <  linked_mesh().nb_points_of_convex(cv); ++i) {
-         const base_node &pt = linked_mesh().points_of_convex(cv)[i];
-         for (size_type d = 1; d < bmin.size(); ++d) {
-           bmin[d] = std::min(bmin[d], pt[d]);
-           bmax[d] = std::max(bmax[d], pt[d]);
-         }
-       }       
-       elt_car_sizes[cv] = gmm::vect_dist2_sqr(bmin, bmax);
+        gmm::copy(linked_mesh().points_of_convex(cv)[0], bmin);
+        gmm::copy(bmin, bmax);
+        for (size_type i = 0; i < linked_mesh().nb_points_of_convex(cv); ++i) {
+          const base_node &pt = linked_mesh().points_of_convex(cv)[i];
+          for (size_type d = 1; d < bmin.size(); ++d) {
+            bmin[d] = std::min(bmin[d], pt[d]);
+            bmax[d] = std::max(bmax[d], pt[d]);
+          }
+        }
+        elt_car_sizes[cv] = gmm::vect_dist2_sqr(bmin, bmax);
       }
     }
 
     dal::bit_vector cv_done;
-    
+
     for (dal::bv_visitor cv(linked_mesh().convex_index());
-        !cv.finished(); ++cv) { // Loop on elements
+         !cv.finished(); ++cv) { // Loop on elements
       if (!fe_convex.is_in(cv)) continue;
       pfem pf = fem_of_element(cv);
       if (pf != first_pf) is_uniform_ = false;
@@ -400,34 +400,34 @@ namespace getfem {
           pgp->transform(linked_mesh().points_of_convex(cv), i, P);
           size_type idof = nbdof;
 
-         if (dof_nodes[cv].nb_points() > 0) {
-           scalar_type dist = dof_nodes[cv].nearest_neighbor(ipt, P);
-           if (gmm::abs(dist) <= 1e-6*elt_car_sizes[cv]) {
-             fd.ind_node=ipt.i;
-             auto it = dof_sorts[cv].find(fd);
-             if (it != dof_sorts[cv].end()) idof = it->second;
-           }
-         }
-         
+          if (dof_nodes[cv].nb_points() > 0) {
+            scalar_type dist = dof_nodes[cv].nearest_neighbor(ipt, P);
+            if (gmm::abs(dist) <= 1e-6*elt_car_sizes[cv]) {
+              fd.ind_node=ipt.i;
+              auto it = dof_sorts[cv].find(fd);
+              if (it != dof_sorts[cv].end()) idof = it->second;
+            }
+          }
+
           if (idof == nbdof) {
-           nbdof += Qdim / pf->target_dim();
-
-           linked_mesh().neighbours_of_convex(cv, pf->faces_of_dof(cv, i), s);
-           for (size_type ncv : s) { // For each unscanned neighbour
-             if (!cv_done[ncv] && fe_convex.is_in(ncv)) { // add the dof
-
-               fd.ind_node = size_type(-1);
-               if (dof_nodes[ncv].nb_points() > 0) {
-                 scalar_type dist = dof_nodes[ncv].nearest_neighbor(ipt, P);
-                 if (gmm::abs(dist) <= 1e-6*elt_car_sizes[ncv])
-                   fd.ind_node=ipt.i;
-               }
-               if (fd.ind_node == size_type(-1))
-                 fd.ind_node = dof_nodes[ncv].add_point(P);
-               dof_sorts[ncv][fd] = idof;
-             }
-           }
-         }
+            nbdof += Qdim / pf->target_dim();
+
+            linked_mesh().neighbours_of_convex(cv, pf->faces_of_dof(cv, i), s);
+            for (size_type ncv : s) { // For each unscanned neighbour
+              if (!cv_done[ncv] && fe_convex.is_in(ncv)) { // add the dof
+
+                fd.ind_node = size_type(-1);
+                if (dof_nodes[ncv].nb_points() > 0) {
+                  scalar_type dist = dof_nodes[ncv].nearest_neighbor(ipt, P);
+                  if (gmm::abs(dist) <= 1e-6*elt_car_sizes[ncv])
+                    fd.ind_node=ipt.i;
+                }
+                if (fd.ind_node == size_type(-1))
+                  fd.ind_node = dof_nodes[ncv].add_point(P);
+                dof_sorts[ncv][fd] = idof;
+              }
+            }
+          }
           itab[i] = idof;
         }
       }
@@ -571,8 +571,8 @@ namespace getfem {
         } else if (bgeot::casecmp(tmp, "DOF_ENUMERATION") == 0) {
           dal::bit_vector doflst;
           dof_structure.clear(); dof_enumeration_made = false;
-         is_uniform_ = true;
-         size_type nbdof_unif = size_type(-1);
+          is_uniform_ = true;
+          size_type nbdof_unif = size_type(-1);
           touch(); v_num = act_counter();
           while (true) {
             bgeot::get_token(ist, tmp);
@@ -585,11 +585,11 @@ namespace getfem {
             std::vector<size_type> tab;
             if (convex_index().is_in(ic) && tmp.size() &&
                 isdigit(tmp[0]) && tmp2 == ":") {
-             size_type nbd = nb_basic_dof_of_element(ic);
-             if (nbdof_unif == size_type(-1))
-               nbdof_unif = nbd;
-             else if (nbdof_unif != nbd)
-               is_uniform_ = false;
+              size_type nbd = nb_basic_dof_of_element(ic);
+              if (nbdof_unif == size_type(-1))
+                nbdof_unif = nbd;
+              else if (nbdof_unif != nbd)
+                is_uniform_ = false;
               tab.resize(nb_basic_dof_of_element(ic));
               for (size_type i=0; i < fem_of_element(ic)->nb_dof(ic);
                    i++) {

Reply via email to