getfem-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Getfem-commits] (no subject)


From: Konstantinos Poulios
Subject: [Getfem-commits] (no subject)
Date: Mon, 11 Jun 2018 09:03:36 -0400 (EDT)

branch: master
commit 8d65fe72becd9e819b23c4637a5b561f2323cc46
Author: Konstantinos Poulios <address@hidden>
Date:   Mon Jun 11 15:03:03 2018 +0200

    White space, file enconding and other minor changes
---
 src/bgeot_mesh_structure.cc |   2 +-
 src/bgeot_poly_composite.cc | 198 ++++++++++++------------
 src/getfem/bgeot_poly.h     | 368 ++++++++++++++++++++++----------------------
 src/getfem_export.cc        |   3 +-
 src/getfem_fem.cc           |  34 ++--
 src/getfem_import.cc        |  14 +-
 src/getfem_mesh_fem.cc      |  35 ++---
 7 files changed, 331 insertions(+), 323 deletions(-)

diff --git a/src/bgeot_mesh_structure.cc b/src/bgeot_mesh_structure.cc
index c1feaa0..721c3b4 100644
--- a/src/bgeot_mesh_structure.cc
+++ b/src/bgeot_mesh_structure.cc
@@ -387,7 +387,7 @@ namespace bgeot {
         ncs--;
         cvs = cvstab[ncs];
         std::vector< size_type > pts = indpttab[ncs];
-        if (cvs->dim() == 1) { // il faudrait �tendre aux autres cas 
classiques.
+        if (cvs->dim() == 1) { // il faudrait étendre aux autres cas 
classiques.
 
           for (size_type j = 1; j < cvs->nb_points(); ++j) {
             //cerr << "ncs=" << ncs << "j=" << j << ", ajout de " << 
(indpttab[ncs])[j-1] << "," << (indpttab[ncs])[j] << endl;
diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index 5a2fb8c..888788e 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -31,23 +31,23 @@ namespace bgeot {
 
   int imbricated_box_less::operator()(const base_node &x,
                                       const base_node &y) const {
-    size_type s = x.size(); 
+    size_type s = x.size();
     scalar_type c1 = c_max, c2 = c_max * scalar_type(base);
     GMM_ASSERT2(y.size() == s, "dimension error");
-    
+
     base_node::const_iterator itx=x.begin(), itex=x.end(), ity=y.begin();
     int ret = 0;
     for (; itx != itex; ++itx, ++ity) {
       long a = long(sfloor((*itx) * c1)), b = long(sfloor((*ity) * c1));
       if ((scalar_type(gmm::abs(a)) > scalar_type(base))
-          || (scalar_type(gmm::abs(b)) > scalar_type(base))) { 
+          || (scalar_type(gmm::abs(b)) > scalar_type(base))) {
         exp_max++; c_max /= scalar_type(base);
         return (*this)(x,y);
       }
       if (ret == 0) { if (a < b) ret = -1; else if (a > b) ret = 1; }
     }
     if (ret) return ret;
-    
+
     for (int e = exp_max; e >= exp_min; --e, c1 *= scalar_type(base),
            c2 *= scalar_type(base)) {
       itx = x.begin(), itex = x.end(), ity = y.begin();
@@ -99,23 +99,23 @@ namespace bgeot {
     mesh_structure::ind_cv_ct::const_iterator itc, itce;
 
     mesh_precomposite::PTAB::const_sorted_iterator
-      it1 = mp->vertexes.sorted_ge(pt), it2 = it1;    
+      it1 = mp->vertexes.sorted_ge(pt), it2 = it1;
     size_type i1 = it1.index(), i2;
 
     --it2; i2 = it2.index();
 
 
-    while (i1 != size_type(-1) || i2 != size_type(-1)) 
+    while (i1 != size_type(-1) || i2 != size_type(-1))
     {
-      if (i1 != size_type(-1)) 
+      if (i1 != size_type(-1))
       {
         const mesh_structure::ind_cv_ct &tc
           = mp->linked_mesh().convex_to_point(i1);
         itc = tc.begin(); itce = tc.end();
-        for (; itc != itce; ++itc) 
+        for (; itc != itce; ++itc)
         {
           size_type ii = *itc;
-          if (elt[ii]) 
+          if (elt[ii])
           {
             elt[ii] = false;
             p0 = pt; p0 -= mp->orgs[ii];
@@ -128,15 +128,15 @@ namespace bgeot {
         ++it1; i1 = it1.index();
       }
 
-      if (i2 != size_type(-1)) 
+      if (i2 != size_type(-1))
       {
         const mesh_structure::ind_cv_ct &tc
           = mp->linked_mesh().convex_to_point(i2);
         itc = tc.begin(); itce = tc.end();
-        for (; itc != itce; ++itc) 
+        for (; itc != itce; ++itc)
         {
           size_type ii = *itc;
-          if (elt[ii]) 
+          if (elt[ii])
           {
             elt[ii] = false;
             p0 = pt; p0 -= mp->orgs[ii];
@@ -151,7 +151,7 @@ namespace bgeot {
     }
     GMM_ASSERT1(false, "Element not found in composite polynomial: " << pt);
   }
-  
+
   DAL_TRIPLE_KEY(base_poly_key, short_type, short_type, 
std::vector<opt_long_scalar_type>);
 
   polynomial_composite::polynomial_composite(
@@ -204,10 +204,11 @@ namespace bgeot {
   /* build a regularly refined mesh for a simplex of dimension <= 3.
   All simplexes have the same orientation as the original simplex.
   */
-  static void structured_mesh_for_simplex_(pconvex_structure cvs, 
-    pgeometric_trans opt_gt,
-    const std::vector<base_node> *opt_gt_pts,
-    short_type k, pbasic_mesh pm) {
+  static void
+  structured_mesh_for_simplex_(pconvex_structure cvs,
+                               pgeometric_trans opt_gt,
+                               const std::vector<base_node> *opt_gt_pts,
+                               short_type k, pbasic_mesh pm) {
       scalar_type h = 1./k;
       switch (cvs->dim()) {
       case 1 :
@@ -234,11 +235,11 @@ namespace bgeot {
               B[0] = b; B[1] = c;
               C[0] = a; C[1] = d;
               D[0] = b; D[1] = d;
-              if (opt_gt) { 
-                A = opt_gt->transform(A, *opt_gt_pts); 
-                B = opt_gt->transform(B, *opt_gt_pts); 
-                C = opt_gt->transform(C, *opt_gt_pts); 
-                D = opt_gt->transform(D, *opt_gt_pts); 
+              if (opt_gt) {
+                A = opt_gt->transform(A, *opt_gt_pts);
+                B = opt_gt->transform(B, *opt_gt_pts);
+                C = opt_gt->transform(C, *opt_gt_pts);
+                D = opt_gt->transform(D, *opt_gt_pts);
               }
               /* add triangle with correct orientation (det [B-A;C-A] > 0) */
               size_type nA = pm->add_point(A);
@@ -252,9 +253,9 @@ namespace bgeot {
           }
         }
         break;
-      case 3 : 
+      case 3 :
         {
-          /* based on decompositions of small cubes 
+          /* based on decompositions of small cubes
           the number of tetrahedrons is k^3
           */
           base_node A,B,C,D,E,F,G,H;
@@ -273,15 +274,15 @@ namespace bgeot {
                 F = {x+h, y, z+h};
                 G = {x, y+h, z+h};
                 H = {x+h, y+h, z+h};
-                if (opt_gt) { 
-                  A = opt_gt->transform(A, *opt_gt_pts); 
-                  B = opt_gt->transform(B, *opt_gt_pts); 
-                  C = opt_gt->transform(C, *opt_gt_pts); 
-                  D = opt_gt->transform(D, *opt_gt_pts); 
-                  E = opt_gt->transform(E, *opt_gt_pts); 
-                  F = opt_gt->transform(F, *opt_gt_pts); 
-                  G = opt_gt->transform(G, *opt_gt_pts); 
-                  H = opt_gt->transform(H, *opt_gt_pts); 
+                if (opt_gt) {
+                  A = opt_gt->transform(A, *opt_gt_pts);
+                  B = opt_gt->transform(B, *opt_gt_pts);
+                  C = opt_gt->transform(C, *opt_gt_pts);
+                  D = opt_gt->transform(D, *opt_gt_pts);
+                  E = opt_gt->transform(E, *opt_gt_pts);
+                  F = opt_gt->transform(F, *opt_gt_pts);
+                  G = opt_gt->transform(G, *opt_gt_pts);
+                  H = opt_gt->transform(H, *opt_gt_pts);
                 }
                 size_type t[8];
                 t[0] = pm->add_point(A);
@@ -316,7 +317,7 @@ namespace bgeot {
           }
         }
         break;
-      default : 
+      default :
         GMM_ASSERT1(false, "Sorry, not implemented for simplices of "
           "dimension " << int(cvs->dim()));
       }
@@ -347,7 +348,7 @@ namespace bgeot {
         { kcnt[kk]++; if (kcnt[kk] == k+1) { kcnt[kk] = 0; kk++; } else break; 
}
       }
 
-      /* 
+      /*
       insert convexes using node ids stored in 'pids'
       */
       std::vector<size_type> ppts(pow2n);
@@ -367,22 +368,23 @@ namespace bgeot {
       }
   }
 
-  static void structured_mesh_for_convex_(pconvex_structure cvs, 
-    pgeometric_trans opt_gt,
-    const std::vector<base_node> *opt_gt_pts,
-    short_type k, pbasic_mesh pm) {
+  static void
+  structured_mesh_for_convex_(pconvex_structure cvs,
+                              pgeometric_trans opt_gt,
+                              const std::vector<base_node> *opt_gt_pts,
+                              short_type k, pbasic_mesh pm) {
       size_type nbp = basic_structure(cvs)->nb_points();
       dim_type n = cvs->dim();
-      /* Identifying simplexes.                                           */   
 
-      if (nbp == size_type(n+1) && 
+      /* Identifying simplexes.                                           */
+      if (nbp == size_type(n+1) &&
           basic_structure(cvs)==simplex_structure(n)) {
           // smc.pm->write_to_file(cout);
           structured_mesh_for_simplex_(cvs,opt_gt,opt_gt_pts,k,pm);
           /* Identifying parallelepipeds.                                     
*/
-      } else if (nbp == (size_type(1) << n) && 
+      } else if (nbp == (size_type(1) << n) &&
           basic_structure(cvs) == parallelepiped_structure(n)) {
           structured_mesh_for_parallelepiped_(cvs,opt_gt,opt_gt_pts,k,pm);
-      } else if (nbp == size_type(2 * n) && 
+      } else if (nbp == size_type(2 * n) &&
           basic_structure(cvs) == prism_structure(n)) {
           GMM_ASSERT1(false, "Sorry, structured_mesh not implemented for 
prisms.");
       } else {
@@ -434,7 +436,7 @@ namespace bgeot {
     std::vector<std::unique_ptr<mesh_structure>> pfacem; /* array of 
mesh_structures for faces */
     dal::bit_vector nodes_on_edges;
     std::unique_ptr<mesh_precomposite> pmp;
-    str_mesh_cv__(pconvex_structure c, short_type k, bool smesh_) : 
+    str_mesh_cv__(pconvex_structure c, short_type k, bool smesh_) :
       cvs(c), n(k), simplex_mesh(smesh_)
     { DAL_STORED_OBJECT_DEBUG_CREATED(this, "cv mesh"); }
     ~str_mesh_cv__() { DAL_STORED_OBJECT_DEBUG_DESTROYED(this,"cv mesh"); }
@@ -450,66 +452,66 @@ namespace bgeot {
   * TODO: move it into another file and separate the pmesh_precomposite part ?
   **/
   void structured_mesh_for_convex(pconvex_ref cvr, short_type k,
-    pbasic_mesh &pm, pmesh_precomposite &pmp, 
-    bool force_simplexification) {
-      size_type n = cvr->structure()->dim();
-      size_type nbp = basic_structure(cvr->structure())->nb_points();
-
-      force_simplexification = force_simplexification || (nbp == n+1);
-      dal::pstatic_stored_object_key
-        pk = std::make_shared<str_mesh_key>(basic_structure(cvr->structure()),
-                                            k, force_simplexification);
-
-      dal::pstatic_stored_object o = dal::search_stored_object(pk);
-      pstr_mesh_cv__ psmc;
-      if (o)
-        psmc = std::dynamic_pointer_cast<const str_mesh_cv__>(o);
-      else {
-
-        auto ppsmc = std::make_shared<str_mesh_cv__>
-          (basic_structure(cvr->structure()), k, force_simplexification);
-        str_mesh_cv__ &smc(*ppsmc);
-        psmc = ppsmc;
-
-        smc.pm = std::make_unique<basic_mesh>();
-
-        if (force_simplexification) {
-          // cout << "cvr = " << int(cvr->structure()->dim()) << " : "
-          //      << cvr->structure()->nb_points() << endl;
-          const mesh_structure* splx_mesh
-            = basic_convex_ref(cvr)->simplexified_convex();
-          /* splx_mesh is correctly oriented */
-          for (size_type ic=0; ic < splx_mesh->nb_convex(); ++ic) {
-            std::vector<base_node> cvpts(splx_mesh->nb_points_of_convex(ic));
-            pgeometric_trans sgt
-              = simplex_geotrans(cvr->structure()->dim(), 1);
-            for (size_type j=0; j < cvpts.size(); ++j) {
-              cvpts[j] = basic_convex_ref(cvr)->points()
-                [splx_mesh->ind_points_of_convex(ic)[j]];
-              //cerr << "cvpts[" << j << "]=" << cvpts[j] << endl;
-            }
-            structured_mesh_for_convex_(splx_mesh->structure_of_convex(ic),
-                                        sgt, &cvpts, k, smc.pm.get());
+                                  pbasic_mesh &pm, pmesh_precomposite &pmp,
+                                  bool force_simplexification) {
+    size_type n = cvr->structure()->dim();
+    size_type nbp = basic_structure(cvr->structure())->nb_points();
+
+    force_simplexification = force_simplexification || (nbp == n+1);
+    dal::pstatic_stored_object_key
+      pk = std::make_shared<str_mesh_key>(basic_structure(cvr->structure()),
+                                          k, force_simplexification);
+
+    dal::pstatic_stored_object o = dal::search_stored_object(pk);
+    pstr_mesh_cv__ psmc;
+    if (o)
+      psmc = std::dynamic_pointer_cast<const str_mesh_cv__>(o);
+    else {
+
+      auto ppsmc = std::make_shared<str_mesh_cv__>
+        (basic_structure(cvr->structure()), k, force_simplexification);
+      str_mesh_cv__ &smc(*ppsmc);
+      psmc = ppsmc;
+
+      smc.pm = std::make_unique<basic_mesh>();
+
+      if (force_simplexification) {
+        // cout << "cvr = " << int(cvr->structure()->dim()) << " : "
+        //      << cvr->structure()->nb_points() << endl;
+        const mesh_structure* splx_mesh
+          = basic_convex_ref(cvr)->simplexified_convex();
+        /* splx_mesh is correctly oriented */
+        for (size_type ic=0; ic < splx_mesh->nb_convex(); ++ic) {
+          std::vector<base_node> cvpts(splx_mesh->nb_points_of_convex(ic));
+          pgeometric_trans sgt
+            = simplex_geotrans(cvr->structure()->dim(), 1);
+          for (size_type j=0; j < cvpts.size(); ++j) {
+            cvpts[j] = basic_convex_ref(cvr)->points()
+            [splx_mesh->ind_points_of_convex(ic)[j]];
+            //cerr << "cvpts[" << j << "]=" << cvpts[j] << endl;
           }
-        } else {
-          structured_mesh_for_convex_(cvr->structure(), 0, 0, k, smc.pm.get());
-        }
-        smc.pfacem.resize(cvr->structure()->nb_faces());
-        for (dim_type f=0; f < cvr->structure()->nb_faces(); ++f) {
-          smc.pfacem[f] = std::make_unique<mesh_structure>();
-          structured_mesh_of_faces_(cvr, f, *(smc.pm), *(smc.pfacem[f]));
+          structured_mesh_for_convex_(splx_mesh->structure_of_convex(ic),
+                                      sgt, &cvpts, k, smc.pm.get());
         }
-
-        smc.pmp = std::make_unique<mesh_precomposite>(*(smc.pm));
-        dal::add_stored_object(pk, psmc, basic_structure(cvr->structure()));
+      } else {
+        structured_mesh_for_convex_(cvr->structure(), 0, 0, k, smc.pm.get());
       }
-      pm  = psmc->pm.get();
-      pmp = psmc->pmp.get();
+      smc.pfacem.resize(cvr->structure()->nb_faces());
+      for (dim_type f=0; f < cvr->structure()->nb_faces(); ++f) {
+        smc.pfacem[f] = std::make_unique<mesh_structure>();
+        structured_mesh_of_faces_(cvr, f, *(smc.pm), *(smc.pfacem[f]));
+      }
+
+      smc.pmp = std::make_unique<mesh_precomposite>(*(smc.pm));
+      dal::add_stored_object(pk, psmc, basic_structure(cvr->structure()));
+    }
+    pm  = psmc->pm.get();
+    pmp = psmc->pmp.get();
   }
 
   const basic_mesh *
     refined_simplex_mesh_for_convex(pconvex_ref cvr, short_type k) {
-      pbasic_mesh pm; pmesh_precomposite pmp; 
+      pbasic_mesh pm; pmesh_precomposite pmp;
       structured_mesh_for_convex(cvr, k, pm, pmp, true);
       return pm;
   }
@@ -523,7 +525,7 @@ namespace bgeot {
     if (o) {
       pstr_mesh_cv__ psmc = std::dynamic_pointer_cast<const str_mesh_cv__>(o);
       return psmc->pfacem;
-    } 
+    }
     else GMM_ASSERT1(false,
                      "call refined_simplex_mesh_for_convex first (or fix me)");
   }
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index a0dc167..d7559e6 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -51,13 +51,13 @@ namespace bgeot
   /// used as the common short type integer in the library
   typedef gmm::uint16_type short_type;
   ///
-  
+
   /** Return the value of @f$ \frac{(n+p)!}{n!p!} @f$ which
    * is the number of monomials of a polynomial of @address@hidden
    * variables and degree @address@hidden
    */
   size_type alpha(short_type n, short_type d);
-  
+
   /** Vector of integer (16 bits type) which represent the powers
    *  of a monomial
    */
@@ -86,31 +86,31 @@ namespace bgeot
     const_reverse_iterator rend() const { return v.rend(); }
 
     size_type size() const { return v.size(); }
-    /// Gives the next power index 
+    /// Gives the next power index
     const power_index &operator ++();
-    /// Gives the next power index 
+    /// Gives the next power index
     const power_index operator ++(int)
-      { power_index res = *this; ++(*this); return res; } 
-    /// Gives the next previous index 
+      { power_index res = *this; ++(*this); return res; }
+    /// Gives the next previous index
     const power_index &operator --();
-    /// Gives the next previous index 
+    /// Gives the next previous index
     const power_index operator --(int)
       { power_index res = *this; --(*this); return res; }
     /**  Gives the global number of the index (i.e. the position of
      *   the corresponding monomial
      */
-    size_type global_index(void) const;
+    size_type global_index() const;
     /// Gives the degree.
-    short_type degree(void) const; 
+    short_type degree() const;
     /// Constructor
     power_index(short_type nn);
     /// Constructor
-    power_index(void) { dirty(); }
+    power_index() { dirty(); }
   };
-  
+
   /**
    * This class deals with plain polynomials with
-   * several variables. 
+   * several variables.
    *
    * A polynomial of @address@hidden variables and degree @address@hidden is 
stored in a vector
    * of @address@hidden components.
@@ -123,9 +123,9 @@ namespace bgeot
    *   bgeot::polynomial<double> P, Q;
    *   P = bgeot::polynomial<double>(2,2,1); // P = x
    *   Q = bgeot::polynomial<double>(2,2,2); // Q = y
-   *   P += Q; // P is equal to x+y. 
+   *   P += Q; // P is equal to x+y.
    *   P *= Q; // P is equal to xy + y^2
-   *   bgeot::power_index pi(P.dim()); 
+   *   bgeot::power_index pi(P.dim());
    *   bgeot::polynomial<double>::const_iterator ite = Q.end();
    *   bgeot::polynomial<double>::const_iterator itb = Q.begin();
    *   for ( ; itb != ite; ++itb, ++pi)
@@ -146,7 +146,7 @@ namespace bgeot
    *       have the same degree. The index of the monomial
    *       @f$ x_0^{i_0}x_1^{i_1} ... x_{n-1}^{i_{n-1}} @f$
    *       is then
-   *       @f$ \alpha_{d-1}^{n} + \alpha_{d-i_0-1}^{n-1} 
+   *       @f$ \alpha_{d-1}^{n} + \alpha_{d-i_0-1}^{n-1}
    *          + \alpha_{d-i_0-i_1-1}^{n-2} + ... + \alpha_{i_{n-1}-1}^{1}, @f$
    *       where @f$d = \sum_{l=0}^{n-1} address@hidden is the degree of the 
monomial.
    *       (by convention @f$\alpha_{-1}^{n} = address@hidden).
@@ -162,10 +162,10 @@ namespace bgeot
    *        make the operations @f$a = i_{n-1}; i_{n-1} = 0; i_{l+1} = a+1;
    *        \mbox{ if } l \ge 0 \mbox{ then } i_l = i_l - address@hidden
    *
-   *        To take the previous coefficient, let @address@hidden be the last 
index 
+   *        To take the previous coefficient, let @address@hidden be the last 
index
    *        between 0 and @address@hidden such that @f$i_l \ne address@hidden 
(if there is not, there
    *        is no previous monomial) then make the operations @f$a = i_l;
-   *        i_l = 0; i_{n-1} = a - 1; \mbox{ if } l \ge 1 \mbox{ then } 
+   *        i_l = 0; i_{n-1} = a - 1; \mbox{ if } l \ge 1 \mbox{ then }
    *        i_{l-1} = i_{l-1} + address@hidden
    *
    *  <h3>Direct product multiplication.</h3>
@@ -179,24 +179,24 @@ namespace bgeot
    */
   template<typename T> class polynomial : public std::vector<T>, public 
dal::static_stored_object {
   protected :
-    
+
     short_type n, d;
-    
+
   public :
-    
+
     typedef typename std::vector<T>::iterator iterator;
     typedef typename std::vector<T>::const_iterator const_iterator;
     typedef typename std::vector<T>::reverse_iterator reverse_iterator;
     typedef typename std::vector<T>::const_reverse_iterator 
const_reverse_iterator;
 
     /// Gives the degree of the polynomial
-    short_type degree(void) const { return d; }
+    short_type degree() const { return d; }
     /**  gives the degree of the polynomial, considering only non-zero
      * coefficients
      */
-    short_type real_degree(void) const;
+    short_type real_degree() const;
     ///     Gives the dimension (number of variables)
-    short_type dim(void) const { return n; }
+    short_type dim() const { return n; }
     /// Change the degree of the polynomial to d.
     void change_degree(short_type dd);
     /** Add to the polynomial a monomial of coefficient a and
@@ -213,10 +213,10 @@ namespace bgeot
     /// Subtract Q from P.
     polynomial operator -(const polynomial &Q) const
       { polynomial R = *this; R -= Q; return R; }
-    polynomial operator -(void) const;
+    polynomial operator -() const;
     /// Multiply P with Q. P contains the result.
     polynomial &operator *=(const polynomial &Q);
-    /// Multiply P with Q. 
+    /// Multiply P with Q.
     polynomial operator *(const polynomial &Q) const;
     /** Product of P and Q considering that variables of Q come after
      * variables of P. P contains the result
@@ -230,12 +230,12 @@ namespace bgeot
     polynomial &operator /=(const T &e);
     /// Divide P with the scalar a.
     polynomial operator /(const T &e) const
-      { polynomial res = *this; res /= e; return res; }   
+      { polynomial res = *this; res /= e; return res; }
     /// operator ==.
-    bool operator ==(const polynomial &Q) const; 
+    bool operator ==(const polynomial &Q) const;
     /// operator !=.
     bool operator !=(const polynomial &Q) const
-    { return !(operator ==(*this,Q)); }   
+    { return !(operator ==(*this,Q)); }
     /// Derivative of P with respect to the variable k. P contains the result.
     void derivative(short_type k);
     /// Makes P = 1.
@@ -244,14 +244,14 @@ namespace bgeot
     bool is_zero()
     { return(this->real_degree()==0) && (this->size()==0 || (*this)[0]==T(0)); 
}
     template <typename ITER> T horner(power_index &mi, short_type k,
-                                  short_type de, const ITER &it) const;
+                                      short_type de, const ITER &it) const;
     /** Evaluate the polynomial. "it" is an iterator pointing to the list
      * of variables. A Horner scheme is used.
      */
     template <typename ITER> T eval(const ITER &it) const;
 
     /// Constructor.
-    polynomial(void) : std::vector<T>(1)
+    polynomial() : std::vector<T>(1)
       { n = 0; d = 0; (*this)[0] = 0.0; }
     /// Constructor.
     polynomial(short_type dim_, short_type degree_);
@@ -265,7 +265,7 @@ namespace bgeot
   { n = nn; d = dd; std::fill(this->begin(), this->end(), T(0)); }
 
   template<typename T> polynomial<T>::polynomial(short_type nn,
-                                                short_type dd, short_type k)
+                                                 short_type dd, short_type k)
     : std::vector<T>(alpha(nn,dd)) {
     n = nn; d = std::max(short_type(1), dd);
     std::fill(this->begin(), this->end(), T(0));
@@ -277,7 +277,7 @@ namespace bgeot
   { polynomial res = *this; res *= Q; return res; }
 
   template<typename T>
-  bool polynomial<T>::operator ==(const polynomial &Q) const { 
+  bool polynomial<T>::operator ==(const polynomial &Q) const {
     if (dim() != Q.dim()) return false;
     const_iterator it1 = this->begin(), ite1 = this->end();
     const_iterator it2 = Q.begin(), ite2 = Q.end();
@@ -287,12 +287,12 @@ namespace bgeot
     for ( ; it2 != ite2; ++it2) if (*it2 != T(0)) return false;
     return true;
   }
-  
+
   template<typename T>
   polynomial<T> polynomial<T>::operator *(const T &e) const
   { polynomial res = *this; res *= e; return res; }
-  
-  template<typename T> short_type polynomial<T>::real_degree(void) const {
+
+  template<typename T> short_type polynomial<T>::real_degree() const {
     const_reverse_iterator it = this->rbegin(), ite = this->rend();
     size_type l = this->size();
     for ( ; it != ite; ++it, --l) { if (*it != T(0)) break; }
@@ -300,13 +300,13 @@ namespace bgeot
     while (dd > 0 && alpha(n, short_type(dd-1)) > l) --dd;
     return dd;
   }
-  
+
   template<typename T> void polynomial<T>::change_degree(short_type dd) {
     this->resize(alpha(n,dd));
     if (dd > d) std::fill(this->begin() + alpha(n,d), this->end(), T(0));
     d = dd;
   }
-  
+
   template<typename T>
   void polynomial<T>::add_monomial(const T &coeff, const power_index &power) {
     size_type i = power.global_index();
@@ -314,92 +314,92 @@ namespace bgeot
     if (i >= this->size()) { change_degree(power.degree()); }
     ((*this)[i]) += coeff;
   }
-  
-  template<typename T> 
+
+  template<typename T>
   polynomial<T> &polynomial<T>::operator +=(const polynomial &Q) {
     GMM_ASSERT2(Q.dim() == dim(), "dimensions mismatch");
-    
+
     if (Q.degree() > degree()) change_degree(Q.degree());
     iterator it = this->begin();
     const_iterator itq = Q.begin(), ite = Q.end();
     for ( ; itq != ite; ++itq, ++it) *it += *itq;
     return *this;
   }
-  
-  template<typename T> 
+
+  template<typename T>
   polynomial<T> &polynomial<T>::operator -=(const polynomial &Q) {
     GMM_ASSERT2(Q.dim() == dim() && dim() != 0, "dimensions mismatch");
-    
+
     if (Q.degree() > degree()) change_degree(Q.degree());
     iterator it = this->begin();
     const_iterator itq = Q.begin(), ite = Q.end();
     for ( ; itq != ite; ++itq, ++it) *it -= *itq;
     return *this;
   }
-  
-  template<typename T> 
-  polynomial<T> polynomial<T>::operator -(void) const {
+
+  template<typename T>
+  polynomial<T> polynomial<T>::operator -() const {
     polynomial<T> Q = *this;
     iterator itq = Q.begin(), ite = Q.end();
     for ( ; itq != ite; ++itq) *itq = -(*itq);
     return Q;
   }
-  
+
   template<typename T>
   polynomial<T> &polynomial<T>::operator *=(const polynomial &Q) {
     GMM_ASSERT2(Q.dim() == dim(), "dimensions mismatch");
-    
+
     polynomial aux = *this;
     change_degree(0); (*this)[0] = T(0);
-    
+
     power_index miq(Q.dim()), mia(dim()), mitot(dim());
     if (dim() > 0) miq[dim()-1] = Q.degree();
     const_reverse_iterator itq = Q.rbegin(), ite = Q.rend();
     for ( ; itq != ite; ++itq, --miq) {
       if (*itq != T(0)) {
-       reverse_iterator ita = aux.rbegin(), itae = aux.rend();
-       std::fill(mia.begin(), mia.end(), 0);
-       if (dim() > 0) mia[dim()-1] = aux.degree();
-       for ( ; ita != itae; ++ita, --mia)
-         if (*ita != T(0)) {
-           power_index::iterator mita = mia.begin(), mitq = miq.begin();
-           power_index::iterator mit = mitot.begin(), mite = mia.end();
-           for ( ; mita != mite; ++mita, ++mitq, ++mit)
-             *mit = short_type((*mita) + (*mitq)); /* on pourrait calculer
-                                          directement l'index global. */
-           //       cerr << "*= : " << *this << ", itq*ita="
-           //            << (*itq) * (*ita) << endl;
-           //       cerr << " itq = " << *itq << endl;
-           //       cerr << " ita = " << *ita << endl;
-           add_monomial((*itq) * (*ita), mitot);
-          
-         }
+        reverse_iterator ita = aux.rbegin(), itae = aux.rend();
+        std::fill(mia.begin(), mia.end(), 0);
+        if (dim() > 0) mia[dim()-1] = aux.degree();
+        for ( ; ita != itae; ++ita, --mia)
+          if (*ita != T(0)) {
+            power_index::iterator mita = mia.begin(), mitq = miq.begin();
+            power_index::iterator mit = mitot.begin(), mite = mia.end();
+            for ( ; mita != mite; ++mita, ++mitq, ++mit)
+              *mit = short_type((*mita) + (*mitq)); /* on pourrait calculer
+                                           directement l'index global. */
+            //             cerr << "*= : " << *this << ", itq*ita="
+            //            << (*itq) * (*ita) << endl;
+            //             cerr << " itq = " << *itq << endl;
+            //             cerr << " ita = " << *ita << endl;
+            add_monomial((*itq) * (*ita), mitot);
+
+          }
       }
     }
     return *this;
   }
 
   template<typename T>
-    void polynomial<T>::direct_product(const polynomial &Q) { 
+    void polynomial<T>::direct_product(const polynomial &Q) {
     polynomial aux = *this;
 
     change_degree(0); n = short_type(n+Q.dim()); (*this)[0] = T(0);
-    
+
     power_index miq(Q.dim()), mia(aux.dim()), mitot(dim());
     if (Q.dim() > 0) miq[Q.dim()-1] = Q.degree();
     const_reverse_iterator itq = Q.rbegin(), ite = Q.rend();
     for ( ; itq != ite; ++itq, --miq)
       if (*itq != T(0)) {
-       reverse_iterator ita = aux.rbegin(), itae = aux.rend();
-       std::fill(mia.begin(), mia.end(), 0); 
-       if (aux.dim() > 0) mia[aux.dim()-1] = aux.degree();
-       for ( ; ita != itae; ++ita, --mia)
-         if (*ita != T(0)) {
-           std::copy(mia.begin(), mia.end(), mitot.begin());
-           std::copy(miq.begin(), miq.end(), mitot.begin() + aux.dim());
-           add_monomial((*itq) * (*ita), mitot); /* on pourrait calculer
-                                          directement l'index global. */
-         }
+        reverse_iterator ita = aux.rbegin(), itae = aux.rend();
+        std::fill(mia.begin(), mia.end(), 0);
+        if (aux.dim() > 0) mia[aux.dim()-1] = aux.degree();
+        for ( ; ita != itae; ++ita, --mia)
+          if (*ita != T(0)) {
+            std::copy(mia.begin(), mia.end(), mitot.begin());
+            std::copy(miq.begin(), miq.end(), mitot.begin() + aux.dim());
+            add_monomial((*itq) * (*ita), mitot); /* on pourrait calculer
+                                           directement l'index global. */
+          }
       }
   }
 
@@ -420,9 +420,9 @@ namespace bgeot
   template<typename T>
   inline void polynomial<T>::derivative(short_type k) {
     GMM_ASSERT2(k < n, "index out of range");
-    
+
      iterator it = this->begin(), ite = this->end();
-     power_index mi(dim()); 
+     power_index mi(dim());
      for ( ; it != ite; ++it) {
        if ((*it) != T(0) && mi[k] > 0)
          { mi[k]--; (*this)[mi.global_index()] = (*it) * T(mi[k] + 1); 
mi[k]++; }
@@ -433,16 +433,16 @@ namespace bgeot
   }
 
    template<typename T> template<typename ITER>
-  inline T polynomial<T>::horner(power_index &mi, short_type k, 
-                                short_type de, const ITER &it) const {
+  inline T polynomial<T>::horner(power_index &mi, short_type k,
+                                 short_type de, const ITER &it) const {
     if (k == 0)
       return (*this)[mi.global_index()];
     else {
       T v = (*(it+k-1)), res = T(0);
       for (mi[k-1] = short_type(degree() - de); mi[k-1] != short_type(-1);
-          (mi[k-1])--)
-       res = horner(mi, short_type(k-1), short_type(de + mi[k-1]), it)
-         + v * res;
+           (mi[k-1])--)
+        res = horner(mi, short_type(k-1), short_type(de + mi[k-1]), it)
+          + v * res;
       mi[k-1] = 0;
       return res;
     }
@@ -460,73 +460,73 @@ namespace bgeot
       for (size_type i=0; i < dim(); ++i) s += it[i]*P[i+1];
       return s;
     }
- 
+
     switch (dim()) {
       case 1: {
-       T x = it[0];
-       if (deg == 2)     return P[0] + x*(P[1] + x*(P[2]));
-       if (deg == 3)     return P[0] + x*(P[1] + x*(P[2] + x*(P[3])));
-       if (deg == 4)     return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + 
x*(P[4]))));
-       if (deg == 5)     return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4] + 
x*(P[5])))));
-       if (deg == 6)     return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4] + 
x*(P[5] + x*(P[6]))))));
+        T x = it[0];
+        if (deg == 2)     return P[0] + x*(P[1] + x*(P[2]));
+        if (deg == 3)     return P[0] + x*(P[1] + x*(P[2] + x*(P[3])));
+        if (deg == 4)     return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + 
x*(P[4]))));
+        if (deg == 5)     return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4] 
+ x*(P[5])))));
+        if (deg == 6)     return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4] 
+ x*(P[5] + x*(P[6]))))));
       } break;
       case 2: {
-       T x = it[0];
-       T y = it[1];
-       if (deg == 2)     return P[0] + x*(P[1] + x*(P[3])) + y*(P[2] + 
x*(P[4]) + y*(P[5]));
-       if (deg == 3)     return P[0] + x*(P[1] + x*(P[3] + x*(P[6]))) + 
y*(P[2] + x*(P[4] + x*(P[7])) + y*(P[5] + x*(P[8]) + y*(P[9])));
-       if (deg == 4)     return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + 
x*(P[10])))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11]))) + y*(P[5] + x*(P[8] + 
x*(P[12])) + y*(P[9] + x*(P[13]) + y*(P[14]))));
-       if (deg == 5)     return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10] 
+ x*(P[15]))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16])))) + 
y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17]))) + y*(P[9] + x*(P[13] + x*(P[18])) + 
y*(P[14] + x*(P[19]) + y*(P[20])))));
-       if (deg == 6)     return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10] 
+ x*(P[15] + x*(P[21])))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16] 
+ x*(P[22]))))) + y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17] + x*(P[23])))) + 
y*(P[9] + x*(P[13] + x*(P[18] + x*(P[24]))) + y*(P[14] + x*(P[19] + x*(P[25])) 
+ y*(P[20] + x*(P[26]) + y*(P[27]))))));
+        T x = it[0];
+        T y = it[1];
+        if (deg == 2)     return P[0] + x*(P[1] + x*(P[3])) + y*(P[2] + 
x*(P[4]) + y*(P[5]));
+        if (deg == 3)     return P[0] + x*(P[1] + x*(P[3] + x*(P[6]))) + 
y*(P[2] + x*(P[4] + x*(P[7])) + y*(P[5] + x*(P[8]) + y*(P[9])));
+        if (deg == 4)     return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + 
x*(P[10])))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11]))) + y*(P[5] + x*(P[8] + 
x*(P[12])) + y*(P[9] + x*(P[13]) + y*(P[14]))));
+        if (deg == 5)     return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10] 
+ x*(P[15]))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16])))) + 
y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17]))) + y*(P[9] + x*(P[13] + x*(P[18])) + 
y*(P[14] + x*(P[19]) + y*(P[20])))));
+        if (deg == 6)     return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10] 
+ x*(P[15] + x*(P[21])))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16] 
+ x*(P[22]))))) + y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17] + x*(P[23])))) + 
y*(P[9] + x*(P[13] + x*(P[18] + x*(P[24]))) + y*(P[14] + x*(P[19] + x*(P[25])) 
+ y*(P[20] + x*(P[26]) + y*(P[27]))))));
       } break;
       case 3: {
-       T x = it[0];
-       T y = it[1];
-       T z = it[2];
-       if (deg == 2)     return P[0] + x*(P[1] + x*(P[4])) + y*(P[2] + 
x*(P[5]) + y*(P[7])) + z*(P[3] + x*(P[6]) + y*(P[8]) + z*(P[9]));
-       if (deg == 3)     return P[0] + x*(P[1] + x*(P[4] + x*(P[10]))) + 
y*(P[2] + x*(P[5] + x*(P[11])) + y*(P[7] + x*(P[13]) + y*(P[16]))) + z*(P[3] + 
x*(P[6] + x*(P[12])) + y*(P[8] + x*(P[14]) + y*(P[17])) + z*(P[9] + x*(P[15]) + 
y*(P[18]) + z*(P[19])));
-       if (deg == 4)     return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + 
x*(P[20])))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21]))) + y*(P[7] + x*(P[13] 
+ x*(P[23])) + y*(P[16] + x*(P[26]) + y*(P[30])))) + z*(P[3] + x*(P[6] + 
x*(P[12] + x*(P[22]))) + y*(P[8] + x*(P[14] + x*(P[24])) + y*(P[17] + x*(P[27]) 
+ y*(P[31]))) + z*(P[9] + x*(P[15] + x*(P[25])) + y*(P[18] + x*(P[28]) + 
y*(P[32])) + z*(P[19] + x*(P[29]) + y*(P[33]) + z*(P[34]))));
-       if (deg == 5)     return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + x*(P[20] 
+ x*(P[35]))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] + x*(P[36])))) + 
y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38]))) + y*(P[16] + x*(P[26] + x*(P[41])) 
+ y*(P[30] + x*(P[45]) + y*(P[50]))))) + z*(P[3] + x*(P[6] + x*(P[12] + 
x*(P[22] + x*(P[37])))) + y*(P[8] + x*(P[14] + x*(P[24] + x*(P[39]))) + 
y*(P[17] + x*(P[27] + x*(P[42])) + y*(P[31] + x*(P[46]) + y*(P[51])))) + 
z*(P[9] + x*(P[15] + x*(P[25] + x*(P[40]))) + y* [...]
-       if (deg == 6)     return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + x*(P[20] 
+ x*(P[35] + x*(P[56])))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] + 
x*(P[36] + x*(P[57]))))) + y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38] + 
x*(P[59])))) + y*(P[16] + x*(P[26] + x*(P[41] + x*(P[62]))) + y*(P[30] + 
x*(P[45] + x*(P[66])) + y*(P[50] + x*(P[71]) + y*(P[77])))))) + z*(P[3] + 
x*(P[6] + x*(P[12] + x*(P[22] + x*(P[37] + x*(P[58]))))) + y*(P[8] + x*(P[14] + 
x*(P[24] + x*(P[39] + x*(P[60])))) + y*(P[17] + x* [...]
+        T x = it[0];
+        T y = it[1];
+        T z = it[2];
+        if (deg == 2)     return P[0] + x*(P[1] + x*(P[4])) + y*(P[2] + 
x*(P[5]) + y*(P[7])) + z*(P[3] + x*(P[6]) + y*(P[8]) + z*(P[9]));
+        if (deg == 3)     return P[0] + x*(P[1] + x*(P[4] + x*(P[10]))) + 
y*(P[2] + x*(P[5] + x*(P[11])) + y*(P[7] + x*(P[13]) + y*(P[16]))) + z*(P[3] + 
x*(P[6] + x*(P[12])) + y*(P[8] + x*(P[14]) + y*(P[17])) + z*(P[9] + x*(P[15]) + 
y*(P[18]) + z*(P[19])));
+        if (deg == 4)     return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + 
x*(P[20])))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21]))) + y*(P[7] + x*(P[13] 
+ x*(P[23])) + y*(P[16] + x*(P[26]) + y*(P[30])))) + z*(P[3] + x*(P[6] + 
x*(P[12] + x*(P[22]))) + y*(P[8] + x*(P[14] + x*(P[24])) + y*(P[17] + x*(P[27]) 
+ y*(P[31]))) + z*(P[9] + x*(P[15] + x*(P[25])) + y*(P[18] + x*(P[28]) + 
y*(P[32])) + z*(P[19] + x*(P[29]) + y*(P[33]) + z*(P[34]))));
+        if (deg == 5)     return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + 
x*(P[20] + x*(P[35]))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] + 
x*(P[36])))) + y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38]))) + y*(P[16] + 
x*(P[26] + x*(P[41])) + y*(P[30] + x*(P[45]) + y*(P[50]))))) + z*(P[3] + 
x*(P[6] + x*(P[12] + x*(P[22] + x*(P[37])))) + y*(P[8] + x*(P[14] + x*(P[24] + 
x*(P[39]))) + y*(P[17] + x*(P[27] + x*(P[42])) + y*(P[31] + x*(P[46]) + 
y*(P[51])))) + z*(P[9] + x*(P[15] + x*(P[25] + x*(P[40]) [...]
+        if (deg == 6)     return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + 
x*(P[20] + x*(P[35] + x*(P[56])))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] 
+ x*(P[36] + x*(P[57]))))) + y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38] + 
x*(P[59])))) + y*(P[16] + x*(P[26] + x*(P[41] + x*(P[62]))) + y*(P[30] + 
x*(P[45] + x*(P[66])) + y*(P[50] + x*(P[71]) + y*(P[77])))))) + z*(P[3] + 
x*(P[6] + x*(P[12] + x*(P[22] + x*(P[37] + x*(P[58]))))) + y*(P[8] + x*(P[14] + 
x*(P[24] + x*(P[39] + x*(P[60])))) + y*(P[1 [...]
       } break;
     }
 
     /*
     switch (deg) {
       case 0: return (*this)[0];
-      case 1: { 
-       T s = (*this)[0];
-       for (size_type i=0; i < dim(); ++i) s += it[i]*(*this)[i+1];
-       return s; 
+      case 1: {
+        T s = (*this)[0];
+        for (size_type i=0; i < dim(); ++i) s += it[i]*(*this)[i+1];
+        return s;
       }
       case 2:
-      case 3: { 
-       if (dim() == 1) {
-         const T &x = it[0]; 
-         if      (deg == 2) return p[0] + x*(p[1] + x*p[2]);
-         else if (deg == 3) return p[0] + x*(p[1] + x*(p[2]+x*p[3]));
-       } else if (dim() == 2) {
-         const T &x = it[0]; 
-         const T &y = it[1]; 
-         if      (deg == 2) 
-           return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y;
-         else if (deg == 3)
-           return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y + 
-             p[6]*x*x*x + p[7]*x*x*y + p[8]*x*y*y + p[9]*y*y*y;
-       } else if (dim() == 3) {
-         const T &x = it[0]; 
-         const T &y = it[1]; 
-         const T &z = it[2]; 
-         if (deg == 2)
-           return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y + 
p[6]*x*z +
-             p[7]*y*y + p[8]*y*z + p[9]*z*z;
-         else if (deg == 3)
-           return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y + 
p[6]*x*z +
-             p[7]*y*y + p[8]*y*z + p[9]*z*z + 
-             p[10]*x*x*x + p[11]*x*x*y + p[12]*x*x*z + p[13]*x*y*y + 
p[14]*x*y*z + p[15]*x*z*z +
-             p[16]*y*y*y + p[17]*y*y*z + p[18]*y*z*z + 
-             p[19]*z*z*z;
-       }
+      case 3: {
+        if (dim() == 1) {
+          const T &x = it[0];
+          if      (deg == 2) return p[0] + x*(p[1] + x*p[2]);
+          else if (deg == 3) return p[0] + x*(p[1] + x*(p[2]+x*p[3]));
+        } else if (dim() == 2) {
+          const T &x = it[0];
+          const T &y = it[1];
+          if      (deg == 2)
+            return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y;
+          else if (deg == 3)
+            return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y +
+              p[6]*x*x*x + p[7]*x*x*y + p[8]*x*y*y + p[9]*y*y*y;
+        } else if (dim() == 3) {
+          const T &x = it[0];
+          const T &y = it[1];
+          const T &z = it[2];
+          if (deg == 2)
+            return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y + 
p[6]*x*z +
+              p[7]*y*y + p[8]*y*z + p[9]*z*z;
+          else if (deg == 3)
+            return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y + 
p[6]*x*z +
+              p[7]*y*y + p[8]*y*z + p[9]*z*z +
+              p[10]*x*x*x + p[11]*x*x*y + p[12]*x*x*z + p[13]*x*y*y + 
p[14]*x*y*z + p[15]*x*z*z +
+              p[16]*y*y*y + p[17]*y*y*z + p[18]*y*z*z +
+              p[19]*z*z*z;
+        }
       }
       }*/
     /* for other polynomials, Horner evaluation (quite slow..) */
@@ -542,14 +542,14 @@ namespace bgeot
     power_index::const_iterator mit = mi.begin(), mite = mi.end();
     for ( ; mit != mite; ++mit, ++it)
       for (short_type l = 0; l < *mit; ++l)
-       res *= *it;
+        res *= *it;
     return res;
   }
 
 
   /// Print P to the output stream o. for instance cout << P;
   template<typename T>  std::ostream &operator <<(std::ostream &o,
-                                                    const polynomial<T>& P) { 
+                                                  const polynomial<T>& P) {
     bool first = true; size_type n = 0;
     typename polynomial<T>::const_iterator it = P.begin(), ite = P.end();
     power_index mi(P.dim());
@@ -557,19 +557,19 @@ namespace bgeot
       { o << *it; first = false; ++it; ++n; ++mi; }
     for ( ; it != ite ; ++it, ++mi ) {
       if (*it != T(0)) {
-       bool first_var = true;
-       if (!first) { if (*it < T(0)) o << " - "; else o << " + "; }
-       else if (*it < T(0)) o << "-";
-       if (gmm::abs(*it)!=T(1)) { o << gmm::abs(*it); first_var = false; }
-       for (short_type j = 0; j < P.dim(); ++j)
-         if (mi[j] != 0) {
-           if (!first_var) o << "*";
-           first_var = false;
+        bool first_var = true;
+        if (!first) { if (*it < T(0)) o << " - "; else o << " + "; }
+        else if (*it < T(0)) o << "-";
+        if (gmm::abs(*it)!=T(1)) { o << gmm::abs(*it); first_var = false; }
+        for (short_type j = 0; j < P.dim(); ++j)
+          if (mi[j] != 0) {
+            if (!first_var) o << "*";
+            first_var = false;
             if (P.dim() <= 7) o << "xyzwvut"[j];
-            else o << "x_" << j; 
-           if (mi[j] > 1) o << "^" << mi[j];
-         }
-       first = false; ++n;
+            else o << "x_" << j;
+            if (mi[j] > 1) o << "^" << mi[j];
+          }
+        first = false; ++n;
       }
     }
     if (n == 0) o << "0";
@@ -583,30 +583,30 @@ namespace bgeot
      @param subs_dim : which variable is substituted
      example: poly_subs(x+y*x^2, x+1, 0) = x+1 + y*(x+1)^2
   */
-  template<typename T>    
+  template<typename T>
   polynomial<T> poly_substitute_var(const polynomial<T>& P,
-                                   const polynomial<T>& S,
-                                   size_type subs_dim) {
+                                    const polynomial<T>& S,
+                                    size_type subs_dim) {
     GMM_ASSERT2(S.dim() == 1 && subs_dim < P.dim(),
-               "wrong arguments for polynomial substitution");
+                "wrong arguments for polynomial substitution");
     polynomial<T> res(P.dim(),0);
     bgeot::power_index pi(P.dim());
     std::vector< polynomial<T> > Spow(1);
     Spow[0] = polynomial<T>(1, 0); Spow[0].one(); // Spow stores powers of S
     for (size_type k=0; k < P.size(); ++k, ++pi) {
       if (P[k] == T(0)) continue;
-      while (pi[subs_dim] >= Spow.size()) 
-       Spow.push_back(S*Spow.back());
+      while (pi[subs_dim] >= Spow.size())
+        Spow.push_back(S*Spow.back());
       const polynomial<T>& p = Spow[pi[subs_dim]];
       bgeot::power_index pi2(pi);
       for (short_type i=0; i < p.size(); ++i) {
-       pi2[subs_dim] = i;
-       res.add_monomial(p[i]*P[k],pi2);
+        pi2[subs_dim] = i;
+        res.add_monomial(p[i]*P[k],pi2);
       }
     }
     return res;
   }
-  
+
   template<typename U, typename T>
   polynomial<T> operator *(T c, const polynomial<T> &p)
   { polynomial<T> q = p; q *= c; return q; }
@@ -614,7 +614,7 @@ namespace bgeot
   typedef polynomial<opt_long_scalar_type> base_poly;
 
   /* usual constant polynomials  */
-  
+
   inline base_poly null_poly(short_type n) { return base_poly(n, 0); }
   inline base_poly one_poly(short_type n)
   { base_poly res=base_poly(n, 0); res.one(); return res;  }
@@ -634,16 +634,16 @@ namespace bgeot
 
   template<typename T> class rational_fraction : public std::vector<T> {
   protected :
-    
+
     polynomial<T> numerator_, denominator_;
-    
+
   public :
 
     const polynomial<T> &numerator() const { return numerator_; }
     const polynomial<T> &denominator() const { return denominator_; }
 
-    short_type dim(void) const { return numerator_.dim(); }
-    
+    short_type dim() const { return numerator_.dim(); }
+
     /// Add Q to P. P contains the result.
     rational_fraction &operator +=(const rational_fraction &Q) {
       numerator_ = numerator_*Q.denominator() + Q.numerator()*denominator_;
@@ -668,7 +668,7 @@ namespace bgeot
     /// Subtract Q from P.
     rational_fraction operator -(const polynomial<T> &Q) const
     { rational_fraction R(numerator_-Q*denominator_, denominator_); return R; }
-    rational_fraction operator -(void) const
+    rational_fraction operator -() const
     { rational_fraction R(-numerator_, denominator_); return R; }
     /// Multiply P with Q. P contains the result.
     rational_fraction &operator *=(const rational_fraction &Q)
@@ -676,10 +676,10 @@ namespace bgeot
     /// Divide P by Q. P contains the result.
     rational_fraction &operator /=(const rational_fraction &Q)
     { numerator_*=Q.denominator(); denominator_*=Q.numerator(); return *this; }
-    /// Multiply P with Q. 
+    /// Multiply P with Q.
     rational_fraction operator *(const rational_fraction &Q) const
     { rational_fraction R = *this; R *= Q; return R; }
-    /// Divide P by Q. 
+    /// Divide P by Q.
     rational_fraction operator /(const rational_fraction &Q) const
     { rational_fraction R = *this; R /= Q; return R; }
     /// Multiply P with the scalar a. P contains the result.
@@ -693,22 +693,23 @@ namespace bgeot
     { denominator_ *= e; return *this; }
     /// Divide P with the scalar a.
     rational_fraction operator /(const T &e) const
-    { rational_fraction res = *this; res /= e; return res; }   
+    { rational_fraction res = *this; res /= e; return res; }
     /// operator ==.
     bool operator ==(const rational_fraction &Q) const
     { return  (numerator_==Q.numerator() && denominator_==Q.denominator()); }
     /// operator !=.
     bool operator !=(const rational_fraction  &Q) const
-    { return !(operator ==(*this,Q)); }   
+    { return !(operator ==(*this,Q)); }
     /// Derivative of P with respect to the variable k. P contains the result.
     void derivative(short_type k) {
       polynomial<T> der_num = numerator_;   der_num.derivative(k);
       polynomial<T> der_den = denominator_; der_den.derivative(k);
       if (der_den.is_zero()) {
-       if (der_num.is_zero()) this->clear(); else numerator_ = der_num;
+        if (der_num.is_zero()) this->clear();
+        else numerator_ = der_num;
       } else {
-       numerator_ = der_num * denominator_ - der_den * numerator_;
-       denominator_ =  denominator_ * denominator_;
+        numerator_ = der_num * denominator_ - der_den * numerator_;
+        denominator_ =  denominator_ * denominator_;
       }
     }
     /// Makes P = 1.
@@ -718,12 +719,13 @@ namespace bgeot
       typedef typename gmm::number_traits<T>::magnitude_type R;
       T a = numerator_.eval(it), b = denominator_.eval(it);
       if (b == T(0)) { // The better should be to evaluate the derivatives ...
-       std::vector<T> p(it, it+dim()), q(dim(), T(1));
-       R no = gmm::vect_norm2(p);
-       if (no == R(0)) no = R(1E-35); else no*=gmm::default_tol(R())*R(100000);
-       gmm::add(gmm::scaled(q, T(no)), p);
-       a = numerator_.eval(p.begin());
-       b = denominator_.eval(p.begin());
+        std::vector<T> p(it, it+dim()), q(dim(), T(1));
+        R no = gmm::vect_norm2(p);
+        if (no == R(0)) no = R(1E-35);
+        else no*=gmm::default_tol(R())*R(100000);
+        gmm::add(gmm::scaled(q, T(no)), p);
+        a = numerator_.eval(p.begin());
+        b = denominator_.eval(p.begin());
       }
       if (a != T(0)) a /= b;
       return a;
@@ -746,18 +748,18 @@ namespace bgeot
   /// Add Q to P.
   template<typename T>
   rational_fraction<T> operator +(const polynomial<T> &P,
-                                 const rational_fraction<T> &Q) {
+                                  const rational_fraction<T> &Q) {
     rational_fraction<T> R(P*Q.denominator()+Q.numerator(), Q.denominator());
     return R;
   }
   /// Subtract Q from P.
   template<typename T>
   rational_fraction<T> operator -(const polynomial<T> &P,
-                                 const rational_fraction<T> &Q) {
+                                  const rational_fraction<T> &Q) {
     rational_fraction<T> R(P*Q.denominator()-Q.numerator(), Q.denominator());
     return R;
   }
-  
+
   template<typename T>  std::ostream &operator <<
   (std::ostream &o, const rational_fraction<T>& P) {
     o << "[" << P.numerator() << "]/[" << P.denominator() << "]";
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index b039b2d..ef21e68 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -577,7 +577,8 @@ namespace getfem
         /* could be a better test for discontinuity .. */
         if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
       }
-      pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 1) 
: classical_fem(pgt, 1);
+      pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 1)
+                                         : classical_fem(pgt, 1);
       pmf->set_finite_element(cv, classical_pf1);
     }
     pmf_dof_used.add(0, pmf->nb_basic_dof());
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index 29e1665..772bc9b 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -895,8 +895,8 @@ namespace getfem {
     }
   }
 
-  static pfem gen_hierarchical_fem(fem_param_list &params,
-        std::vector<dal::pstatic_stored_object> &dependencies) {
+  static pfem gen_hierarchical_fem
+  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
                 << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
@@ -910,8 +910,8 @@ namespace getfem {
                 "Bad parameters");
     pfem p  = std::make_shared<thierach_femi_comp>(ppolycompfem(pf1.get()),
                                                    ppolycompfem(pf2.get()));
-    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;
   }
 
@@ -920,7 +920,7 @@ namespace getfem {
   /* ******************************************************************** */
 
   static pfem PK_hierarch_fem(fem_param_list &params,
-        std::vector<dal::pstatic_stored_object> &) {
+                              std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
                 << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
@@ -942,7 +942,7 @@ namespace getfem {
   }
 
   static pfem QK_hierarch_fem(fem_param_list &params,
-        std::vector<dal::pstatic_stored_object> &) {
+                              std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
                 << params.size() << " should be 2.");
     GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
@@ -1014,12 +1014,14 @@ namespace getfem {
            << k << alpha << ")," << fempk << "(1," << k << alpha << "))";
     return fem_descriptor(name.str());
   }
+
   static pfem QK_fem(fem_param_list &params,
-        std::vector<dal::pstatic_stored_object> &) {
+                     std::vector<dal::pstatic_stored_object> &) {
     return QK_fem_(params, false);
   }
+
   static pfem QK_discontinuous_fem(fem_param_list &params,
-        std::vector<dal::pstatic_stored_object> &) {
+                                   std::vector<dal::pstatic_stored_object> &) {
     return QK_fem_(params, true);
   }
 
@@ -1048,8 +1050,9 @@ namespace getfem {
     return fem_descriptor(name.str());
   }
 
-  static pfem PK_prism_discontinuous_fem(fem_param_list &params,
-        std::vector<dal::pstatic_stored_object> &) {
+  static pfem
+  PK_prism_discontinuous_fem(fem_param_list &params,
+                             std::vector<dal::pstatic_stored_object> &) {
     GMM_ASSERT1(params.size() == 2 || params.size() == 3,
                 "Bad number of parameters : "
                 << params.size() << " should be 2.");
@@ -1131,8 +1134,9 @@ namespace getfem {
   // |/        |/
   // 0----1----2
 
-   static pfem Q2_incomplete_fem(fem_param_list &params,
-        std::vector<dal::pstatic_stored_object> &dependencies) {
+   static pfem
+   Q2_incomplete_fem(fem_param_list &params,
+                     std::vector<dal::pstatic_stored_object> &deps) {
     GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
     dim_type n = 2;
     if (params.size() > 0) {
@@ -1223,8 +1227,8 @@ namespace getfem {
       p->add_node(lagrange_dof(3), base_small_vector(0.5, 1.0, 1.0));
       p->add_node(lagrange_dof(3), base_small_vector(1.0, 1.0, 1.0));
     }
-    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 pfem(p);
   }
@@ -1578,7 +1582,7 @@ namespace getfem {
     base_[nc+1] = base_[1]; base_[nc+1] *= scalar_type(1 << nc);
     for (int i = 2; i <= nc; ++i) base_[nc+1] *= base_[i];
     // Le raccord assure la continuite
-    // des possibilit�s de raccord avec du P2 existent mais il faudrait
+    // des possibilités de raccord avec du P2 existent mais il faudrait
     // modifier qlq chose (transformer les fct de base P1)
   }
 
diff --git a/src/getfem_import.cc b/src/getfem_import.cc
index 5e638bc..29f78cb 100644
--- a/src/getfem_import.cc
+++ b/src/getfem_import.cc
@@ -664,7 +664,7 @@ namespace getfem {
       } break;
       case TRI: {
         if (nnode == 3) pgt = bgeot::simplex_geotrans(2,1);
-        else if (nnode == 6) { // valid�
+        else if (nnode == 6) { // validé
           static size_type lorder[6] = {0,3,1,5,4,2};
           pgt = bgeot::simplex_geotrans(2,2);
           std::copy(lorder,lorder+nnode,order.begin());
@@ -680,7 +680,7 @@ namespace getfem {
       } break;
       case TETR: {
         if (nnode == 4) pgt = bgeot::simplex_geotrans(3,1);
-        else if (nnode == 10) { // valid�
+        else if (nnode == 10) { // validé
           static size_type lorder[10] = {0,4,1, 7,8, 3, 6, 5, 9, 2};
           pgt = bgeot::simplex_geotrans(3,2);
           std::copy(lorder,lorder+nnode,order.begin());
@@ -1194,7 +1194,7 @@ namespace getfem {
     using namespace std;
     gmm::stream_standard_locale sl(f);
 
-    ofstream fichier_GiD("noboite_to_GiD.gid",    ios::out | ios::trunc );  
//d�claration du flux et ouverture du fichier
+    ofstream fichier_GiD("noboite_to_GiD.gid",    ios::out | ios::trunc );  
//déclaration du flux et ouverture du fichier
 
     fichier_GiD << "MESH    dimension 3 ElemType Tetrahedra  Nnode 4"<<endl;
 
@@ -1257,21 +1257,21 @@ namespace getfem {
     }
     fichier_GiD << "end elements" <<endl<<endl;
 
-    if(fichier_GiD)  // si l'ouverture a r�ussi
+    if(fichier_GiD)  // si l'ouverture a réussi
       {
         // instructions
         fichier_GiD.close();  // on referme le fichier
       }
     else  // sinon
-      cerr << "Erreur � l'ouverture !" << endl;
+      cerr << "Erreur à l'ouverture !" << endl;
 
-    if(f)  // si l'ouverture a r�ussi
+    if(f)  // si l'ouverture a réussi
       {
         // instructions
         //f.close();  // on referme le fichier
       }
     else  // sinon
-      cerr << "Erreur � l'ouverture !" << endl;
+      cerr << "Erreur à l'ouverture !" << endl;
 
     // appeler sunroutine import_gid_mesh_file
     //import_mesh(const std::string& "noboite_to_GiD.gid", mesh& msh)
diff --git a/src/getfem_mesh_fem.cc b/src/getfem_mesh_fem.cc
index 0816453..8217a5a 100644
--- a/src/getfem_mesh_fem.cc
+++ b/src/getfem_mesh_fem.cc
@@ -28,43 +28,42 @@
 namespace getfem {
 
   void mesh_fem::update_from_context() const {
-    for (dal::bv_visitor i(fe_convex); !i.finished(); ++i) {
-      if (linked_mesh_->convex_index().is_in(i)) {
-        if (v_num_update < linked_mesh_->convex_version_number(i)) {
+    for (dal::bv_visitor cv(fe_convex); !cv.finished(); ++cv) {
+      if (linked_mesh_->convex_index().is_in(cv)) {
+        if (v_num_update < linked_mesh_->convex_version_number(cv)) {
           if (auto_add_elt_pf != 0)
             const_cast<mesh_fem *>(this)
-              ->set_finite_element(i, auto_add_elt_pf);
+              ->set_finite_element(cv, auto_add_elt_pf);
           else if (auto_add_elt_K != dim_type(-1)) {
             if (auto_add_elt_disc)
               const_cast<mesh_fem *>(this)
-                ->set_classical_discontinuous_finite_element(i,auto_add_elt_K,
-                                                         auto_add_elt_alpha);
+                ->set_classical_discontinuous_finite_element
+                  (cv, auto_add_elt_K, auto_add_elt_alpha);
             else
               const_cast<mesh_fem *>(this)
-                ->set_classical_finite_element(i, auto_add_elt_K);
+                ->set_classical_finite_element(cv, auto_add_elt_K);
           }
           else
-            const_cast<mesh_fem *>(this)
-              ->set_finite_element(i, 0);
+            const_cast<mesh_fem *>(this)->set_finite_element(cv, 0);
         }
       }
-      else const_cast<mesh_fem *>(this)->set_finite_element(i, 0);
+      else const_cast<mesh_fem *>(this)->set_finite_element(cv, 0);
     }
-    for (dal::bv_visitor i(linked_mesh_->convex_index());
-         !i.finished(); ++i) {
-      if (!fe_convex.is_in(i)
-          && v_num_update < linked_mesh_->convex_version_number(i)) {
+    for (dal::bv_visitor cv(linked_mesh_->convex_index());
+         !cv.finished(); ++cv) {
+      if (!fe_convex.is_in(cv)
+          && v_num_update < linked_mesh_->convex_version_number(cv)) {
         if (auto_add_elt_pf != 0)
           const_cast<mesh_fem *>(this)
-            ->set_finite_element(i, auto_add_elt_pf);
+            ->set_finite_element(cv, auto_add_elt_pf);
         else if (auto_add_elt_K != dim_type(-1)) {
           if (auto_add_elt_disc)
             const_cast<mesh_fem *>(this)
-              ->set_classical_discontinuous_finite_element(i, auto_add_elt_K,
-                                                           auto_add_elt_alpha);
+              ->set_classical_discontinuous_finite_element
+                (cv, auto_add_elt_K, auto_add_elt_alpha);
           else
             const_cast<mesh_fem *>(this)
-              ->set_classical_finite_element(i, auto_add_elt_K);
+              ->set_classical_finite_element(cv, auto_add_elt_K);
         }
       }
     }



reply via email to

[Prev in Thread] Current Thread [Next in Thread]