getfem_fem.cc

Go to the documentation of this file.

/*===========================================================================
 
 Copyright (C) 1999-2024 Yves Renard
 
 This file is a part of GetFEM
 
 GetFEM  is  free software;  you  can  redistribute  it  and/or modify it
 under  the  terms  of the  GNU  Lesser General Public License as published
 by  the  Free Software Foundation;  either version 3 of the License,  or
 (at your option) any later version along with the GCC Runtime Library
 Exception either version 3.1 or (at your option) any later version.
 This program  is  distributed  in  the  hope  that it will be useful,  but
 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 License and GCC Runtime Library Exception for more details.
 You  should  have received a copy of the GNU Lesser General Public License
 along  with  this program;  if not, write to the Free Software Foundation,
 Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
 
===========================================================================*/
 
/** @file getfem_fem.cc
    @author Yves Renard <Yves.Renard@insa-lyon.fr>
    @date December 21, 1999.
    @brief implementation of various finite elements.
 */
 
#include "getfem/bgeot_torus.h"
#include "getfem/dal_singleton.h"
#include "getfem/dal_tree_sorted.h"
#include "gmm/gmm_algobase.h"
#include "getfem/dal_naming_system.h"
#include "getfem/getfem_fem.h"
#include "getfem/getfem_integration.h"
#include "getfem/getfem_omp.h"
#include "getfem/getfem_torus.h"
#include "getfem/getfem_assembling.h"
 
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) {
      GMM_ASSERT2(have_pgt() && have_G() && have_pf(), "cannot compute M");
      M_.resize(pf_->nb_dof(convex_num()), pf_->nb_base(convex_num()));
      pf_->mat_trans(M_,G(),pgt());
    }
    return M_;
  }
 
  size_type fem_interpolation_context::convex_num() const {
    GMM_ASSERT3(convex_num_ != size_type(-1), "");
    return convex_num_;
  }
 
  bool fem_interpolation_context::is_convex_num_valid() const {
    return convex_num_ != size_type(-1);
  }
 
  short_type fem_interpolation_context::face_num() const {
    GMM_ASSERT3(face_num_ != short_type(-1),
                "Face number is asked but not defined");
    return face_num_;
  }
 
  bool fem_interpolation_context::is_on_face() const {
    return (face_num_ != short_type(-1));
  }
 
  // Comment regarding fem defined on the real element:
  // Precomputed values, gradients and hessians of fem defined on the real
  // element are not dealt with within fem_interpolation_context.
  // In that case, any precomputations can be performed within the fem
  // itself, which has to store the corresponding values internally.
  // The methods real_base_value, real_grad_base_value or real_hess_base_value
  // of the fem have the possibility to check if the passed ctx has a pfp
  // and extract the corresponding internally stored results based on
  // ctx.convex_num(), ctx.pfp()->get_ppoint_tab() and ctx.ii().
  // In that case, the storage available in ctx.pfp()->c, ctx.pfp()->pc
  // and ctx.pfp()->hpc is not used.
 
 
  // Specific multiplication for fem_interpolation_context use.
  static inline void spec_mat_tmult_(const base_tensor &g, const base_matrix &B,
                                     base_tensor &t) {
    size_type P = B.nrows(), N = B.ncols();
    size_type M = t.adjust_sizes_changing_last(g, P);
    bgeot::mat_tmult(&(*(g.begin())), &(*(B.begin())), &(*(t.begin())),M,N,P);
  }
 
  void fem_interpolation_context::pfp_base_value(base_tensor& t,
                                                 const pfem_precomp &pfp__) {
    const pfem &pf__ = pfp__->get_pfem();
    GMM_ASSERT1(ii_ != size_type(-1), "Internal error");
 
    if (pf__->is_standard())
      t = pfp__->val(ii());
    else {
      if (pf__->is_on_real_element())
        pf__->real_base_value(*this, t);
      else {
        switch(pf__->vectorial_type()) {
        case virtual_fem::VECTORIAL_NOTRANSFORM_TYPE:
          t = pfp__->val(ii()); break;
        case virtual_fem::VECTORIAL_PRIMAL_TYPE:
          t.mat_transp_reduction(pfp__->val(ii()), K(), 1); break;
        case virtual_fem::VECTORIAL_DUAL_TYPE:
          t.mat_transp_reduction(pfp__->val(ii()), B(), 1); break;
        }
        if (!(pf__->is_equivalent())) {
          set_pfp(pfp__);
          { base_tensor u = t; t.mat_transp_reduction(u, M(), 0); }
        }
      }
    }
  }
 
 
  void fem_interpolation_context::base_value(base_tensor& t,
                                             bool withM) const {
    if (pfp_ && ii_ != size_type(-1) && pf_->is_standard())
      t = pfp_->val(ii());
    else {
      if (pf_->is_on_real_element())
        pf_->real_base_value(*this, t);
      else {
        if (pfp_ && ii_ != size_type(-1)) {
          switch(pf_->vectorial_type()) {
          case virtual_fem::VECTORIAL_NOTRANSFORM_TYPE:
            t = pfp_->val(ii()); break;
          case virtual_fem::VECTORIAL_PRIMAL_TYPE:
            t.mat_transp_reduction(pfp_->val(ii()), K(), 1); break;
          case virtual_fem::VECTORIAL_DUAL_TYPE:
            t.mat_transp_reduction(pfp_->val(ii()), B(), 1); break;
          }
        }
        else {
          switch(pf_->vectorial_type()) {
          case virtual_fem::VECTORIAL_NOTRANSFORM_TYPE:
            pf_->base_value(xref(), t); break;
          case virtual_fem::VECTORIAL_PRIMAL_TYPE:
            {
              base_tensor u; pf_->base_value(xref(), u);
              t.mat_transp_reduction(u,K(),1);
            } break;
          case virtual_fem::VECTORIAL_DUAL_TYPE:
            {
              base_tensor u; pf_->base_value(xref(), u);
              t.mat_transp_reduction(u,B(),1);
            } break;
          }
        }
        if (withM && !(pf_->is_equivalent()))
          { base_tensor u = t; t.mat_transp_reduction(u, M(), 0); }
      }
    }
  }
 
  void fem_interpolation_context::pfp_grad_base_value
  (base_tensor& t, const pfem_precomp &pfp__) {
    const pfem &pf__ = pfp__->get_pfem();
    GMM_ASSERT1(ii_ != size_type(-1), "Internal error");
 
    if (pf__->is_standard()) {
      // t.mat_transp_reduction(pfp__->grad(ii()), B(), 2);
      spec_mat_tmult_(pfp__->grad(ii()), B(), t);
    } else {
      if (pf__->is_on_real_element())
        pf__->real_grad_base_value(*this, t);
      else {
        switch(pf__->vectorial_type()) {
        case virtual_fem::VECTORIAL_PRIMAL_TYPE:
          {
            base_tensor u;
            // u.mat_transp_reduction(pfp__->grad(ii()), B(), 2);
            spec_mat_tmult_(pfp__->grad(ii()), B(), u);
            t.mat_transp_reduction(u, K(), 1);
          }
          break;
        case virtual_fem::VECTORIAL_DUAL_TYPE:
          {
            base_tensor u;
            // u.mat_transp_reduction(pfp__->grad(ii()), B(), 2);
            spec_mat_tmult_(pfp__->grad(ii()), B(), u);
            t.mat_transp_reduction(u, B(), 1);
          }
          break;
        default:
          // t.mat_transp_reduction(pfp__->grad(ii()), B(), 2);
          spec_mat_tmult_(pfp__->grad(ii()), B(), t);
        }
        if (!(pf__->is_equivalent())) {
          set_pfp(pfp__);
          base_tensor u = t; t.mat_transp_reduction(u, M(), 0);
        }
      }
    }
  }
 
 
  void fem_interpolation_context::grad_base_value(base_tensor& t,
                                                  bool withM) const {
    if (pfp_ && ii_ != size_type(-1) && pf_->is_standard()) {
      // t.mat_transp_reduction(pfp_->grad(ii()), B(), 2);
      spec_mat_tmult_(pfp_->grad(ii()), B(), t);
    } else {
      if (pf()->is_on_real_element())
        pf()->real_grad_base_value(*this, t);
      else {
        if (have_pfp() && ii() != size_type(-1)) {
          switch(pf()->vectorial_type()) {
          case virtual_fem::VECTORIAL_PRIMAL_TYPE:
            {
              base_tensor u;
              // u.mat_transp_reduction(pfp_->grad(ii()), B(), 2);
              spec_mat_tmult_(pfp_->grad(ii()), B(), u);
              t.mat_transp_reduction(u, K(), 1);
            }
            break;
          case virtual_fem::VECTORIAL_DUAL_TYPE:
            {
              base_tensor u;
              // u.mat_transp_reduction(pfp_->grad(ii()), B(), 2);
              spec_mat_tmult_(pfp_->grad(ii()), B(), u);
              t.mat_transp_reduction(u, B(), 1);
            }
            break;
          default:
            // t.mat_transp_reduction(pfp_->grad(ii()), B(), 2);
            spec_mat_tmult_(pfp_->grad(ii()), B(), t);
          }
 
        } else {
          base_tensor u;
          pf()->grad_base_value(xref(), u);
          if (u.size()) { /* only if the FEM can provide grad_base_value */
            // t.mat_transp_reduction(u, B(), 2);
            spec_mat_tmult_(u, B(), t);
            switch(pf()->vectorial_type()) {
            case virtual_fem::VECTORIAL_PRIMAL_TYPE:
              u = t; t.mat_transp_reduction(u, K(), 1); break;
            case virtual_fem::VECTORIAL_DUAL_TYPE:
              u = t; t.mat_transp_reduction(u, B(), 1);  break;
            default: break;
            }
          }
        }
        if (withM && !(pf()->is_equivalent()))
          { base_tensor u = t; t.mat_transp_reduction(u, M(), 0); }
      }
    }
  }
 
  void fem_interpolation_context::hess_base_value(base_tensor& t,
                                                  bool withM) const {
    if (pf()->is_on_real_element())
      pf()->real_hess_base_value(*this, t);
    else {
      base_tensor tt;
      if (have_pfp() && ii() != size_type(-1))
        tt = pfp()->hess(ii());
      else
        pf()->hess_base_value(xref(), tt);
 
      switch(pf()->vectorial_type()) {
      case virtual_fem::VECTORIAL_PRIMAL_TYPE:
        { base_tensor u = tt; tt.mat_transp_reduction(u, K(), 1); } break;
      case virtual_fem::VECTORIAL_DUAL_TYPE:
        { base_tensor u = tt; tt.mat_transp_reduction(u, B(), 1); } break;
      default: break;
      }
 
      if (tt.size()) { /* only if the FEM can provide hess_base_value */
        tt.adjust_sizes(tt.sizes()[0], tt.sizes()[1], gmm::sqr(tt.sizes()[2]));
        t.mat_transp_reduction(tt, B3(), 2);
        if (!pgt()->is_linear()) {
          if (have_pfp()) {
            tt.mat_transp_reduction(pfp()->grad(ii()), B32(), 2);
          } else {
            base_tensor u;
            pf()->grad_base_value(xref(), u);
            tt.mat_transp_reduction(u, B32(), 2);
          }
          t -= tt;
        }
        if (!(pf()->is_equivalent()) && withM)
          { tt = t; t.mat_transp_reduction(tt, M(), 0); }
      }
    }
  }
 
  void fem_interpolation_context::set_pfp(pfem_precomp newpfp) {
    if (pfp_ != newpfp) {
      pfp_ = newpfp;
      if (pfp_) { pf_ = pfp()->get_pfem(); }
      else pf_ = 0;
      M_.resize(0,0);
    }
  }
 
  void fem_interpolation_context::set_pf(pfem newpf) {
    if (pf_ != newpf || have_pfp()) {
      set_pfp(0);
      pf_ = newpf;
    }
  }
 
  void virtual_fem::real_base_value(const fem_interpolation_context &c,
                                    base_tensor &t, bool withM) const
  { c.base_value(t, withM); }
 
  void virtual_fem::real_grad_base_value(const fem_interpolation_context &c,
                                    base_tensor &t, bool withM) const
  { c.grad_base_value(t, withM);}
 
  void virtual_fem::real_hess_base_value(const fem_interpolation_context &c,
                                         base_tensor &t, bool withM) const
  { c.hess_base_value(t, withM); }
 
  /* ******************************************************************** */
  /*        Class for description of an interpolation dof.                */
  /* ******************************************************************** */
 
  enum ddl_type { LAGRANGE, NORMAL_DERIVATIVE, DERIVATIVE, MEAN_VALUE,
                  BUBBLE1, LAGRANGE_NONCONFORMING, GLOBAL_DOF,
                  SECOND_DERIVATIVE, NORMAL_COMPONENT, EDGE_COMPONENT,
                  IPK_CENTER};
 
  struct ddl_elem {
    ddl_type t;
    gmm::int16_type hier_degree;
    short_type hier_raff;
    size_type spec;
    bool operator < (const ddl_elem &l) const {
      if (t < l.t) return true;
      if (t > l.t) return false;
      if (hier_degree < l.hier_degree) return true;
      if (hier_degree > l.hier_degree) return false;
      if (hier_raff < l.hier_raff) return true;
      if (hier_raff > l.hier_raff) return false;
      if (spec < l.spec) return true;
      return false;
    }
    ddl_elem(ddl_type s = LAGRANGE, gmm::int16_type k = -1, short_type l = 0,
             size_type spec_ = 0)
      : t(s), hier_degree(k), hier_raff(l), spec(spec_) {}
  };
 
  struct dof_description {
    std::vector<ddl_elem> ddl_desc;
    bool linkable;
    dim_type coord_index;
    size_type xfem_index;
    bool all_faces;
 
    dof_description()
    { linkable = true; all_faces = false; coord_index = 0; xfem_index = 0; }
  };
 
  struct dof_description_comp__ {
    int operator()(const dof_description &m, const dof_description &n) const;
  };
 
  // ATTENTION : en cas de modif, changer aussi dof_description_compare,
  //             product_dof, et dof_hierarchical_compatibility.
  int dof_description_comp__::operator()(const dof_description &m,
                                         const dof_description &n) const {
    int nn = gmm::lexicographical_less<std::vector<ddl_elem> >()
      (m.ddl_desc, n.ddl_desc);
    if (nn < 0) return -1;
    if (nn > 0) return 1;
    nn = int(m.linkable) - int(n.linkable);
    if (nn < 0) return -1;
    if (nn > 0) return 1;
    nn = int(m.coord_index) - int(n.coord_index);
    if (nn < 0) return -1;
    if (nn > 0) return 1;
    nn = int(m.xfem_index) - int(n.xfem_index);
    if (nn < 0) return -1;
    if (nn > 0) return 1;
    nn = int(m.all_faces) - int(n.all_faces);
    if (nn < 0) return -1;
    if (nn > 0) return 1;
    return 0;
  }
 
  typedef dal::dynamic_tree_sorted<dof_description, dof_description_comp__> dof_d_tab;
 
  pdof_description lagrange_dof(dim_type n) {
    THREAD_SAFE_STATIC dim_type n_old = dim_type(-2);
    THREAD_SAFE_STATIC pdof_description p_old = nullptr;
    if (n != n_old) {
      dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
      dof_description l;
      l.ddl_desc.resize(n);
      std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
      p_old = &(tab[tab.add_norepeat(l)]);
      n_old = n;
    }
    return p_old;
  }
 
  pdof_description lagrange_0_dof(dim_type n) {
    THREAD_SAFE_STATIC dim_type n_old = dim_type(-2);
    THREAD_SAFE_STATIC pdof_description p_old = nullptr;
    if (n != n_old) {
      dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
      dof_description l;
      l.all_faces = true;
      l.ddl_desc.resize(n);
      l.linkable = false;
      std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
      p_old = &(tab[tab.add_norepeat(l)]);
      n_old = n;
    }
    return p_old;
  }
 
  pdof_description deg_hierarchical_dof(pdof_description p, int deg) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l = *p;
    for (size_type i = 0; i < l.ddl_desc.size(); ++i)
      l.ddl_desc[i].hier_degree = gmm::int16_type(deg);
    return &(tab[tab.add_norepeat(l)]);
  }
 
  size_type reserve_xfem_index() {
    static size_type ind = 100;
    return ind += 1000;
  }
 
  pdof_description xfem_dof(pdof_description p, size_type ind) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l = *p; l.xfem_index = ind;
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description ipk_center_dof(dim_type n, size_type k_ind) {
    THREAD_SAFE_STATIC dim_type n_old = dim_type(-2);
    THREAD_SAFE_STATIC pdof_description p_old = nullptr;
    if (n != n_old) {
      dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
      dof_description l;
      l.ddl_desc.resize(n);
      std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),
                ddl_elem(IPK_CENTER, -1, 0, k_ind));
      p_old = &(tab[tab.add_norepeat(l)]);
      n_old = n;
    }
    return p_old;
  }
  
 
 
  pdof_description to_coord_dof(pdof_description p, dim_type ct) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l = *p;
    l.coord_index = ct;
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description raff_hierarchical_dof(pdof_description p, short_type deg) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l = *p;
    for (size_type i = 0; i < l.ddl_desc.size(); ++i)
      l.ddl_desc[i].hier_raff = deg;
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description lagrange_nonconforming_dof(dim_type n) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l; l.linkable = false;
    l.ddl_desc.resize(n);
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description bubble1_dof(dim_type n) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l;
    l.ddl_desc.resize(n);
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(BUBBLE1));
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description derivative_dof(dim_type n, dim_type num_der) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l;
    l.ddl_desc.resize(n);
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
    l.ddl_desc.at(num_der) = ddl_elem(DERIVATIVE);
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description second_derivative_dof(dim_type n, dim_type num_der1,
                                         dim_type num_der2) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l;
    l.ddl_desc.resize(n);
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
    l.ddl_desc[num_der1] = ddl_elem(SECOND_DERIVATIVE);
    l.ddl_desc[num_der2] = ddl_elem(SECOND_DERIVATIVE);
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description normal_derivative_dof(dim_type n) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l;
    l.ddl_desc.resize(n);
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),
              ddl_elem(NORMAL_DERIVATIVE));
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description normal_component_dof(dim_type n) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l;
    l.ddl_desc.resize(n);
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),
              ddl_elem(NORMAL_COMPONENT));
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description edge_component_dof(dim_type n) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l;
    l.ddl_desc.resize(n);
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),
              ddl_elem(EDGE_COMPONENT));
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description mean_value_dof(dim_type n) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l;
    l.ddl_desc.resize(n);
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(MEAN_VALUE));
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description global_dof(dim_type n) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    dof_description l;
    l.all_faces = true;
    l.ddl_desc.resize(n);
    l.linkable = false;
    std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),ddl_elem(GLOBAL_DOF));
    return &(tab[tab.add_norepeat(l)]);
  }
 
  pdof_description product_dof(pdof_description a, pdof_description b) {
    dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
    size_type nb1 = a->ddl_desc.size(), nb2 = b->ddl_desc.size();
    dof_description l;
    l.linkable = a->linkable && b->linkable;
    l.coord_index = std::max(a->coord_index, b->coord_index); // logique ?
    l.xfem_index = a->xfem_index;
    l.all_faces = a->all_faces || b->all_faces;
    GMM_ASSERT2(a->xfem_index == b->xfem_index, "Invalid product of dof");
    l.ddl_desc.resize(nb1+nb2);
    std::copy(a->ddl_desc.begin(), a->ddl_desc.end(), l.ddl_desc.begin());
    std::copy(b->ddl_desc.begin(), b->ddl_desc.end(), l.ddl_desc.begin()+nb1);
 
    {
      gmm::int16_type deg = -1;
      for (size_type i = 0; i < l.ddl_desc.size(); ++i)
        deg = std::max(deg, l.ddl_desc[i].hier_degree);
      for (size_type i = 0; i < l.ddl_desc.size(); ++i)
        l.ddl_desc[i].hier_degree = deg;
    }
    {
      short_type deg = 0;
      for (size_type i = 0; i < l.ddl_desc.size(); ++i)
        deg = std::max(deg, l.ddl_desc[i].hier_raff);
      for (size_type i = 0; i < l.ddl_desc.size(); ++i)
        l.ddl_desc[i].hier_raff = deg;
    }
    return &(tab[tab.add_norepeat(l)]);
  }
 
  // ATTENTION : en cas de modif, changer aussi
  //             dof_description_comp__::operator,
  //             product_dof, et dof_hierarchical_compatibility.
  //             et dof_description_compare
  int dof_weak_compatibility(pdof_description a, pdof_description b) {
    int nn;
    std::vector<ddl_elem>::const_iterator
      ita = a->ddl_desc.begin(), itae = a->ddl_desc.end(),
      itb = b->ddl_desc.begin(), itbe = b->ddl_desc.end();
    for (; ita != itae && itb != itbe; ++ita, ++itb)
    {
      if ((nn = int(ita->t) - int (itb->t)) != 0) return nn;
      if ((nn = int(ita->hier_degree) - int (itb->hier_degree)) != 0)
        return nn;
    }
    for (; ita != itae; ++ita) if (ita->t != LAGRANGE) return 1;
    for (; itb != itbe; ++itb) if (itb->t != LAGRANGE) return -1;
    return 0;
  }
 
 
  // ATTENTION : en cas de modif, changer aussi
  //             dof_description_comp__::operator,
  //             product_dof, et dof_hierarchical_compatibility.
  int dof_description_compare(pdof_description a, pdof_description b) {
    if (a == b) return 0;
    int nn;
    if ((nn = int(a->coord_index) - int(b->coord_index)) != 0) return nn;
    if ((nn = int(a->linkable) - int(b->linkable)) != 0) return nn;
    if ((nn = int(a->xfem_index) - int(b->xfem_index)) != 0) return nn;
    return dof_weak_compatibility(a,b);
  }
 
  bool dof_linkable(pdof_description a)
  { return a->linkable; }
 
  bool dof_compatibility(pdof_description a, pdof_description b)
  { return (dof_linkable(a) && dof_description_compare(a, b) == 0); }
 
  size_type dof_xfem_index(pdof_description a)
  { return a->xfem_index; }
 
  dim_type coord_index_of_dof(pdof_description a)
  { return a->coord_index; }
 
  bool dof_hierarchical_compatibility(pdof_description a, pdof_description b)
  {
    if (a->coord_index != b->coord_index) return false;
    if (a->linkable != b->linkable) return false;
    if (a->xfem_index != b->xfem_index) return false;
    std::vector<ddl_elem>::const_iterator
      ita = a->ddl_desc.begin(), itae = a->ddl_desc.end(),
      itb = b->ddl_desc.begin(), itbe = b->ddl_desc.end();
    for (; ita != itae && itb != itbe; ++ita, ++itb)
    { if (ita->t != itb->t) return false; }
    for (; ita != itae; ++ita) if (ita->t != LAGRANGE) return false;
    for (; itb != itbe; ++itb) if (itb->t != LAGRANGE) return false;
    return true;
  }
 
 
  /* ******************************************************************** */
  /*        Members methods of virtual_fem           .                    */
  /* ******************************************************************** */
 
 
  void virtual_fem::add_node(const pdof_description &d, const base_node &pt,
                             const dal::bit_vector &faces) {
    short_type nb = cv_node.nb_points();
    cv_node.points().resize(nb+1);
    cv_node.points()[nb] = pt;
    dof_types_.resize(nb+1);
    face_tab.resize(nb+1);
    dof_types_[nb] = d;
    cvs_node->add_point_adaptative(nb, short_type(-1));
    for (dal::bv_visitor f(faces); !f.finished(); ++f) {
      cvs_node->add_point_adaptative(nb, short_type(f));
      face_tab[nb].push_back(short_type(f));
    }
    pspt_valid = false;
  }
 
 
  void virtual_fem::add_node(const pdof_description &d, const base_node &pt) {
    dal::bit_vector faces;
    for (short_type f = 0; f < cvs_node->nb_faces(); ++f)
      if (d->all_faces || gmm::abs(cvr->is_in_face(f, pt)) < 1.0E-7)
        faces.add(f);
     add_node(d, pt, faces);
  }
 
  void virtual_fem::init_cvs_node() {
    cvs_node->init_for_adaptative(cvr->structure());
    cv_node = bgeot::convex<base_node>(cvs_node);
    face_tab.resize(0);
    pspt_valid = false;
  }
 
  void virtual_fem::unfreeze_cvs_node() {
    cv_node.structure() = cvs_node;
    pspt_valid = false;
  }
 
  const std::vector<short_type> &
  virtual_fem::faces_of_dof(size_type /*cv*/, size_type i) const {
    static const std::vector<short_type> no_faces;
    return (i < face_tab.size()) ? face_tab[i] : no_faces;
  }
 
  void virtual_fem::copy(const virtual_fem &f) {
    dof_types_ = f.dof_types_;
    cvs_node = bgeot::new_convex_structure();
    *cvs_node = *f.cvs_node; // deep copy
    cv_node = f.cv_node;
    cv_node.structure() = cvs_node;
    pspt = 0;
    pspt_valid = false;
    cvr = f.cvr;
    dim_ = f.dim_;
    ntarget_dim = f.ntarget_dim;
    vtype = f.vtype;
    is_equiv = f.is_equiv;
    is_lag = f.is_lag;
    is_pol = f.is_pol;
    is_polycomp = f.is_polycomp;
    real_element_defined = f.real_element_defined;
    is_standard_fem = f.is_standard_fem;
    es_degree = f.es_degree;
    hier_raff = f.hier_raff;
    debug_name_ = f.debug_name_;
    face_tab = f.face_tab;
  }
 
  /* ******************************************************************** */
  /*        PK class.                                                         */
  /* ******************************************************************** */
 
  class PK_fem_ : public fem<base_poly> {
  public :
    void calc_base_func(base_poly &p, size_type i, short_type K) const;
    PK_fem_(dim_type nc, short_type k);
    ~PK_fem_() {}
  };
 
  void PK_fem_::calc_base_func(base_poly &p, size_type i, short_type K) const {
    dim_type N = dim();
    base_poly l0(N, 0), l1(N, 0);
    bgeot::power_index w(short_type(N+1));
    l0.one(); l1.one(); p = l0;
 
    if (K != 0) {
      for (short_type nn = 0; nn < N; ++nn) l0 -= base_poly(N, 1, nn);
 
      w[0] = K;
      for (short_type nn = 1; nn <= N; ++nn) {
        w[nn]=short_type(floor(0.5+bgeot::to_scalar((cv_node.points()[i])[nn-1]*opt_long_scalar_type(K))));
        w[0]=short_type(w[0] - w[nn]);
      }
 
      for (int nn = 0; nn <= N; ++nn)
        for (int j = 0; j < w[nn]; ++j) {
          if (nn == 0)
            p *= (l0 * (opt_long_scalar_type(K) / opt_long_scalar_type(j+1)))
              - (l1 * (opt_long_scalar_type(j) / opt_long_scalar_type(j+1)));
          else
            p *= (base_poly(N,1,short_type(nn-1)) * (opt_long_scalar_type(K)/opt_long_scalar_type(j+1)))
              - (l1 * (opt_long_scalar_type(j) / opt_long_scalar_type(j+1)));
        }
    }
  }
 
  PK_fem_::PK_fem_(dim_type nc, short_type k) {
    cvr = bgeot::simplex_of_reference(nc);
    dim_ = cvr->structure()->dim();
    is_standard_fem = is_equiv = is_pol = is_lag = true;
    es_degree = k;
 
    init_cvs_node();
    bgeot::pconvex_ref cvn = bgeot::simplex_of_reference(nc, k);
    size_type R = cvn->nb_points();
    for (size_type i = 0; i < R; ++i)
      add_node(k==0 ? lagrange_0_dof(nc) : lagrange_dof(nc), cvn->points()[i]);
 
    base_.resize(R);
    for (size_type r = 0; r < R; r++) calc_base_func(base_[r], r, k);
  }
 
  static pfem PK_fem(fem_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() == 0 && params[1].type() == 0,
                "Bad type of parameters");
    int n = dim_type(::floor(params[0].num() + 0.01));
    int k = short_type(::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");
    pfem p = std::make_shared<PK_fem_>(dim_type(n), short_type(k));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
 
  /* ******************************************************************** */
  /*        Tensorial product of fem (for polynomial fem).                */
  /* ******************************************************************** */
 
  struct tproduct_femi : public fem<base_poly> {
    tproduct_femi(ppolyfem fi1, ppolyfem fi2);
  };
 
  tproduct_femi::tproduct_femi(ppolyfem fi1, ppolyfem fi2) {
    if (fi2->target_dim() != 1) std::swap(fi1, fi2);
    GMM_ASSERT1(fi2->target_dim() == 1, "dimensions mismatch");
 
    is_pol = true;
    is_equiv = fi1->is_equivalent() && fi2->is_equivalent();
    GMM_ASSERT1(is_equiv,
                "Product of non equivalent elements not available, sorry.");
    is_lag = fi1->is_lagrange() && fi2->is_lagrange();
    is_standard_fem = fi1->is_standard() && fi2->is_standard();
    es_degree = short_type(fi1->estimated_degree() + fi2->estimated_degree());
 
    bgeot::convex<base_node> cv
      = bgeot::convex_direct_product(fi1->node_convex(0), fi2->node_convex(0));
    cvr = bgeot::convex_ref_product(fi1->ref_convex(0), fi2->ref_convex(0));
    dim_ = cvr->structure()->dim();
    init_cvs_node();
 
    ntarget_dim = fi2->target_dim();
    base_.resize(cv.nb_points() * ntarget_dim);
    size_type i, j, r;
    for (j = 0, r = 0; j < fi2->nb_dof(0); ++j)
      for (i = 0; i < fi1->nb_dof(0); ++i, ++r)
        add_node(product_dof(fi1->dof_types()[i], fi2->dof_types()[j]),
                 cv.points()[r]);
 
    for (j = 0, r = 0; j < fi2->nb_base_components(0); j++)
      for (i = 0; i < fi1->nb_base_components(0); i++, ++r) {
        base_[r] = fi1->base()[i];
        base_[r].direct_product(fi2->base()[j]);
      }
  }
 
  static pfem product_fem(fem_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");
    pfem pf1 = params[0].method();
    pfem pf2 = params[1].method();
    GMM_ASSERT1(pf1->is_polynomial() && pf2->is_polynomial(),
                "Both arguments to FEM_PRODUCT must be polynomial FEM");
    pfem p = std::make_shared<tproduct_femi>(ppolyfem(pf1.get()),
                                             ppolyfem(pf2.get()));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*    Generic Hierarchical fem (for polynomial fem). To be interfaced.  */
  /* ******************************************************************** */
 
  struct thierach_femi : public fem<base_poly> {
    thierach_femi(ppolyfem fi1, ppolyfem fi2);
  };
 
  thierach_femi::thierach_femi(ppolyfem fi1, ppolyfem fi2)
    : fem<base_poly>(*fi1) {
    grad_computed_ = false;
    hess_computed_ = false;
    GMM_ASSERT1(fi2->target_dim()==fi1->target_dim(), "dimensions mismatch.");
    GMM_ASSERT1(fi2->basic_structure(0) == fi1->basic_structure(0),
                "Incompatible elements.");
    GMM_ASSERT1(fi1->is_equivalent() &&  fi2->is_equivalent(), "Sorry, "
                "no hierachical construction for non tau-equivalent fems.");
    es_degree = fi2->estimated_degree();
    is_lag = false;
    unfreeze_cvs_node();
    for (size_type i = 0; i < fi2->nb_dof(0); ++i) {
      bool found = false;
      for (size_type j = 0; j < fi1->nb_dof(0); ++j) {
        if ( gmm::vect_dist2(fi2->node_of_dof(0,i),
                             fi1->node_of_dof(0,j)) < 1e-13
             && dof_hierarchical_compatibility(fi2->dof_types()[i],
                                               fi1->dof_types()[j]))
            { found = true; break; }
      }
      if (!found) {
        add_node(deg_hierarchical_dof(fi2->dof_types()[i],
                                      fi1->estimated_degree()),
                 fi2->node_of_dof(0,i));
        base_.resize(nb_dof(0));
        base_[nb_dof(0)-1] = (fi2->base())[i];
      }
    }
  }
 
  struct thierach_femi_comp : public fem<bgeot::polynomial_composite> {
    thierach_femi_comp(ppolycompfem fi1, ppolycompfem fi2);
  };
 
  thierach_femi_comp::thierach_femi_comp(ppolycompfem fi1, ppolycompfem fi2)
    : fem<bgeot::polynomial_composite>(*fi1) {
    GMM_ASSERT1(fi2->target_dim()==fi1->target_dim(), "dimensions mismatch.");
    GMM_ASSERT1(fi2->basic_structure(0) == fi1->basic_structure(0),
                "Incompatible elements.");
    GMM_ASSERT1(fi1->is_equivalent() &&  fi2->is_equivalent(), "Sorry, "
                "no hierachical construction for non tau-equivalent fems.");
 
    es_degree = std::max(fi2->estimated_degree(), fi1->estimated_degree());
 
    is_lag = false;
    hier_raff = short_type(fi1->hierarchical_raff() + 1);
    unfreeze_cvs_node();
    for (size_type i = 0; i < fi2->nb_dof(0); ++i) {
      bool found = false;
      for (size_type j = 0; j < fi1->nb_dof(0); ++j) {
        if ( gmm::vect_dist2(fi2->node_of_dof(0,i),
                             fi1->node_of_dof(0,j)) < 1e-13
             && dof_hierarchical_compatibility(fi2->dof_types()[i],
                                               fi1->dof_types()[j]))
          { found = true; break; }
      }
      if (!found) {
        add_node(raff_hierarchical_dof(fi2->dof_types()[i], hier_raff),
                 fi2->node_of_dof(0,i));
        base_.resize(nb_dof(0));
        base_[nb_dof(0)-1] = (fi2->base())[i];
      }
    }
  }
 
  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,
                "Bad type of parameters");
    pfem pf1 = params[0].method();
    pfem pf2 = params[1].method();
    if (pf1->is_polynomial() && pf2->is_polynomial())
      return std::make_shared<thierach_femi>(ppolyfem(pf1.get()),
                                             ppolyfem(pf2.get()));
    GMM_ASSERT1(pf1->is_polynomialcomp() && pf2->is_polynomialcomp(),
                "Bad parameters");
    pfem p  = std::make_shared<thierach_femi_comp>(ppolycompfem(pf1.get()),
                                                   ppolycompfem(pf2.get()));
    deps.push_back(p->ref_convex(0));
    deps.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /* PK hierarchical fem.                                                 */
  /* ******************************************************************** */
 
  static pfem PK_hierarch_fem(fem_param_list &params,
                              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,
                "Bad type of parameters");
    int n = int(::floor(params[0].num() + 0.01));
    int k = int(::floor(params[1].num() + 0.01)), s;
    GMM_ASSERT1(n > 0 && n < 100 && k > 0 && k <= 150 &&
                double(n) == params[0].num() && double(k) == params[1].num(),
                "Bad parameters");
    std::stringstream name;
    if (k == 1)
      name << "FEM_PK(" << n << ",1)";
    else {
      for (s = 2; s <= k; ++s) if ((k % s) == 0) break;
      name << "FEM_GEN_HIERARCHICAL(FEM_PK_HIERARCHICAL(" << n << ","
           << k/s << "), FEM_PK(" << n << "," << k << "))";
    }
    return fem_descriptor(name.str());
  }
 
  static pfem QK_hierarch_fem(fem_param_list &params,
                              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,
                "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");
    std::stringstream name;
    if (n == 1)
      name << "FEM_PK_HIERARCHICAL(1," << k << ")";
    else
      name << "FEM_PRODUCT(FEM_PK_HIERARCHICAL(" << n-1 << "," << k
           << "),FEM_PK_HIERARCHICAL(1," << k << "))";
    return fem_descriptor(name.str());
  }
 
  static pfem prism_PK_hierarch_fem(fem_param_list &params,
        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,
                "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");
    std::stringstream name;
    if (n == 2)
      name << "FEM_QK_HIERARCHICAL(1," << k << ")";
    else
      name << "FEM_PRODUCT(FEM_PK_HIERARCHICAL(" << n-1 << "," << k
           << "),FEM_PK_HIERARCHICAL(1," << k << "))";
    return fem_descriptor(name.str());
  }
 
 
  /* ******************************************************************** */
  /* parallelepiped fems.                                                 */
  /* ******************************************************************** */
 
  static pfem QK_fem_(fem_param_list &params, bool discontinuous) {
    const char *fempk = discontinuous ? "FEM_PK_DISCONTINUOUS" : "FEM_PK";
    const char *femqk = discontinuous ? "FEM_QK_DISCONTINUOUS" : "FEM_QK";
    GMM_ASSERT1(params.size() == 2 || (discontinuous && params.size() == 3),
                "Bad number of parameters : "
                << params.size() << " should be 2.");
    GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0 &&
                (params.size() != 3 || params[2].type() == 0),
                "Bad type of parameters");
    int n = int(::floor(params[0].num() + 0.01));
    int k = int(::floor(params[1].num() + 0.01));
    char alpha[128]; alpha[0] = 0;
    if (discontinuous && params.size() == 3) {
      scalar_type v = params[2].num();
      GMM_ASSERT1(v >= 0 && v <= 1, "Bad value for alpha: " << v);
      snprintf(alpha, 127, ",%g", v);
    }
    GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
                double(n) == params[0].num() && double(k) == params[1].num(),
                "Bad parameters");
    std::stringstream name;
    if (n == 1)
      name << fempk << "(1," << k << alpha << ")";
    else
      name << "FEM_PRODUCT(" << femqk << "(" << n-1 << ","
           << 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> &) {
    return QK_fem_(params, false);
  }
 
  static pfem QK_discontinuous_fem(fem_param_list &params,
                                   std::vector<dal::pstatic_stored_object> &) {
    return QK_fem_(params, true);
  }
 
 
  /* ******************************************************************** */
  /* prims fems.                                                          */
  /* ******************************************************************** */
 
  static pfem prism_PK_fem(fem_param_list &params,
                           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,
                "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");
    std::stringstream name;
    if (n == 2)
      name << "FEM_QK(1," << k << ")";
    else
      name << "FEM_PRODUCT(FEM_PK(" << n-1 << "," << k << "),FEM_PK(1,"
           << k << "))";
    return fem_descriptor(name.str());
  }
 
  static pfem
  prism_PK_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 or 3.");
    GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0 &&
                (params.size() != 3 || params[2].type() == 0),
                "Bad type of parameters");
    int n = int(::floor(params[0].num() + 0.01));
    int k = int(::floor(params[1].num() + 0.01));
    char alpha[128]; alpha[0] = 0;
    if (params.size() == 3) {
      scalar_type v = params[2].num();
      GMM_ASSERT1(v >= 0 && v <= 1, "Bad value for alpha: " << v);
      snprintf(alpha, 127, ",%g", v);
    }
    GMM_ASSERT1(n > 1 && n < 100 && k >= 0 && k <= 150 &&
                double(n) == params[0].num() && double(k) == params[1].num(),
                "Bad parameters");
    std::stringstream name;
    if (n == 2)
      name << "FEM_QK_DISCONTINUOUS(1," << k << alpha << ")";
    else
      name << "FEM_PRODUCT(FEM_PK_DISCONTINUOUS(" << n-1 << "," << k << alpha
           << "),FEM_PK_DISCONTINUOUS(1,"
           << k << alpha << "))";
    return fem_descriptor(name.str());
  }
 
  /* ******************************************************************** */
  /*        P1 NON CONFORMING (dim 2)                                         */
  /* ******************************************************************** */
 
   static pfem P1_nonconforming_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    GMM_ASSERT1(params.size() == 0, "Bad number of parameters ");
    auto p = std::make_shared<fem<base_poly>>();
    p->mref_convex() = bgeot::simplex_of_reference(2);
    p->dim() = 2;
    p->is_standard() = p->is_equivalent() = true;
    p->is_polynomial() = p->is_lagrange() = true;
    p->estimated_degree() = 1;
    p->init_cvs_node();
    p->base().resize(3);
 
    p->add_node(lagrange_dof(2), base_small_vector(0.5, 0.5));
    read_poly(p->base()[0], 2, "2*x + 2*y - 1");
    p->add_node(lagrange_dof(2), base_small_vector(0.0, 0.5));
    read_poly(p->base()[1], 2, "1 - 2*x");
    p->add_node(lagrange_dof(2), base_small_vector(0.5, 0.0));
    read_poly(p->base()[2], 2, "1 - 2*y");
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
 
    return pfem(p);
  }
 
 
  /* ******************************************************************** */
  /*     Quad8/Hexa20 SERENDIPITY ELEMENT (dim 2 or 3) (incomplete Q2)    */
  /* ******************************************************************** */
 
  // local dof numeration for 2D:
  // 5--6--7
  // |     |
  // 3     4
  // |     |
  // 0--1--2
  //
  // local dof numeration for 3D:
  //
  //       17---18---19
  //      /|        /|
  //     / 10      / 11
  //   15  |      16 |
  //   /   5----6/---7
  //  /   /     /   /
  // 12---13---14  /
  // |  3      |  4
  // 8 /       9 /
  // |/        |/
  // 0----1----2
 
   static pfem
   build_Q2_incomplete_fem(fem_param_list &params,
                           std::vector<dal::pstatic_stored_object> &deps,
                           bool discontinuous) {
    GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
    dim_type n = 2;
    if (params.size() > 0) {
      GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
      n = dim_type(::floor(params[0].num() + 0.01));
       GMM_ASSERT1(n == 2 || n == 3, "Bad parameter, expected value 2 or 3");
    }
    auto p = std::make_shared<fem<base_poly>>();
    p->mref_convex() = bgeot::parallelepiped_of_reference(n);
    p->dim() = n;
    p->is_standard() = p->is_equivalent() = true;
    p->is_polynomial() = p->is_lagrange() = true;
    p->estimated_degree() = 2;
    p->init_cvs_node();
    p->base().resize(n == 2 ? 8 : 20);
    auto lag_dof = discontinuous ? lagrange_nonconforming_dof(n) 
                                 : lagrange_dof(n);
 
    if (n == 2) {
      std::stringstream s
        ( "1 - 2*x^2*y - 2*x*y^2 + 2*x^2 + 5*x*y + 2*y^2 - 3*x - 3*y;"
          "4*(x^2*y - x^2 - x*y + x);"
          "2*x*y*y - 2*x*x*y + 2*x*x - x*y - x;"
          "4*(x*y*y - x*y - y*y + y);"
          "4*(x*y - x*y*y);"
          "2*x*x*y - 2*x*y*y - x*y + 2*y*y - y;"
          "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] = read_base_poly(2, s);
 
      p->add_node(lag_dof, base_small_vector(0.0, 0.0));
      p->add_node(lag_dof, base_small_vector(0.5, 0.0));
      p->add_node(lag_dof, base_small_vector(1.0, 0.0));
      p->add_node(lag_dof, base_small_vector(0.0, 0.5));
      p->add_node(lag_dof, base_small_vector(1.0, 0.5));
      p->add_node(lag_dof, base_small_vector(0.0, 1.0));
      p->add_node(lag_dof, base_small_vector(0.5, 1.0));
      p->add_node(lag_dof, base_small_vector(1.0, 1.0));
    } else {
      std::stringstream s
        ("1 + 2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2"
           " - 2*x^2*y - 2*x^2*z - 2*x*y^2 - 2*y^2*z - 2*y*z^2 - 2*x*z^2 - 7*x*y*z"
           " + 2*x^2 + 2*y^2 + 2*z^2 + 5*y*z + 5*x*z + 5*x*y - 3*x - 3*y - 3*z;"
         "4*( - x^2*y*z + x*y*z + x^2*z - x*z + x^2*y - x*y - x^2 + x);"
         "2*x^2*y*z - 2*x*y^2*z - 2*x*y*z^2"
           " - 2*x^2*y - 2*x^2*z + 2*x*y^2 + 2*x*z^2 + 3*x*y*z + 2*x^2 - x*y - x*z - x;"
         "4*( - x*y^2*z + x*y^2 + y^2*z + x*y*z - x*y - y^2 - y*z + y);"
         "4*(x*y^2*z - x*y^2 - x*y*z + x*y);"
         " - 2*x^2*y*z + 2*x*y^2*z - 2*x*y*z^2"
           " + 2*x^2*y - 2*x*y^2 - 2*y^2*z + 2*y*z^2 + 3*x*y*z - x*y + 2*y^2 - y*z - y;"
         "4*(x^2*y*z - x^2*y - x*y*z + x*y);"
         " - 2*x^2*y*z - 2*x*y^2*z + 2*x*y*z^2 + 2*x^2*y + 2*x*y^2 + x*y*z - 3*x*y;"
         "4*( - x*y*z^2 + x*z^2 + y*z^2 + x*y*z - x*z - y*z - z^2 + z);"
         "4*(x*y*z^2 - x*y*z - x*z^2 + x*z);"
         "4*(x*y*z^2 - x*y*z - y*z^2 + y*z);"
         "4*( - x*y*z^2 + x*y*z);"
         " - 2*x^2*y*z - 2*x*y^2*z + 2*x*y*z^2"
           " + 2*x^2*z + 2*y^2*z - 2*x*z^2 - 2*y*z^2 + 3*x*y*z - x*z - y*z + 2*z^2 - z;"
         "4*(x^2*y*z - x^2*z - x*y*z + x*z);"
         " - 2*x^2*y*z + 2*x*y^2*z - 2*x*y*z^2 + 2*x^2*z + 2*x*z^2 + x*y*z - 3*x*z;"
         "4*(x*y^2*z - y^2*z - x*y*z + y*z);"
         "4*( - x*y^2*z + x*y*z);"
         "2*x^2*y*z - 2*x*y^2*z - 2*x*y*z^2 + 2*y^2*z + 2*y*z^2 + x*y*z - 3*y*z;"
         "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] = read_base_poly(3, s);
 
      p->add_node(lag_dof, base_small_vector(0.0, 0.0, 0.0));
      p->add_node(lag_dof, base_small_vector(0.5, 0.0, 0.0));
      p->add_node(lag_dof, base_small_vector(1.0, 0.0, 0.0));
      p->add_node(lag_dof, base_small_vector(0.0, 0.5, 0.0));
      p->add_node(lag_dof, base_small_vector(1.0, 0.5, 0.0));
      p->add_node(lag_dof, base_small_vector(0.0, 1.0, 0.0));
      p->add_node(lag_dof, base_small_vector(0.5, 1.0, 0.0));
      p->add_node(lag_dof, base_small_vector(1.0, 1.0, 0.0));
 
      p->add_node(lag_dof, base_small_vector(0.0, 0.0, 0.5));
      p->add_node(lag_dof, base_small_vector(1.0, 0.0, 0.5));
      p->add_node(lag_dof, base_small_vector(0.0, 1.0, 0.5));
      p->add_node(lag_dof, base_small_vector(1.0, 1.0, 0.5));
 
      p->add_node(lag_dof, base_small_vector(0.0, 0.0, 1.0));
      p->add_node(lag_dof, base_small_vector(0.5, 0.0, 1.0));
      p->add_node(lag_dof, base_small_vector(1.0, 0.0, 1.0));
      p->add_node(lag_dof, base_small_vector(0.0, 0.5, 1.0));
      p->add_node(lag_dof, base_small_vector(1.0, 0.5, 1.0));
      p->add_node(lag_dof, base_small_vector(0.0, 1.0, 1.0));
      p->add_node(lag_dof, base_small_vector(0.5, 1.0, 1.0));
      p->add_node(lag_dof, base_small_vector(1.0, 1.0, 1.0));
    }
    deps.push_back(p->ref_convex(0));
    deps.push_back(p->node_tab(0));
 
    return pfem(p);
  }
 
  static pfem Q2_incomplete_fem
  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps)
  { return build_Q2_incomplete_fem(params, deps, false); }
 
  static pfem Q2_incomplete_discontinuous_fem
  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps)
  { return build_Q2_incomplete_fem(params, deps, true); }
 
  /* ******************************************************************** */
  /*    Lagrange Pyramidal element of degree 0, 1 and 2                   */
  /* ******************************************************************** */
 
  // local dof numeration for K=1:
  //    4
  //   /|||
  //  / || |
  // 2-|--|-3
  // | |  | |
  // ||    ||
  // ||    ||
  // 0------1
  //
  // local dof numeration for K=2:
  //
  //    13
  //   /  |
  //  11---12
  //  |    |
  //  9----10
  //  /     |
  // 6---7---8
  // |       |
  // 3   4   5
  // |       |
  // 0---1---2
 
  static pfem
  build_pyramid_QK_fem(short_type k, bool disc, scalar_type alpha=0) {
    auto p = std::make_shared<fem<base_rational_fraction>>();
    p->mref_convex() = bgeot::pyramid_QK_of_reference(1);
    p->dim() = 3;
    p->is_standard() = p->is_equivalent() = true;
    p->is_polynomial() = false;
    p->is_lagrange() = true;
    p->estimated_degree() = k;
    p->init_cvs_node();
    auto lag_dof = disc ? lagrange_nonconforming_dof(3) : lagrange_dof(3);
 
    if (k == 0) {
      p->base().resize(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);
      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));
      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));
      p->add_node(lag_dof, base_small_vector( 0.0,  0.0, 1.0));
 
    } else if (k == 2) {
      p->base().resize(14);
 
      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");
 
      std::vector<base_node> points = { base_node(-1.0, -1.0, 0.0),
                                        base_node( 0.0, -1.0, 0.0),
                                        base_node( 1.0, -1.0, 0.0),
                                        base_node(-1.0,  0.0, 0.0),
                                        base_node( 0.0,  0.0, 0.0),
                                        base_node( 1.0,  0.0, 0.0),
                                        base_node(-1.0,  1.0, 0.0),
                                        base_node( 0.0,  1.0, 0.0),
                                        base_node( 1.0,  1.0, 0.0),
                                        base_node(-0.5, -0.5, 0.5),
                                        base_node( 0.5, -0.5, 0.5),
                                        base_node(-0.5,  0.5, 0.5),
                                        base_node( 0.5,  0.5, 0.5),
                                        base_node( 0.0,  0.0, 1.0) };
 
      if (disc && alpha != scalar_type(0)) {
        base_node G =
          gmm::mean_value(points.begin(), points.end());
        for (auto && pt : points)
          pt = (1-alpha)*pt + alpha*G;
        for (size_type d = 0; d < 3; ++d) {
          base_poly S(1,2);
          S[0] = -alpha * G[d] / (1-alpha);
          S[1] = 1. / (1-alpha);
          xi0 = bgeot::poly_substitute_var(xi0, S, d);
          xi1 = bgeot::poly_substitute_var(xi1, S, d);
          xi2 = bgeot::poly_substitute_var(xi2, S, d);
          xi3 = bgeot::poly_substitute_var(xi3, S, d);
          x = bgeot::poly_substitute_var(x, S, d);
          y = bgeot::poly_substitute_var(y, S, d);
          z = bgeot::poly_substitute_var(z, S, d);
          un_z = bgeot::poly_substitute_var(un_z, S, d);
        }
      }
      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.;
      p->base()[ 2] = Q*Q*xi1*xi2*(-x*y-z*un_z);
      p->base()[ 3] = -Q*Q*xi3*xi0*xi1*x*4.;
      p->base()[ 4] = Q*Q*xi0*xi1*xi2*xi3*16.;
      p->base()[ 5] = Q*Q*xi1*xi2*xi3*x*4.;
      p->base()[ 6] = Q*Q*xi3*xi0*(-x*y-z*un_z);
      p->base()[ 7] = Q*Q*xi2*xi3*xi0*y*4.;
      p->base()[ 8] = Q*Q*xi2*xi3*(x*y-z*un_z);
      p->base()[ 9] = Q*z*xi0*xi1*4.;
      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] = z*(z*2-ones);
 
      for (const auto &pt : points)
        p->add_node(lag_dof, pt);
 
    } else GMM_ASSERT1(false, "Sorry, pyramidal Lagrange fem "
                       "implemented only for degree 0, 1 or 2");
 
    return pfem(p);
  }
 
 
  static pfem pyramid_QK_fem
  (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) {
      GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
      k = dim_type(::floor(params[0].num() + 0.01));
    }
    pfem p = build_pyramid_QK_fem(k, false);
    deps.push_back(p->ref_convex(0));
    deps.push_back(p->node_tab(0));
    return p;
  }
 
  static pfem pyramid_QK_disc_fem
  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps) {
    GMM_ASSERT1(params.size() <= 2, "Bad number of parameters");
    short_type k = 2;
    if (params.size() > 0) {
      GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
      k = dim_type(::floor(params[0].num() + 0.01));
    }
    scalar_type alpha(0);
    if (params.size() > 1) {
      GMM_ASSERT1(params[1].type() == 0, "Bad type of parameters");
      alpha = params[1].num();
    }
    pfem p = build_pyramid_QK_fem(k, true, alpha);
    deps.push_back(p->ref_convex(0));
    deps.push_back(p->node_tab(0));
    return p;
  }
 
 
  /* ******************************************************************** */
  /*    Incomplete quadratic Lagrange pyramidal element.                  */
  /* ******************************************************************** */
 
  // local dof numeration:
  //
  //    12
  //   /  |
  //  10---11
  //  |    |
  //  8----9
  //  /     |
  // 5---6---7
  // |       |
  // 3       4
  // |       |
  // 0---1---2
 
  static pfem build_pyramid_Q2_incomplete_fem(bool disc) {
    auto p = std::make_shared<fem<base_rational_fraction>>();
    p->mref_convex() = bgeot::pyramid_QK_of_reference(1);
    p->dim() = 3;
    p->is_standard() = p->is_equivalent() = true;
    p->is_polynomial() = false;
    p->is_lagrange() = true;
    p->estimated_degree() = 2;
    p->init_cvs_node();
    auto lag_dof = disc ? lagrange_nonconforming_dof(3) : lagrange_dof(3);
 
    p->base().resize(13);
 
    base_poly xy = read_base_poly(3, "x*y");
    base_poly z = read_base_poly(3, "z");
 
    base_poly p00m = read_base_poly(3, "1-z");
 
    base_poly pmmm = read_base_poly(3, "1-x-y-z");
    base_poly ppmm = read_base_poly(3, "1+x-y-z");
    base_poly pmpm = read_base_poly(3, "1-x+y-z");
    base_poly pppm = read_base_poly(3, "1+x+y-z");
 
    base_poly mmm0_ = read_base_poly(3, "-1-x-y")*0.25;
    base_poly mpm0_ = read_base_poly(3, "-1+x-y")*0.25;
    base_poly mmp0_ = read_base_poly(3, "-1-x+y")*0.25;
    base_poly mpp0_ = read_base_poly(3, "-1+x+y")*0.25;
 
    base_poly pp0m = read_base_poly(3, "1+x-z");
    base_poly pm0m = read_base_poly(3, "1-x-z");
    base_poly p0mm = read_base_poly(3, "1-y-z");
    base_poly p0pm = read_base_poly(3, "1+y-z");
 
    // (Bedrosian, 1992)
    p->base()[ 0] = mmm0_ * pmmm + base_rational_fraction(mmm0_ * xy, p00m); // N5
    p->base()[ 1] = base_rational_fraction(pp0m * pm0m * p0mm * 0.5, p00m);  // N6
    p->base()[ 2] = mpm0_ * ppmm - base_rational_fraction(mpm0_ * xy, p00m); // N7
    p->base()[ 3] = base_rational_fraction(p0pm * p0mm * pm0m * 0.5, p00m);  // N4
    p->base()[ 4] = base_rational_fraction(p0pm * p0mm * pp0m * 0.5, p00m);  // N8
    p->base()[ 5] = mmp0_ * pmpm - base_rational_fraction(mmp0_ * xy, p00m); // N3
    p->base()[ 6] = base_rational_fraction(pp0m * pm0m * p0pm * 0.5, p00m);  // N2
    p->base()[ 7] = mpp0_ * pppm + base_rational_fraction(mpp0_ * xy, p00m); // N1
    p->base()[ 8] = base_rational_fraction(pm0m * p0mm * z, p00m);           // N11
    p->base()[ 9] = base_rational_fraction(pp0m * p0mm * z, p00m);           // N12
    p->base()[10] = base_rational_fraction(pm0m * p0pm * z, p00m);           // N10
    p->base()[11] = base_rational_fraction(pp0m * p0pm * z, p00m);           // N9
    p->base()[12] = read_base_poly(3, "2*z^2-z");                            // N13
 
    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));
    p->add_node(lag_dof, base_small_vector( 1.0, -1.0, 0.0));
    p->add_node(lag_dof, base_small_vector(-1.0,  0.0, 0.0));
    p->add_node(lag_dof, base_small_vector( 1.0,  0.0, 0.0));
    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));
    p->add_node(lag_dof, base_small_vector( 1.0,  1.0, 0.0));
    p->add_node(lag_dof, base_small_vector(-0.5, -0.5, 0.5));
    p->add_node(lag_dof, base_small_vector( 0.5, -0.5, 0.5));
    p->add_node(lag_dof, base_small_vector(-0.5,  0.5, 0.5));
    p->add_node(lag_dof, base_small_vector( 0.5,  0.5, 0.5));
    p->add_node(lag_dof, base_small_vector( 0.0,  0.0, 1.0));
 
    return pfem(p);
  }
 
 
  static pfem pyramid_Q2_incomplete_fem
  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps) {
    GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
    pfem p = build_pyramid_Q2_incomplete_fem(false);
    deps.push_back(p->ref_convex(0));
    deps.push_back(p->node_tab(0));
    return p;
  }
 
  static pfem pyramid_Q2_incomplete_disc_fem
  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps) {
    GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
    pfem p = build_pyramid_Q2_incomplete_fem(true);
    deps.push_back(p->ref_convex(0));
    deps.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*    Incomplete quadratic Lagrange prism element.                      */
  /* ******************************************************************** */
 
  // local dof numeration:
  //
  //    14
  //    /|`
  //  12 | 13
  //  /  8  `
  // 9--10--11
  // |   |   |
  // |   5   |
  // 6  / `  7
  // | 3   4 |
  // |/     `|
  // 0---1---2
 
  static pfem build_prism_incomplete_P2_fem(bool disc) {
    auto p = std::make_shared<fem<base_rational_fraction>>();
    p->mref_convex() = bgeot::prism_of_reference(3);
    p->dim() = 3;
    p->is_standard() = p->is_equivalent() = true;
    p->is_polynomial() = false;
    p->is_lagrange() = true;
    p->estimated_degree() = 2;
    p->init_cvs_node();
    auto lag_dof = disc ? lagrange_nonconforming_dof(3) : lagrange_dof(3);
 
    p->base().resize(15);
 
    std::stringstream s
      ( "-2*y*z^2-2*x*z^2+2*z^2-2*y^2*z-4*x*y*z+5*y*z-2*x^2*z+5*x*z"
          "-3*z+2*y^2+4*x*y-3*y+2*x^2-3*x+1;"
        "4*(x*y*z+x^2*z-x*z-x*y-x^2+x);"
        "2*x*z^2-2*x^2*z-x*z+2*x^2-x;"
        "4*(y^2*z+x*y*z-y*z-y^2-x*y+y);"
        "4*(x*y-x*y*z);"
        "2*y*z^2-2*y^2*z-y*z+2*y^2-y;"
        "4*(y*z^2+x*z^2-z^2-y*z-x*z+z);"
        "4*(x*z-x*z^2);"
        "4*(y*z-y*z^2);"
        "-2*y*z^2-2*x*z^2+2*z^2+2*y^2*z+4*x*y*z-y*z+2*x^2*z-x*z-z;"
        "4*(-x*y*z-x^2*z+x*z);"
        "2*x*z^2+2*x^2*z-3*x*z;"
        "4*(-y^2*z-x*y*z+y*z);"
        "4*x*y*z;"
        "2*y*z^2+2*y^2*z-3*y*z;");
 
    for (int i = 0; i < 15; ++i)
      p->base()[i] = read_base_poly(3, s);
 
    p->add_node(lag_dof, base_small_vector(0.0, 0.0, 0.0));
    p->add_node(lag_dof, base_small_vector(0.5, 0.0, 0.0));
    p->add_node(lag_dof, base_small_vector(1.0, 0.0, 0.0));
    p->add_node(lag_dof, base_small_vector(0.0, 0.5, 0.0));
    p->add_node(lag_dof, base_small_vector(0.5, 0.5, 0.0));
    p->add_node(lag_dof, base_small_vector(0.0, 1.0, 0.0));
    p->add_node(lag_dof, base_small_vector(0.0, 0.0, 0.5));
    p->add_node(lag_dof, base_small_vector(1.0, 0.0, 0.5));
    p->add_node(lag_dof, base_small_vector(0.0, 1.0, 0.5));
    p->add_node(lag_dof, base_small_vector(0.0, 0.0, 1.0));
    p->add_node(lag_dof, base_small_vector(0.5, 0.0, 1.0));
    p->add_node(lag_dof, base_small_vector(1.0, 0.0, 1.0));
    p->add_node(lag_dof, base_small_vector(0.0, 0.5, 1.0));
    p->add_node(lag_dof, base_small_vector(0.5, 0.5, 1.0));
    p->add_node(lag_dof, base_small_vector(0.0, 1.0, 1.0));
 
    return pfem(p);
  }
 
 
  static pfem prism_incomplete_P2_fem
  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps) {
    GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
    pfem p = build_prism_incomplete_P2_fem(false);
    deps.push_back(p->ref_convex(0));
    deps.push_back(p->node_tab(0));
    return p;
  }
 
  static pfem prism_incomplete_P2_disc_fem
  (fem_param_list &params, std::vector<dal::pstatic_stored_object> &deps) {
    GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
    pfem p = build_prism_incomplete_P2_fem(true);
    deps.push_back(p->ref_convex(0));
    deps.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*        P1 element with a bubble base function on a face              */
  /* ******************************************************************** */
 
  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;
    base_node pt(nc); pt.fill(0.5);
    unfreeze_cvs_node();
    add_node(bubble1_dof(nc), pt);
    base_.resize(nb_dof(0));
    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
    // modifier qlq chose (transformer les fct de base P1)
  }
 
  static pfem P1_with_bubble_on_a_face(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    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 n = int(::floor(params[0].num() + 0.01));
    GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
                "Bad parameter");
    pfem p  = std::make_shared<P1_wabbfoaf_>(dim_type(n));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*    Element Pk on a square with internal nodes (with possible node    */
  /*  correspondance for surface in 3D (build for HHO elements).          */
  /* ******************************************************************** */
  // Not tested. To be tested
 
  struct CIPK_SQUARE_ : public fem<base_poly> {
    dim_type k;   // 
    mutable base_matrix K;
    base_small_vector norient;
    mutable bgeot::pgeotrans_precomp pgp;
    mutable bgeot::pgeometric_trans pgt_stored;
    // mutable pfem_precomp pfp;
 
    bgeot::pstored_point_tab pC;
 
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    CIPK_SQUARE_(dim_type nc_);
  };
 
  void CIPK_SQUARE_::mat_trans(base_matrix &M,
                              const base_matrix &G,
                              bgeot::pgeometric_trans pgt) const {
    // All the dof of this element are concentrated at the center of the element
    // This is a PK element on a square. The base functions are the monomials
    // The trans_mat is here to eliminate a potential rotation of angle k PI/2
    // The idea is to compute the gradient at the center and to sort and
    // orientate the derivative in the two directions, and exchange the base
    // functions accordingly.
    dim_type N = dim_type(G.nrows());
    gmm::copy(gmm::identity_matrix(), M);
    if (pgt != pgt_stored) {
      pgt_stored = pgt;
      pgp = bgeot::geotrans_precomp(pgt, pC, 0);
    }
    gmm::resize(K, N, 2);
    gmm::mult(G, pgp->grad(0), K);
    scalar_type s0(0), m0(M_PI);
    for (size_type i = 0; i < N; ++i, m0*=M_PI) s0 += m0 * K(i, 0);
    scalar_type s1(0), m1(M_PI);
    for (size_type i = 0; i < N; ++i, m1*=M_PI) s1 += m1 * K(i, 1);
    scalar_type a0 = gmm::sgn(s0), a1 = gmm::sgn(s1);
 
    bool inv = false;
    for (size_type i = 0; i < N; ++i) {
      if (K(i, 0) * a0 < K(i, 1) * a1 - 1e-6) break;
      if (K(i, 0) * a0 > K(i, 1) * a1 + 1e-6) { inv = true; break; }
    }
 
    if (a0 < 0.) {
      for (size_type i = 1, l = 1; i < k; ++i)
        for (size_type j = 0; j <= i; ++j, ++l)
          { if (((i-j) % 2) == 1) M(l, l) = -1.0; } // X^(i-j) Y^j
    }
 
    if (a1 < 0.) {
      for (size_type i = 1, l = 1; i < k; ++i)
        for (size_type j = 0; j <= i; ++j, ++l)
          { if ((j % 2) == 1) M(l, l) = -1.0; } // X^(i-j) Y^j
    }
 
    if (inv) {
      // std::swap(a0, a1);
      for (size_type i = 1, l = 1; i < k; ++i)
        for (size_type j = 0; j <= i; ++j, ++l)
          { M(l, l) = 0.0; M(l, l+i-2*j) = 1.0; }
    }
  }
 
  CIPK_SQUARE_::CIPK_SQUARE_(dim_type k_) {
    k = k_;
    pgt_stored = 0;
 
    cvr = bgeot::parallelepiped_of_reference(2);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = k;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
    base_.resize((k+1)*(k+2)/2);
 
    std::vector<base_node> C(1);
    C[0] = base_node(0.5, 0.5);
    pC = bgeot::store_point_tab(C);
 
    base_poly X, Y;
    read_poly(X, 2, "x-0.5"); read_poly(Y, 2, "y-0.5");
 
    base_[0] = bgeot::one_poly(2);
    add_node(ipk_center_dof(2,0), C[0]);
 
    for (size_type i = 1, l = 1; i < k; ++i)
      for (size_type j = 0; j <= i; ++j, ++l) { // X^(i-j) Y^j
        base_[l] = base_[0];
        for (size_type n = 0; n < i-j; ++n) base_[l] *= X;
        for (size_type n = 0; n < j; ++n) base_[l] *= Y;
        
        add_node(ipk_center_dof(2,l), C[0]);
      }
  }
 
  static pfem CIPK_SQUARE(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    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));
    GMM_ASSERT1(k >= 0 && k < 50 && double(k) == params[0].num(),
                "Bad parameter");
    pfem p = std::make_shared<CIPK_SQUARE_>(dim_type(k));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
 
 
 
 
  /* ******************************************************************** */
  /*    Element RTk on simplices.                                         */
  /* ******************************************************************** */
 
  struct RTk_ : public fem<base_poly> {
    dim_type nc, k;
    mutable base_matrix K;
    base_small_vector norient;
    mutable bgeot::pgeotrans_precomp pgp;
    mutable bgeot::pgeometric_trans pgt_stored;
    // mutable pfem_precomp pfp;
 
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    RTk_(dim_type nc_, dim_type k_);
  };
 
  void RTk_::mat_trans(base_matrix &M,
                       const base_matrix &G,
                       bgeot::pgeometric_trans pgt) const {
    dim_type N = dim_type(G.nrows());
    gmm::copy(gmm::identity_matrix(), M);
    if (pgt != pgt_stored) {
      pgt_stored = pgt;
      pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
      // pfp = fem_precomp(this, node_tab(0), 0);
    }
    GMM_ASSERT1(N == nc, "Sorry, this element works only in dimension " << nc);
 
    gmm::mult(G, pgp->grad(0), K); gmm::lu_inverse(K);
    for (unsigned j = 0; j <= nb_dof(0); ++j) {
      
      if (faces_of_dof(0, j).size()) {
        unsigned nf = faces_of_dof(0, j)[0];
        if (!(pgt->is_linear()))
        { gmm::mult(G, pgp->grad(j), K); gmm::lu_inverse(K); }
        bgeot::base_small_vector n(nc);
        gmm::mult(gmm::transposed(K), cvr->normals()[nf], n);
 
        M(j,j) = gmm::vect_norm2(n);
        n /= M(j,j);
        scalar_type ps = gmm::vect_sp(n, norient);
 
        if (ps < 0) M(j, j) *= scalar_type(-1);
        if (gmm::abs(ps) < 1E-8)
          GMM_WARNING2("RTk : The normal orientation may be incorrect");
      }
    }
  }
 
  // Raviart-Thomas Element on simplices, dimension d and degree k
  // The dof nodes are on the Lagrange lattice of a Lagrange P_(k+2)
  // The internal nodes (i.e. a P_(k-1) lattice) correspond to vectorial
  // Lagrange dof and the nodes on the boundary which are one oly one face
  // (i.e. internal to one face) correspond to a normal component dof.
  // The orientation of the normal component dof is selected in comparaison to
  // a the fixed direction (pi, pi^2, pi^3) arbitrarily chosen.
  // The number of ddl is (k+d+1)*(k+d-1)!/(k! (d-1)!)
  // The definition of RTk is
  // RTk = (Pk)^d + x(Pk)
  // But this is not a direct sum, in fact a direct sum is given by
  // RTk = (Pk)^d + x(Pk \ P(k-1))
  RTk_::RTk_(dim_type nc_, dim_type k_) {
    nc = nc_; k = k_;
    pgt_stored = 0;
    gmm::resize(K, nc, nc);
    gmm::resize(norient, nc);
    norient[0] = M_PI;
    for (unsigned i = 1; i < nc; ++i) norient[i] = norient[i-1]*M_PI;
 
    cvr = bgeot::simplex_of_reference(nc);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = k;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
    ntarget_dim = nc;
    vtype = VECTORIAL_PRIMAL_TYPE;
 
    bgeot::pconvex_ref cvn
      = bgeot::simplex_of_reference(nc, bgeot::short_type(k+2));
 
    size_type nddl = (k+nc+1);
    for (unsigned i = 1; i < nc; ++i) nddl *= (k+i);
    for (unsigned i = 1; i < nc; ++i) nddl /= i;
    std::vector<bgeot::base_poly> base_RT(nddl*nc, base_poly(nc,0));
 
    // Building a simple base of RTk 
    PK_fem_ PK(nc, k); // First (PK)^d
    size_type nddl_pk = (PK.base()).size();
    for (unsigned i = 0; i < nddl_pk; ++i)
      for (unsigned j = 0; j < nc; ++j) {
        base_RT[i*nc + j * (nddl_pk*nc+1)] = PK.base()[i];
      }
    unsigned l = 0;
    bgeot::power_index pi(nc);
    do { // Second, searching for degree k monomials
      base_poly p0(nc,0); p0.add_monomial(1., pi);
      
      if (pi.degree() == k) {
        for (dim_type j = 0; j < nc; ++j) {
          base_poly p(nc,0); p.add_monomial(1., pi);
          base_RT[nddl_pk*nc*nc + l * nc + j] = p * base_poly(nc, 1, j);
        }
        ++l;
      }
      for (dim_type j = 0; j < nc; ++j)
        { pi[j]++; if (pi[j] == k+1) pi[j] = 0; else break; }
    } while(pi.degree() != 0);
    GMM_ASSERT2(nddl_pk*nc+l == nddl, "Internal error");
 
    // Estimating the base functions on the dofs
    base_matrix A(nddl, nddl);
    unsigned ipoint = 0;
    for (unsigned i = 0; i < cvn->nb_points(); ++i) {
      l = 0; unsigned cpt_found = 0;
      for(dim_type j = 0; j < nc+1; ++j)
        if (gmm::abs(cvn->is_in_face(j, cvn->points()[i])) < 1E-10)
          { l = j; cpt_found++; }
      
      switch (cpt_found) {
      case 0 :
        for (unsigned j = 0; j < nddl; ++j) {
          for (unsigned m = 0; m < nc; ++m)
            A(ipoint+m, j) = base_RT[j*nc+m].eval(cvn->points()[i].begin());
        }
        for (dim_type m = 0; m < nc; ++m)
          add_node(to_coord_dof(lagrange_dof(nc), m), cvn->points()[i]);
        
        ipoint += nc; break;
      case 1 :
         for (unsigned j = 0; j < nddl; ++j) {
          base_small_vector v(nc);
          for (unsigned m = 0; m < nc; ++m)
            v[m] = base_RT[j*nc+m].eval(cvn->points()[i].begin());
          A(ipoint, j) = gmm::vect_sp(v, cvn->normals()[l]);          
        }
        add_node(normal_component_dof(nc), cvn->points()[i]);
        ++ipoint; break;
      default : break;
      }
    }
    GMM_ASSERT2(ipoint == nddl, "Internal error");
 
    // Deducing the base of shape functions
    gmm::lu_inverse(A);
    base_.resize(nddl*nc, base_poly(nc,0));
    for (size_type i = 0; i < nddl; ++i)
      for (size_type j = 0; j < nddl; ++j)
        for (unsigned m = 0; m < nc; ++m)
          if (gmm::abs(A(j, i)) > 1e-14)
            base_[i+m*nddl] += base_RT[j*nc+m] * A(j, i);
  }
 
  static pfem RTk(fem_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() == 0, "Bad type of parameters");
    GMM_ASSERT1(params[1].type() == 0, "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 && double(n) == params[0].num(),
                "Bad parameter n");
    GMM_ASSERT1(k >= 0 && k < 100 && double(k) == params[1].num(),
                "Bad parameter k");
    pfem p = std::make_shared<RTk_>(dim_type(n), dim_type(k));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
  
 
 
  // struct P1_RT0_ : public fem<base_poly> {
  //   dim_type nc;
  //   mutable base_matrix K;
  //   base_small_vector norient;
  //   mutable bgeot::pgeotrans_precomp pgp;
  //   mutable bgeot::pgeometric_trans pgt_stored;
  //   // mutable pfem_precomp pfp;
 
  //   virtual void mat_trans(base_matrix &M, const base_matrix &G,
  //                          bgeot::pgeometric_trans pgt) const;
  //   P1_RT0_(dim_type nc_);
  // };
 
  // void P1_RT0_::mat_trans(base_matrix &M,
  //                         const base_matrix &G,
  //                         bgeot::pgeometric_trans pgt) const {
  //   dim_type N = dim_type(G.nrows());
  //   gmm::copy(gmm::identity_matrix(), M);
  //   if (pgt != pgt_stored) {
  //     pgt_stored = pgt;
  //     pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
  //     // pfp = fem_precomp(this, node_tab(0), 0);
  //   }
  //   GMM_ASSERT1(N == nc, "Sorry, this element works only in dimension " << nc);
 
  //   gmm::mult(G, pgp->grad(0), K); gmm::lu_inverse(K);
  //   for (unsigned i = 0; i <= nc; ++i) {
  //     if (!(pgt->is_linear()))
  //       { gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
  //     bgeot::base_small_vector n(nc);
  //     gmm::mult(gmm::transposed(K), cvr->normals()[i], n);
 
  //     M(i,i) = gmm::vect_norm2(n);
  //     n /= M(i,i);
  //     scalar_type ps = gmm::vect_sp(n, norient);
  //     if (ps < 0) M(i, i) *= scalar_type(-1);
  //     if (gmm::abs(ps) < 1E-8)
  //       GMM_WARNING2("RT0 : The normal orientation may be incorrect");
  //   }
  // }
 
  // // The dof of this RT0 are not the integral on the edges or faces of the
  // // normal component but the normal component on the midpoint of each edge/face
  // // The reason : easier to deal in case of nonlinear transformation (otherwise
  // // an integral with several Gauss points would be necessary to compute on
  // // edges / faces)
  // // Shape fonctions on the reference element for nc = 2
  // // sqrt(2)(x, y) ; (x-1, y) ; (x, y-1)
  // // The shape functions on the real element are K \phi ||K^{-T}n_i||, where
  // //   K is the gradient of the transformation. 
  // P1_RT0_::P1_RT0_(dim_type nc_) {
  //   nc = nc_;
  //   pgt_stored = 0;
  //   gmm::resize(K, nc, nc);
  //   gmm::resize(norient, nc);
  //   norient[0] = M_PI;
  //   for (unsigned i = 1; i < nc; ++i) norient[i] = norient[i-1]*M_PI;
 
  //   cvr = bgeot::simplex_of_reference(nc);
  //   dim_ = cvr->structure()->dim();
  //   init_cvs_node();
  //   es_degree = 1;
  //   is_pol = true;
  //   is_standard_fem = is_lag = is_equiv = false;
  //   ntarget_dim = nc;
  //   vtype = VECTORIAL_PRIMAL_TYPE;
  //   base_.resize(nc*(nc+1));
 
 
  //   for (size_type j = 0; j < nc; ++j)
  //     for (size_type i = 0; i <= nc; ++i) {
  //       base_[i+j*(nc+1)] = base_poly(nc, 1, short_type(j));
  //       if (i-1 == j) base_[i+j*(nc+1)] -= bgeot::one_poly(nc);
  //       if (i == 0) base_[i+j*(nc+1)] *= sqrt(opt_long_scalar_type(nc));
  //     }
 
  //   base_node pt(nc);
  //   pt.fill(scalar_type(1)/scalar_type(nc));
  //   add_node(normal_component_dof(nc), pt);
 
  //   for (size_type i = 0; i < nc; ++i) {
  //     pt[i] = scalar_type(0);
  //     add_node(normal_component_dof(nc), pt);
  //     pt[i] = scalar_type(1)/scalar_type(nc);
      
  //   }
  // }
 
  static pfem P1_RT0(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &) {
    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 n = int(::floor(params[0].num() + 0.01));
    GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
                "Bad parameter");
    std::stringstream name;
    name << "FEM_RTK(" << n << ",0)";
    return fem_descriptor(name.str());
    
    // pfem p = std::make_shared<P1_RT0_>(dim_type(n));
    // dependencies.push_back(p->ref_convex(0));
    // dependencies.push_back(p->node_tab(0));
    // return p;
  }
 
 
  /* ******************************************************************** */
  /*    Element RT0 on parallelepideds.                                   */
  /* ******************************************************************** */
 
  struct P1_RT0Q_ : public fem<base_poly> {
    dim_type nc;
    mutable base_matrix K;
    base_small_vector norient;
    mutable bgeot::pgeotrans_precomp pgp;
    mutable bgeot::pgeometric_trans pgt_stored;
    // mutable pfem_precomp pfp;
 
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    P1_RT0Q_(dim_type nc_);
  };
 
  void P1_RT0Q_::mat_trans(base_matrix &M,
                          const base_matrix &G,
                          bgeot::pgeometric_trans pgt) const {
    dim_type N = dim_type(G.nrows());
    gmm::copy(gmm::identity_matrix(), M);
    if (pgt != pgt_stored) {
      pgt_stored = pgt;
      pgp = bgeot::geotrans_precomp(pgt, node_tab(0),0);
      // pfp = fem_precomp(this, node_tab(0), 0);
    }
    GMM_ASSERT1(N == nc, "Sorry, this element works only in dimension " << nc);
 
    gmm::mult(G, pgp->grad(0), K); gmm::lu_inverse(K);
    for (unsigned i = 0; i < unsigned(2*nc); ++i) {
      if (!(pgt->is_linear()))
        { gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
      bgeot::base_small_vector n(nc);
      gmm::mult(gmm::transposed(K), cvr->normals()[i], n);
 
      M(i,i) = gmm::vect_norm2(n);
      n /= M(i,i);
      scalar_type ps = gmm::vect_sp(n, norient);
      if (ps < 0) M(i, i) *= scalar_type(-1);
      if (gmm::abs(ps) < 1E-8)
        GMM_WARNING2("RT0Q : The normal orientation may be incorrect");
    }
  }
 
  P1_RT0Q_::P1_RT0Q_(dim_type nc_) {
    nc = nc_;
    pgt_stored = 0;
    gmm::resize(K, nc, nc);
    gmm::resize(norient, nc);
    norient[0] = M_PI;
    for (unsigned i = 1; i < nc; ++i) norient[i] = norient[i-1]*M_PI;
 
    cvr = bgeot::parallelepiped_of_reference(nc);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = 1;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
    ntarget_dim = nc;
    vtype = VECTORIAL_PRIMAL_TYPE;
    base_.resize(nc*2*nc);
 
    for (size_type j = 0; j < size_type(nc*2*nc); ++j)
      base_[j] = bgeot::null_poly(nc);
 
    for (short_type i = 0; i < nc; ++i) {
      base_[2*i+i*2*nc] = base_poly(nc, 1, i);
      base_[2*i+1+i*2*nc] = base_poly(nc, 1, i) - bgeot::one_poly(nc);
    }
 
    base_node pt(nc); pt.fill(0.5);
 
    for (size_type i = 0; i < nc; ++i) {
      pt[i] = 1.0;
      add_node(normal_component_dof(nc), pt);
      pt[i] = 0.0;
      add_node(normal_component_dof(nc), pt);
      pt[i] = 0.5;
    }
  }
 
  static pfem P1_RT0Q(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    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 n = int(::floor(params[0].num() + 0.01));
    GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
                "Bad parameter");
    pfem p = std::make_shared<P1_RT0Q_>(dim_type(n));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
 
  /* ******************************************************************** */
  /*    Element BDMk on simplices.                                        */
  /* ******************************************************************** */
 
  struct BDMk_ : public fem<base_poly> {
    dim_type nc, k;
    mutable base_matrix K;
    base_small_vector norient;
    mutable bgeot::pgeotrans_precomp pgp;
    mutable bgeot::pgeometric_trans pgt_stored;
    // mutable pfem_precomp pfp;
 
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    BDMk_(dim_type nc_, dim_type k_);
  };
 
  void BDMk_::mat_trans(base_matrix &M,
                       const base_matrix &G,
                       bgeot::pgeometric_trans pgt) const {
    dim_type N = dim_type(G.nrows());
    gmm::copy(gmm::identity_matrix(), M);
    if (pgt != pgt_stored) {
      pgt_stored = pgt;
      pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
      // pfp = fem_precomp(this, node_tab(0), 0);
    }
    GMM_ASSERT1(N == nc, "Sorry, this element works only in dimension " << nc);
 
    gmm::mult(G, pgp->grad(0), K); gmm::lu_inverse(K);
    for (unsigned j = 0; j <= nb_dof(0); ++j) {
      
      if (faces_of_dof(0, j).size()) {
        unsigned nf = faces_of_dof(0, j)[0];
        if (!(pgt->is_linear()))
        { gmm::mult(G, pgp->grad(j), K); gmm::lu_inverse(K); }
        bgeot::base_small_vector n(nc);
        gmm::mult(gmm::transposed(K), cvr->normals()[nf], n);
 
        M(j,j) = gmm::vect_norm2(n);
        n /= M(j,j);
        scalar_type ps = gmm::vect_sp(n, norient);
 
        if (ps < 0) M(j, j) *= scalar_type(-1);
        if (gmm::abs(ps) < 1E-8)
          GMM_WARNING2("RTk : The normal orientation may be incorrect");
      }
    }
  }
 
  // Brezzi-Douglas-Marini Element on simplices, dimension d and degree k=1 or 2
  // The orientation of the normal component dof is selected in comparison to
  // the fixed direction (pi, pi^2, pi^3) arbitrarily chosen.
  // The definition of RTk is
  // BDMk = (Pk)^d
 
  BDMk_::BDMk_(dim_type nc_, dim_type k_) {
    nc = nc_; k = k_;
    pgt_stored = 0;
    gmm::resize(K, nc, nc);
    gmm::resize(norient, nc);
    norient[0] = M_PI;
    for (unsigned i = 1; i < nc; ++i) norient[i] = norient[i-1]*M_PI;
 
    cvr = bgeot::simplex_of_reference(nc);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = k;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
    ntarget_dim = nc;
    vtype = VECTORIAL_PRIMAL_TYPE;
 
    bgeot::pconvex_ref cvn
      = bgeot::simplex_of_reference(nc, bgeot::short_type(k+2));
 
    // Building a simple base of BDMk 
    PK_fem_ PK(nc, k);
    size_type nddl_pk = (PK.base()).size(), nddl = nddl_pk * nc;
    std::vector<bgeot::base_poly> base_BDM(nddl*nc, base_poly(nc,0));
    for (unsigned i = 0; i < nddl_pk; ++i)
      for (unsigned j = 0; j < nc; ++j) {
        base_BDM[i*nc + j * (nddl_pk*nc+1)] = PK.base()[i];
      }
    GMM_ASSERT2(nddl_pk*nc == nddl, "Internal error");
 
    // Estimating the base functions on the dofs
    base_matrix A(nddl, nddl);
    unsigned ipoint = 0;
    for (unsigned i = 0; i < cvn->nb_points(); ++i) {
      unsigned l = 0; unsigned cpt_found = 0;
      for(dim_type j = 0; j < nc+1; ++j)
        if (gmm::abs(cvn->is_in_face(j, cvn->points()[i])) < 1E-10)
          { l = j; cpt_found++; }
      
      if (cpt_found == 1) { 
        for (unsigned j = 0; j < nddl; ++j) {
          base_small_vector v(nc);
          for (unsigned m = 0; m < nc; ++m)
            v[m] = base_BDM[j*nc+m].eval(cvn->points()[i].begin());
          A(ipoint, j) = gmm::vect_sp(v, cvn->normals()[l]);          
        }
        add_node(normal_component_dof(nc), cvn->points()[i]);
        ++ipoint;
      }
    }
 
    base_node barycenter(nc); barycenter.fill(1./(nc+1));
    
    if (k > 1) {
      bgeot::pbasic_mesh pm;
      bgeot::pmesh_precomposite pmp;
      
      structured_mesh_for_convex(PK.ref_convex(0), 1, pm, pmp);
      mesh m(*pm);
      mesh_fem mf(m);
      mf.set_classical_finite_element(pm->convex_index(), bgeot::dim_type(k-1));
      mesh_fem mfd(m);
      mfd.set_qdim(nc);
      mfd.set_classical_finite_element(pm->convex_index(), k);
      mesh_im mim(m);
      mim.set_integration_method(bgeot::dim_type(k*(k-1)));
      
      gmm::sub_interval Iu(0, mfd.nb_dof()), Ip(Iu.last(), mf.nb_dof());
      base_vector u(mfd.nb_dof()), p(mf.nb_dof());
      ga_workspace workspace;
      workspace.add_fem_variable("u", mfd, Iu, u);
      workspace.add_fem_variable("p", mf, Ip, p);
      workspace.add_expression("Div(u).Test_p", mim);
      workspace.assembly(2);
      gmm::sub_interval IA1(ipoint, mf.nb_dof()), IA2(0, nddl);
      if (gmm::mat_nrows(workspace.assembled_matrix()))
        gmm::add(gmm::sub_matrix(workspace.assembled_matrix(), Ip, Iu),
                 gmm::sub_matrix(A, IA1, IA2));
      
      for (size_type i = 0; i < mf.nb_dof(); ++i) {
        add_node(ipk_center_dof(nc, ipoint), barycenter);
        ++ipoint;
      }
    }
 
    if (k > 2) { // Completing the base with an orthogonal base of the kernel
#     if defined(GMM_USES_LAPACK)
      base_vector EIG(nddl);
      base_matrix U(nddl, nddl), V(nddl, nddl), AA(A);
      gmm::svd(AA, U, V, EIG);
      gmm::clean(V, 1E-14);
 
      for (size_type i = 0; i < nddl; ++i)
        if (gmm::abs(EIG[i]) < 1E-16) {
          gmm::copy(gmm::mat_row(V, i), gmm::mat_row(A, ipoint));
          add_node(ipk_center_dof(nc, ipoint), barycenter); ++ipoint;
        }
#     else
      GMM_ASSERT2(false, "You need Lapack to useget BDMk with k > 2");
#     endif
    }
    
    GMM_ASSERT2(ipoint == nddl, "Internal error");
 
    // Deducing the base of shape functions
    gmm::lu_inverse(A);
    base_.resize(nddl*nc, base_poly(nc,0));
    for (size_type i = 0; i < nddl; ++i)
      for (size_type j = 0; j < nddl; ++j)
        for (unsigned m = 0; m < nc; ++m)
          if (gmm::abs(A(j, i)) > 1e-14)
            base_[i+m*nddl] += base_BDM[j*nc+m] * A(j, i);
 
    // for (size_type i = 0; i < nddl; ++i)
    //   for (unsigned m = 0; m < nc; ++m)
    //     cout << base_[i+m*nddl] << endl;
  }
 
  static pfem BDMk(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    GMM_ASSERT1(params.size() == 2, "Bad number of parameters for BDM element: "
                << params.size() << " should be 2.");
    GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
    GMM_ASSERT1(params[1].type() == 0, "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 && double(n) == params[0].num(),
                "Bad parameter n for BDM element");
    GMM_ASSERT1(k >= 1 && k < 100 && double(k) == params[1].num(),
                "Bad parameter k for BDM element");
    pfem p = std::make_shared<BDMk_>(dim_type(n), dim_type(k));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
 
  /* ******************************************************************** */
  /*    Nedelec Element.                                                  */
  /* ******************************************************************** */
 
  struct P1_nedelec_ : public fem<base_poly> {
    dim_type nc;
    base_small_vector norient;
    std::vector<base_small_vector> tangents;
    mutable bgeot::pgeotrans_precomp pgp;
    mutable bgeot::pgeometric_trans pgt_stored;
    mutable pfem_precomp pfp;
 
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    P1_nedelec_(dim_type nc_);
  };
 
  void P1_nedelec_::mat_trans(base_matrix &M, const base_matrix &G,
                              bgeot::pgeometric_trans pgt) const {
    bgeot::base_small_vector t(nc), v(nc);
    GMM_ASSERT1(G.nrows() == nc,
                "Sorry, this element works only in dimension " << nc);
 
    if (pgt != pgt_stored) {
      pgt_stored = pgt;
      pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
      pfp = fem_precomp(std::make_shared<P1_nedelec_>(nc), node_tab(0), 0);
    }
    fem_interpolation_context ctx(pgp,pfp,0,G,0);
 
    for (unsigned i = 0; i < nb_dof(0); ++i) {
      ctx.set_ii(i);
      gmm::mult(ctx.K(), tangents[i], t);
      t /= gmm::vect_norm2(t);
      gmm::mult(gmm::transposed(ctx.B()), t, v);
      scalar_type ps = gmm::vect_sp(t, norient);
      if (ps < 0) v *= scalar_type(-1);
      if (gmm::abs(ps) < 1E-8)
        GMM_WARNING2("Nedelec element: "
                     "The normal orientation may be incorrect");
 
      const bgeot::base_tensor &tt = pfp->val(i);
      for (size_type j = 0; j < nb_dof(0); ++j) {
        scalar_type a = scalar_type(0);
        for (size_type k = 0; k < nc; ++k) a += tt(j, k) * v[k];
        M(j, i) = a;
      }
    }
    // In fact matrix M is diagonal (at least for linear transformations).
    // The computation can be simplified.
    gmm::lu_inverse(M);
  }
 
  P1_nedelec_::P1_nedelec_(dim_type nc_) {
    nc = nc_;
    pgt_stored = 0;
    gmm::resize(norient, nc);
    norient[0] = M_PI;
    for (unsigned i = 1; i < nc; ++i) norient[i] = norient[i-1]*M_PI;
 
    cvr = bgeot::simplex_of_reference(nc);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = 1;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
    ntarget_dim = nc;
    vtype = VECTORIAL_DUAL_TYPE;
    base_.resize(nc*(nc+1)*nc/2);
    tangents.resize(nc*(nc+1)*nc/2);
 
    std::vector<base_poly> lambda(nc+1);
    std::vector<base_vector> grad_lambda(nc+1);
    lambda[0] = bgeot::one_poly(nc);
    gmm::resize(grad_lambda[0], nc);
    gmm::fill(grad_lambda[0], scalar_type(-1));
    for (size_type i = 1; i <= nc; ++i) {
      lambda[i] = base_poly(nc, 1, short_type(i-1));
      lambda[0] -= lambda[i];
      gmm::resize(grad_lambda[i], nc);
      grad_lambda[i][i-1] = 1;
    }
 
    size_type j = 0;
    for (size_type k = 0; k <= nc; ++k)
      for (size_type l = k+1; l <= nc; ++l, ++j) {
        for (size_type i = 0; i < nc; ++i) {
          base_[j+i*(nc*(nc+1)/2)] = lambda[k] * grad_lambda[l][i]
            - lambda[l] * grad_lambda[k][i];
        }
 
        base_node pt = (cvr->points()[k] + cvr->points()[l]) / scalar_type(2);
        add_node(edge_component_dof(nc), pt);
        tangents[j] = cvr->points()[l] - cvr->points()[k];
        tangents[j] /= gmm::vect_norm2(tangents[j]);
      }
  }
 
  static pfem P1_nedelec(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    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 n = int(::floor(params[0].num() + 0.01));
    GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
                "Bad parameter");
    pfem p = std::make_shared<P1_nedelec_>(dim_type(n));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
 
  /* ******************************************************************** */
  /*    P1 element with a bubble base function on a face : type lagrange  */
  /* ******************************************************************** */
 
  struct P1_wabbfoafla_ : public PK_fem_
  { // idem elt prec mais avec raccord lagrange. A faire en dim. quelconque ..
    P1_wabbfoafla_();
  };
 
  P1_wabbfoafla_::P1_wabbfoafla_() : PK_fem_(2, 1) {
    unfreeze_cvs_node();
    es_degree = 2;
    base_node pt(2); pt.fill(0.5);
    add_node(lagrange_dof(2), pt);
    base_.resize(nb_dof(0));
 
    read_poly(base_[0], 2, "1 - y - x");
    read_poly(base_[1], 2, "x*(1 - 2*y)");
    read_poly(base_[2], 2, "y*(1 - 2*x)");
    read_poly(base_[3], 2, "4*x*y");
  }
 
  static pfem P1_with_bubble_on_a_face_lagrange(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
    pfem p = std::make_shared<P1_wabbfoafla_>();
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
 
  /* ******************************************************************** */
  /*    PK Gauss-Lobatto element on the segment                           */
  /* ******************************************************************** */
 
  static const double fem_coef_gausslob_1[4] =
    { 1.000000000000000e+00, -1.000000000000000e+00, 0.000000000000000e-01,
      1.000000000000000e+00 };
 
  static const double fem_coef_gausslob_2[9] =
    { 1.000000000000000e+00, -3.000000000000000e+00, 2.000000000000000e+00,
      0.000000000000000e-01, 4.000000000000000e+00, -4.000000000000000e+00,
      0.000000000000000e-01, -1.000000000000000e+00, 2.000000000000000e+00 };
 
  static const double fem_coef_gausslob_3[16] =
    { 1.000000000000000e+00, -6.000000000000000e+00, 1.000000000000000e+01,
      -5.000000000000000e+00, 0.000000000000000e-01, 8.090169943749474e+00,
      -1.927050983124842e+01, 1.118033988749895e+01, 0.000000000000000e-01,
      -3.090169943749474e+00, 1.427050983124842e+01, -1.118033988749895e+01,
      0.000000000000000e-01, 1.000000000000000e+00, -5.000000000000000e+00,
      5.000000000000000e+00 };
 
  static const double fem_coef_gausslob_4[25] =
    { 1.000000000000000e+00, -1.000000000000000e+01, 3.000000000000000e+01,
      -3.500000000000000e+01, 1.400000000000000e+01, 0.000000000000000e-01,
      1.351300497744848e+01, -5.687234826567877e+01, 7.602600995489696e+01,
      -3.266666666666667e+01, 0.000000000000000e-01, -5.333333333333333e+0,
      4.266666666666667e+01, -7.466666666666667e+01, 3.733333333333333e+01,
      0.000000000000000e-01, 2.820328355884853e+00, -2.479431840098789e+01,
      5.464065671176971e+01, -3.266666666666667e+01, 0.000000000000000e-01,
      -1.000000000000000e+00, 9.000000000000000e+00, -2.100000000000000e+01,
      1.400000000000000e+01 };
 
  static const double fem_coef_gausslob_5[36] =
    { 1.000000000000000e+00, -1.500000000000000e+01, 7.000000000000000e+01,
      -1.400000000000000e+02, 1.260000000000000e+02, -4.200000000000000e+01,
      0.000000000000000e-01, 2.028283187263934e+01, -1.315819844629968e+02,
      2.996878573260126e+02, -2.884609153819731e+02, 1.000722106463179e+02,
      0.000000000000000e-01, -8.072374540610696e+00, 9.849823568607377e+01,
      -2.894574355386068e+02, 3.201990314777874e+02, -1.211674570846437e+02,
      0.000000000000000e-01, 4.489369296352334e+00, -6.035445290945900e+01,
      2.203358804738940e+02, -2.856382539454310e+02, 1.211674570846437e+02,
      0.000000000000000e-01, -2.699826628380976e+00, 3.743820168638206e+01,
      -1.465663022612998e+02, 2.119001378496167e+02, -1.000722106463179e+02,
      0.000000000000000e-01, 1.000000000000000e+00, -1.400000000000000e+01,
      5.600000000000000e+01, -8.400000000000000e+01, 4.200000000000000e+01 };
 
  static const double fem_coef_gausslob_6[49] =
    { 1.000000000000000e+00, -2.100000000000000e+01, 1.400000000000000e+02,
      -4.200000000000000e+02, 6.300000000000000e+02, -4.620000000000000e+02,
      1.320000000000000e+02, 0.000000000000000e-01, 2.840315320583963e+01,
      -2.618707992013254e+02, 8.915455952644434e+02, -1.426720320066430e+03,
      1.086906031987037e+03, -3.182636611895653e+02, 0.000000000000000e-01,
      -1.133797045109102e+01, 1.954053162272040e+02, -8.515354909154440e+02,
      1.555570646517052e+03, -1.285566162567286e+03, 3.974636611895653e+02,
      0.000000000000000e-01, 6.400000000000000e+00, -1.216000000000000e+02,
      6.528000000000000e+02, -1.382400000000000e+03, 1.267200000000000e+03,
      -4.224000000000000e+02, 0.000000000000000e-01, -4.099929626153486e+00,
      8.051601475380118e+01, -4.643586932712080e+02, 1.089694751524101e+03,
      -1.099215804570106e+03, 3.974636611895653e+02, 0.000000000000000e-01,
      2.634746871404869e+00, -5.245053177967978e+01, 3.115485889222086e+02,
      -7.661450779747230e+02, 8.226759351503546e+02, -3.182636611895653e+02,
      0.000000000000000e-01, -1.000000000000000e+00, 2.000000000000000e+01,
      -1.200000000000000e+02, 3.000000000000000e+02, -3.300000000000000e+02,
      1.320000000000000e+02 };
 
  static const double fem_coef_gausslob_7[64] =
    { 1.000000000000000e+00, -2.800000000000000e+01, 2.520000000000000e+02,
      -1.050000000000000e+03, 2.310000000000000e+03, -2.772000000000000e+03,
      1.716000000000000e+03, -4.290000000000000e+02, 0.000000000000000e-01,
      3.787519721423474e+01, -4.699045385400572e+02, 2.217166528901653e+03,
      -5.195916436930171e+03, 6.469992588404291e+03, -4.101225846986446e+03,
      1.042012507936495e+03, 0.000000000000000e-01, -1.513857963869697e+01,
      3.497259983590744e+02, -2.101837798737552e+03, 5.599947843186200e+03,
      -7.539551580657127e+03, 5.032695631593761e+03, -1.325841514105659e+03,
      0.000000000000000e-01, 8.595816328530350e+00, -2.189405837038675e+02,
      1.612357003153179e+03, -4.947308401212060e+03, 7.342604827965170e+03,
      -5.255204822337236e+03, 1.457896159806284e+03, 0.000000000000000e-01,
      -5.620377978515898e+00, 1.480753190084385e+02, -1.171440824431867e+03,
      3.964008996775197e+03, -6.427195249873722e+03, 4.950068296306753e+03,
      -1.457896159806284e+03, 0.000000000000000e-01, 3.883318851088245e+00,
      -1.038534676200790e+02, 8.481025943868683e+02, -3.011828579891082e+03,
      5.186049587313396e+03, -4.248194967145850e+03, 1.325841514105659e+03,
      0.000000000000000e-01, -2.595374776640464e+00, 6.989727249649080e+01,
      -5.793475032722815e+02, 2.106096578071916e+03, -3.744900173152008e+03,
      3.192861708569018e+03, -1.042012507936495e+03, 0.000000000000000e-01,
      1.000000000000000e+00, -2.700000000000000e+01, 2.250000000000000e+02,
      -8.250000000000000e+02, 1.485000000000000e+03, -1.287000000000000e+03,
      4.290000000000000e+02 };
 
  static const double fem_coef_gausslob_8[81] =
    { 1.000000000000000e+00, -3.600000000000000e+01, 4.200000000000000e+02,
      -2.310000000000000e+03, 6.930000000000000e+03, -1.201200000000000e+04,
      1.201200000000000e+04, -6.435000000000000e+03, 1.430000000000000e+03,
      0.000000000000000e-01, 4.869949034318613e+01, -7.815432549965997e+02,
      4.860656931912704e+03, -1.551737634456316e+04, 2.788918324236022e+04,
      -2.854121637661796e+04, 1.553203650368708e+04, -3.490440192125469e+03,
      0.000000000000000e-01, -1.947740331442309e+01, 5.805138088476765e+02,
      -4.583922431826645e+03, 1.659300238583299e+04, -3.217606831061953e+04,
      3.461498040096808e+04, -1.950464316590829e+04, 4.495614716020147e+03,
      0.000000000000000e-01, 1.108992781389876e+01, -3.644117397881921e+02,
      3.513408770397022e+03, -1.458458798126453e+04, 3.105326914181443e+04,
      -3.569574046476449e+04, 2.111700401254369e+04, -5.050031666751820e+03,
      0.000000000000000e-01, -7.314285714285714e+00, 2.486857142857143e+02,
      -2.574628571428571e+03, 1.174674285714286e+04, -2.719451428571429e+04,
      3.347017142857143e+04,  -2.091885714285714e+04, 5.229714285714286e+03,
      0.000000000000000e-01, 5.181491353118710e+00, -1.789312751409961e+02,
      1.913693930879634e+03, -9.161425477258226e+03, 2.246586272145708e+04,
      -2.927759904600966e+04, 1.928324932147088e+04, -5.050031666751820e+03,
      0.000000000000000e-01, -3.748881746893966e+00, 1.304893011815078e+02,
      -1.418925315009559e+03, 6.967886161876575e+03, -1.767073170824305e+04,
      2.395969028817416e+04, -1.646027456225289e+04, 4.495614716020147e+03,
      0.000000000000000e-01, 2.569661265399177e+00, -8.980255438911062e+01,
      9.847166850754155e+02, -4.899241601766511e+03, 1.264999919894514e+04,
      -1.754928623032154e+04, 1.239148503331668e+04, -3.490440192125469e+03,
      0.000000000000000e-01, -1.000000000000000e+00, 3.500000000000000e+01,
      -3.850000000000000e+02, 1.925000000000000e+03, -5.005000000000000e+03,
      7.007000000000000e+03, -5.005000000000000e+03, 1.430000000000000e+03 };
 
  static const double fem_coef_gausslob_9[100] =
    { 1.000000000000000e+00, -4.500000000000000e+01, 6.600000000000000e+02,
      -4.620000000000000e+03, 1.801800000000000e+04, -4.204200000000000e+04,
      6.006000000000000e+04, -5.148000000000000e+04, 2.431000000000000e+04,
      -4.862000000000000e+03, 0.000000000000000e-01, 6.087629005856384e+01,
      -1.226341297551814e+03, 9.697405499616728e+03, -4.021760400071670e+04,
      9.725280603247945e+04, -1.421240159771644e+05, 1.237105621496494e+05,
      -5.906188663911807e+04, 1.190819794274692e+04, 0.000000000000000e-01,
      -2.435589341485964e+01, 9.095415690555422e+02, -9.111255870205338e+03,
      4.276661412893713e+04, -1.114146601604227e+05, 1.709573510663859e+05,
      -1.539309822784481e+05, 7.531473186776967e+04, -1.546698442965720e+04,
      0.000000000000000e-01, 1.388757697026791e+01, -5.717395054991898e+02,
      6.975542885191629e+03, -3.743822964286359e+04, 1.068054892026842e+05,
      -1.747039036097183e+05, 1.648204809285228e+05, -8.352714678107793e+04,
      1.762561894579012e+04, 0.000000000000000e-01, -9.198709522206266e+00,
      3.919017281073292e+02, -5.132147808492991e+03, 3.020136030457121e+04,
      -9.337958558844686e+04, 1.629937209113228e+05, -1.619399017586618e+05,
      8.553993533073839e+04, -1.866608440961585e+04, 0.000000000000000e-01,
      6.589286067498368e+00, -2.852085019651226e+02, 3.859417687766574e+03,
      -2.381847798089535e+04, 7.784545414265617e+04, -1.434184925463668e+05,
      1.495994578589254e+05, -8.245482435580428e+04, 1.866608440961585e+04,
      0.000000000000000e-01, -4.905768350885374e+00, 2.141208512030290e+02,
      -2.948048350518040e+03, 1.867520721718195e+04, -6.311993447254502e+04,
      1.208313444661297e+05, -1.311255887283438e+05, 7.510342373103317e+04,
      -1.762561894579012e+04, 0.000000000000000e-01, 3.659127863806492e+00,
      -1.604518938956052e+02, 2.230466872754075e+03, -1.433960781600324e+04,
      4.943623515122416e+04, -9.697372467640532e+04, 1.082245668039501e+05,
      -6.388812799914516e+04, 1.546698442965720e+04, 0.000000000000000e-01,
      -2.551909672185327e+00, 1.121770505458310e+02, -1.567380916112636e+03,
      1.015673778978859e+04, -3.539780430762936e+04, 7.040572036581639e+04,
      -7.991059497559391e+04, 4.811189484560421e+04, -1.190819794274692e+04,
      0.000000000000000e-01, 1.000000000000000e+00, -4.400000000000000e+01,
      6.160000000000000e+02, -4.004000000000000e+03, 1.401400000000000e+04,
      -2.802800000000000e+04, 3.203200000000000e+04, -1.944800000000000e+04,
      4.862000000000000e+03 };
 
  static const double fem_coef_gausslob_10[121] =
    { 1.000000000000000e+00, -5.500000000000000e+01, 9.900000000000000e+02,
      -8.580000000000000e+03, 4.204200000000000e+04, -1.261260000000000e+05,
      2.402400000000000e+05, -2.917200000000000e+05, 2.187900000000000e+05,
      -9.237800000000000e+04, 1.679600000000000e+04, 0.000000000000000e-01,
      7.440573476392379e+01, -1.837547309495916e+03, 1.797721877249297e+04,
      -9.362519219895049e+04, 2.909773746534557e+05, -5.668103968059107e+05,
      6.987888266998443e+05, -5.297630128145374e+05, 2.254581472606030e+05,
      -4.123982399226536e+04, 0.000000000000000e-01, -2.977479259019735e+01,
      1.361302600480045e+03, -1.684411461834327e+04, 9.915381814787320e+04,
      -3.316413447923031e+05, 6.777335995676564e+05, -8.637075009663581e+05,
      6.706702925312933e+05, -2.905859066780586e+05, 5.388962900035043e+04,
      0.000000000000000e-01, 1.699123898926703e+01, -8.563552206749944e+02,
      1.288193007592445e+04, -8.652550760766087e+04, 3.163119148667456e+05,
      -6.879421746683590e+05, 9.173107022760963e+05, -7.368809317678005e+05,
      3.277210563158369e+05, -6.203762550909711e+04, 0.000000000000000e-01,
      -1.128077599537177e+01, 5.884060270969327e+02, -9.496934299975062e+03,
      6.981839702203256e+04, -2.759867860090975e+05, 6.390150586777585e+05,
      -8.953334100127128e+05, 7.481407085083284e+05, -3.434511859876541e+05,
      6.671702685021839e+04, 0.000000000000000e-01, 8.126984126984127e+00,
      -4.307301587301587e+02, 7.184253968253968e+03, -5.536101587301587e+04,
      2.309526349206349e+05, -5.631187301587302e+05, 8.261892063492063e+05,
      -7.184253968253968e+05, 3.412520634920635e+05, -6.825041269841270e+04,
      0.000000000000000e-01, -6.131078401763724e+00, 3.277460052754944e+02,
      -5.563892337052540e+03, 4.401679638240306e+04, -1.898229640674875e+05,
      4.802970424048780e+05, -7.315927845245721e+05, 6.593462428792687e+05,
      -3.237190825145298e+05, 6.671702685021839e+04, 0.000000000000000e-01,
      4.719539826177156e+00, -2.535442364309725e+02, 4.348374949258725e+03,
      -3.493745849693315e+04, 1.537770968392435e+05, -3.987659746141970e+05,
      6.242937855878344e+05, -5.790845728346388e+05, 2.926551987751342e+05,
      -6.203762550909711e+04, 0.000000000000000e-01, -3.595976071589556e+00,
      1.937484128537626e+02, -3.342868418261306e+03, 2.710687970739982e+04,
      -1.208013807254569e+05, 3.181552127960248e+05, -5.073176789159280e+05,
      4.804304374445346e+05, -2.483103833254456e+05, 5.388962900035043e+04,
      0.000000000000000e-01, 2.539125352570295e+00, -1.370261203741934e+02,
      2.372031907702074e+03, -1.933271708314826e+04, 8.675745431426523e+04,
      -2.305316371991209e+05, 3.716008535065897e+05, -3.564317671210514e+05,
      1.869400926620506e+05, -4.123982399226536e+04, 0.000000000000000e-01,
      -1.000000000000000e+00, 5.400000000000000e+01, -9.360000000000000e+02,
      7.644000000000000e+03, -3.439800000000000e+04, 9.172800000000000e+04,
      -1.485120000000000e+05, 1.432080000000000e+05, -7.558200000000000e+04,
      1.679600000000000e+04 };
 
  static const double fem_coef_gausslob_11[144] =
    { 1.000000000000000e+00, -6.600000000000000e+01, 1.430000000000000e+03,
      -1.501500000000000e+04, 9.009000000000000e+04, -3.363360000000000e+05,
      8.168160000000000e+05, -1.312740000000000e+06, 1.385670000000000e+06,
      -9.237800000000000e+05, 3.527160000000000e+05, -5.878600000000000e+04,
      0.000000000000000e-01, 8.928790437897459e+01, -2.652104227906408e+03,
      3.141784850363868e+04, -2.002791716233576e+05, 7.743819221375186e+05,
      -1.922871126014229e+06, 3.137025618338122e+06, -3.346678943928588e+06,
      2.248625248986712e+06, -8.636670995581690e+05, 1.446085194818801e+05,
      0.000000000000000e-01, -3.573451504477809e+01, 1.963011183403142e+03,
      -2.937610003978132e+04, 2.114542549822465e+05, -8.792000471014727e+05,
      2.288870840565664e+06, -3.858042797257275e+06, 4.213930780912515e+06,
      -2.881507044264936e+06, 1.121761824282757e+06, -1.898189887480762e+05,
      0.000000000000000e-01, 2.040222049137883e+01,  -1.235400296692680e+03,
      2.244501976036562e+04, -1.840644193090062e+05, 8.352985709834925e+05,
      -2.311501023251837e+06, 4.072373742578464e+06, -4.597525083798586e+06,
      3.224564284996090e+06, -1.280533828091591e+06, 2.201577342088092e+05,
      0.000000000000000e-01, -1.356395972416294e+01, 8.500434612026409e+02,
      -1.656519758544980e+04, 1.484886632288921e+05, -7.274015626850215e+05,
      2.139270125221522e+06, -3.953929242730888e+06, 4.636483776329532e+06,
      -3.352298848558998e+06, 1.364514095203156e+06, -2.393982879242230e+05,
      0.000000000000000e-01, 9.803369220897211e+00, -6.243148521745071e+02,
      1.257271925920887e+04, -1.180754347844014e+05, 6.096877374094766e+05,
      -1.885008008173690e+06, 3.641310114569300e+06, -4.434918589444567e+06,
      3.311645240511804e+06, -1.385401912847929e+06, 2.488026449837513e+05,
      0.000000000000000e-01, -7.447686911212630e+00, 4.784415903433910e+02,
      -9.808275331323619e+03, 9.456732952354300e+04, -5.045513347291155e+05,
      1.617062776783016e+06, -3.237833360324115e+06, 4.079238919323811e+06,
      -3.141771586138831e+06, 1.351427181973335e+06, -2.488026449837513e+05,
      0.000000000000000e-01, 5.819619912607334e+00, -3.757783851118444e+02,
      7.785050290867036e+03, -7.625636515182549e+04, 4.153153824802320e+05,
      -1.363841143951918e+06, 2.804561170833446e+06, -3.631789084056234e+06,
      2.874063732359706e+06, -1.268867071963298e+06, 2.393982879242230e+05,
      0.000000000000000e-01, -4.587084978600364e+00, 2.971291972095372e+02,
      -6.195597982053585e+03, 6.128651035627590e+04, -3.381847677955126e+05,
      1.128582073344179e+06, -2.364480249965060e+06, 3.125557361498120e+06,
      -2.527901385564680e+06, 1.141201248205309e+06, -2.201577342088092e+05,
      0.000000000000000e-01, 3.549738472143437e+00, -2.303803824096343e+02,
      4.822860472410032e+03, -4.799737703690675e+04, 2.670306747478879e+05,
      -9.003482951717774e+05, 1.909697516429202e+06, -2.560483668180433e+06,
      2.103933182581560e+06, -9.662470519460815e+05, 1.898189887480762e+05,
      0.000000000000000e-01, -2.529605817247381e+00, 1.643527121363626e+02,
      -3.448327347881922e+03, 3.443600981453996e+04, -1.924805754474854e+05,
      6.528637806490704e+05, -1.394862512471197e+06, 1.886334531344429e+06,
      -1.565422824908427e+06, 7.270266147425120e+05, -1.446085194818801e+05,
      0.000000000000000e-01, 1.000000000000000e+00, -6.500000000000000e+01,
      1.365000000000000e+03, -1.365000000000000e+04, 7.644000000000000e+04,
      -2.598960000000000e+05, 5.569200000000000e+05, -7.558200000000000e+05,
      6.298500000000000e+05, -2.939300000000000e+05, 5.878600000000000e+04 };
 
  static const double fem_coef_gausslob_12[169] =
    { 1.000000000000000e+00, -7.800000000000000e+01, 2.002000000000000e+03,
      -2.502500000000000e+04, 1.801800000000000e+05, -8.168160000000000e+05,
      2.450448000000000e+06, -4.988412000000000e+06, 6.928350000000000e+06,
      -6.466460000000000e+06, 3.879876000000000e+06, -1.352078000000000e+06,
      2.080120000000000e+05, 0.000000000000000e-01, 1.055228477237708e+02,
      -3.710649283521791e+03, 5.230891096204477e+04, -4.000264992878070e+05,
      1.877737871587605e+06, -5.758751500257082e+06, 1.189876582421287e+07,
      -1.670085336595664e+07, 1.570840983751716e+07, -9.480360124877962e+06,
      3.318799039873028e+06, -5.124248673374161e+05, 0.000000000000000e-01,
      -4.223530748650643e+01, 2.744602756239140e+03, -4.883026612748412e+04,
      4.213447894212839e+05, -2.125569662252139e+06, 6.831231541213443e+06,
      -1.457745185889376e+07, 2.094131826427298e+07, -2.004063010887189e+07,
      1.225627209853367e+07, -4.335340118801754e+06, 6.749529540569050e+05,
      0.000000000000000e-01, 2.412126940135046e+01, -1.727728082317750e+03,
      3.727953470399318e+04, -3.660428937429080e+05, 2.013286677665669e+06,
      -6.871455230919223e+06, 1.531439984708661e+07, -2.272429867457359e+07,
      2.229291249432404e+07, -1.390087225188044e+07, 4.993770899110666e+06,
      -7.872767949619026e+05, 0.000000000000000e-01, -1.605014152757175e+01,
      1.189832345018574e+03, -2.753035300171167e+04, 2.951729666043859e+05,
      -1.750245295939125e+06, 6.340418364700963e+06, -1.480658414760223e+07,
      2.279585234765233e+07, -2.303126179585568e+07, 1.470734589654657e+07,
      -5.387526099148482e+06, 8.631843338394749e+05, 0.000000000000000e-01,
      1.162284069412542e+01, -8.756167723813171e+02, 2.093616690467957e+04,
      -2.350848402511507e+05, 1.467905987404880e+06, -5.583024412516426e+06,
      1.360724316266146e+07, -2.172800271110029e+07, 2.264080373496021e+07,
      -1.484050586374976e+07, 5.558088637639386e+06, -9.074958680213037e+05,
      0.000000000000000e-01, -8.865800865800866e+00, 6.738008658008658e+02,
      -1.640173160173160e+04, 1.890632034632035e+05, -1.219313593073593e+06,
      4.803100813852814e+06, -1.211898237229437e+07, 1.998830268398268e+07,
      -2.144876606060606e+07, 1.443281454545455e+07, -5.532578909090909e+06,
      9.220964848484848e+05, 0.000000000000000e-01, 6.984318980773296e+00,
      -5.335955916987455e+02, 1.312836634638491e+04, -1.537652068533838e+05,
      1.012269836848560e+06, -4.084347297774684e+06, 1.057602476941974e+07,
      -1.790936242524428e+07, 1.971847079705798e+07, -1.359625813912256e+07,
      5.331861778616258e+06, -9.074958680213037e+05, 0.000000000000000e-01,
      -5.596679201904262e+00, 4.289927383042951e+02, -1.062596939367885e+04,
      1.257256555151274e+05, -8.388435073376655e+05, 3.440109149730765e+06,
      -9.074697250266149e+06, 1.567950042058767e+07, -1.762881516112804e+07,
      1.241472483931862e+07, -4.970685906925217e+06, 8.631843338394749e+05,
      0.000000000000000e-01, 4.489137849846207e+00, -3.448281543178960e+02,
      8.578250876817155e+03, -1.021659510074640e+05, 6.876731048919487e+05,
      -2.851145716716003e+06, 7.618634882796075e+06, -1.335715271315875e+07,
      1.525930546501226e+07, -1.092966082914868e+07, 4.453550640432165e+06,
      -7.872767949619026e+05, 0.000000000000000e-01, -3.514808358523940e+00,
      2.703477424140389e+02, -6.743800341396192e+03, 8.065306199056715e+04,
      -5.459331868312813e+05, 2.279586137602263e+06, -6.143562568432354e+06,
      1.087848437431968e+07, -1.256803423488744e+07, 9.114425759470104e+06,
      -3.764095329881106e+06, 6.749529540569050e+05, 0.000000000000000e-01,
      2.522322790441051e+00, -1.941585635394145e+02, 4.850890672082830e+03,
      -5.815428585185421e+04, 3.949277670351431e+05, -1.655905848916830e+06,
      4.485333711312109e+06, -7.989838200781784e+06, 9.294715032477446e+06,
      -6.793611930544113e+06, 2.830299368175965e+06, -5.124248673374161e+05,
      0.000000000000000e-01, -1.000000000000000e+00, 7.700000000000000e+01,
      -1.925000000000000e+03, 2.310000000000000e+04, -1.570800000000000e+05,
      6.597360000000000e+05, -1.790712000000000e+06, 3.197700000000000e+06,
      -3.730650000000000e+06, 2.735810000000000e+06, -1.144066000000000e+06,
      2.080120000000000e+05 };
 
  static const double fem_coef_gausslob_13[196] =
    { 1.000000000000000e+00, -9.100000000000000e+01, 2.730000000000000e+03,
      -4.004000000000000e+04, 3.403400000000000e+05, -1.837836000000000e+06,
      6.651216000000000e+06, -1.662804000000000e+07, 2.909907000000000e+07,
      -3.556553000000000e+07, 2.974571600000000e+07, -1.622493600000000e+07,
      5.200300000000000e+06, -7.429000000000000e+05, 0.000000000000000e-01,
      1.231105960095249e+02, -5.057514000718615e+03, 8.362619816879477e+04,
      -7.548172444491492e+05, 4.219784854542762e+06, -1.560990588170226e+07,
      3.960523865471669e+07, -7.003645336672149e+07, 8.625846242015506e+07,
      -7.256273938600024e+07, 3.975792117062384e+07, -1.278833059440045e+07,
      1.832147578471145e+06, 0.000000000000000e-01, -4.927732466948325e+01,
      3.738734011094227e+03, -7.796486012083982e+04, 7.935560370804071e+05,
      -4.765562636392240e+06, 1.846680880578437e+07, -4.837506802694851e+07,
      8.753277424234758e+07, -1.096660444159086e+08, 9.346798694406962e+07,
      -5.173889595224656e+07, 1.677849814552107e+07, -2.419777739872722e+06,
      0.000000000000000e-01, 2.814884111282567e+01, -2.353904695653627e+03,
      5.948276501192433e+04, -6.883051529064705e+05, 4.502895601024915e+06,
      -1.851735708867110e+07, 5.063079101609929e+07, -9.458217148027056e+07,
      1.214199817163488e+08, -1.054743067017811e+08, 5.927651305189011e+07,
      -1.946011726944643e+07, 2.834919298555227e+06, 0.000000000000000e-01,
      -1.874042032746388e+01, 1.621968976936385e+03, -4.394233902947216e+04,
      5.547892506224127e+05, -3.908876121817967e+06, 1.704431633766850e+07,
      -4.878627691196220e+07, 9.448003328782282e+07, -1.248200449226441e+08,
      1.109678546438889e+08, -6.355498649510845e+07, 2.119358718443647e+07,
      -3.128057142433586e+06, 0.000000000000000e-01, 1.358783086373605e+01,
      -1.195146716951513e+03, 3.345811195278515e+04, -4.422483366082136e+05,
      3.278781781106079e+06, -1.499532485039955e+07, 4.474689669452104e+07,
      -8.978037194074775e+07, 1.222039761714010e+08, -1.114085938050346e+08,
      6.517880226738924e+07, -2.213159774622623e+07, 3.317403211532315e+06,
      0.000000000000000e-01, -1.039065134180050e+01, 9.220321824274881e+02,
      -2.627964906979335e+04, 3.565631308351814e+05, -2.729347397416285e+06,
      1.291899991743804e+07, -3.987098284679090e+07, 8.253644504125924e+07,
      -1.155541226713629e+08, 1.080161713852185e+08, -6.460510650895925e+07,
      2.236736005147009e+07, -3.410612094153076e+06, 0.000000000000000e-01,
      8.225051804237637e+00, -7.337438600306013e+02, 2.113982918743374e+04,
      -2.914573383452468e+05, 2.277144431399071e+06, -1.103660550091489e+07,
      3.493361321847434e+07, -7.418006034176242e+07, 1.064417028079657e+08,
      -1.018292957440870e+08, 6.222452923525811e+07, -2.197059717251990e+07,
      3.410612094153076e+06, 0.000000000000000e-01,  -6.651370533890156e+00,
      5.953674398464600e+02, -1.727143617692029e+04, 2.405949101477657e+05,
      -1.905359074321061e+06, 9.386078020968711e+06, -3.025905095154839e+07,
      6.552811535463406e+07, -9.594395490329804e+07, 9.365009838355809e+07,
      -5.835707981219508e+07, 2.099464400369387e+07, -3.317403211532315e+06,
      0.000000000000000e-01, 5.430796129527214e+00, -4.871978584294892e+02,
      1.419769025501944e+04, -1.991370307462581e+05, 1.591472100521582e+06,
      -7.928247009292726e+06, 2.589562028603176e+07, -5.690356974983853e+07,
      8.463743197871011e+07, -8.398458536550263e+07, 5.322039739169053e+07,
      -1.947115566720015e+07, 3.128057142433586e+06, 0.000000000000000e-01,
      -4.414467779036191e+00, 3.966097971621681e+02, -1.159268854505275e+04,
      1.633445673328560e+05, -1.313458735263122e+06, 6.593624701202045e+06,
      -2.173390279168494e+07, 4.826160481317836e+07, -7.262663174126566e+07,
      7.298651647234048e+07, -4.687881110584063e+07, 1.739383361177152e+07,
      -2.834919298555227e+06, 0.000000000000000e-01, 3.487743182287134e+00,
      -3.136500314357888e+02, 9.185689380767778e+03, -1.298134042763866e+05,
      1.048017196494400e+06, -5.287706376855868e+06, 1.753577438278919e+07,
      -3.921741432164137e+07, 5.949694434313328e+07, -6.033542452985016e+07,
      3.913958191606598e+07, -1.467861247282431e+07, 2.419777739872722e+06,
      0.000000000000000e-01, -2.516624450464659e+00, 2.264447557529062e+02,
      -6.639311014646825e+03, 9.399061131310191e+04, -7.605959998781342e+05,
      3.848998924774717e+06, -1.281093272369737e+07, 2.877371846174006e+07,
      -4.386952078323465e+07, 4.473878170318014e+07, -2.920546515856783e+07,
      1.102958792572444e+07, -1.832147578471145e+06, 0.000000000000000e-01,
      1.000000000000000e+00, -9.000000000000000e+01, 2.640000000000000e+03,
      -3.740000000000000e+04, 3.029400000000000e+05, -1.534896000000000e+06,
      5.116320000000000e+06, -1.151172000000000e+07, 1.758735000000000e+07,
      -1.797818000000000e+07, 1.176753600000000e+07, -4.457400000000000e+06,
      7.429000000000000e+05 };
 
  static const double fem_coef_gausslob_14[225] =
    { 1.000000000000000e+00, -1.050000000000000e+02, 3.640000000000000e+03,
      -6.188000000000000e+04, 6.126120000000000e+05, -3.879876000000000e+06,
      1.662804000000000e+07, -4.988412000000000e+07, 1.066965900000000e+08,
      -1.636014380000000e+08, 1.784742960000000e+08, -1.352078000000000e+08,
      6.760390000000000e+07, -2.005830000000000e+07, 2.674440000000000e+06,
      0.000000000000000e-01, 1.420511699552723e+02, -6.740724197514907e+03,
      1.291563810650317e+05, -1.357536885890637e+06, 8.899789739175488e+06,
      -3.898284725729828e+07, 1.186784002318284e+08, -2.564866709310747e+08,
      3.962828733570607e+08, -4.348015522628213e+08, 3.308658899677545e+08,
      -1.660170000478420e+08, 4.939776003229131e+07, -6.601663651220939e+06,
      0.000000000000000e-01, -5.686066782897496e+01, 4.980782943180571e+03,
      -1.202886759366057e+05, 1.425067611910860e+06, -1.003204881513823e+07,
      4.601736532886305e+07, -1.446080878593359e+08, 3.197256124417587e+08,
      -5.024241189222256e+08, 5.584380089699472e+08, -4.292685524711295e+08,
      2.171358326952990e+08, -6.503152653250142e+07, 8.737812306213077e+06,
      0.000000000000000e-01, 3.248522677540712e+01, -3.136209432322820e+03,
      9.172216152088957e+04, -1.234458149638161e+06, 9.460578149039340e+06,
      -4.602710110035357e+07, 1.508977142152744e+08, -3.443000936718750e+08,
      5.541917935770469e+08, -6.276282328622294e+08, 4.896977188070307e+08,
      -2.507043130148833e+08, 7.583045543918066e+07, -1.027267982590785e+07,
      0.000000000000000e-01, -2.163547691954172e+01, 2.161829695966707e+03,
      -6.777232432848026e+04, 9.945601317987818e+05, -8.202377657972944e+06,
      4.227975782999525e+07, -1.449995018490026e+08, 3.427553248466447e+08,
      -5.674380093336526e+08, 6.573468625138306e+08, -5.224446315267705e+08,
      2.715766154508912e+08, -8.319461131269656e+07, 1.139164303704407e+07,
      0.000000000000000e-01, 1.569978564006007e+01, -1.594280678736617e+03,
      5.164364569735884e+04, -7.932250762767995e+05, 6.879605840139093e+06,
      -3.716431671531196e+07, 1.327627119052969e+08, -3.248633554644047e+08,
      5.536614212284358e+08, -6.572276042331931e+08, 5.332102141003966e+08,
      -2.820526436178851e+08, 8.770029066945359e+07, -1.216316370145453e+07,
      0.000000000000000e-01, -1.202518395631915e+01, 1.231993083463710e+03,
      -4.063141778456666e+04, 6.405521495953326e+05, -5.734055805513453e+06,
      3.204057307537486e+07, -1.182863821502541e+08, 2.983631937198568e+08,
      -5.225422090403253e+08, 6.354190974001179e+08, -5.265537919663995e+08,
      2.837551732890861e+08, -8.968009595208465e+07, 1.261735673043107e+07,
      0.000000000000000e-01, 9.547785547785548e+00, -9.834219114219114e+02,
      3.278709557109557e+04, -5.252427785547786e+05, 4.798602442890443e+06,
      -2.744701911421911e+07, 1.038669217715618e+08, -2.685490364568765e+08,
      4.816180870862471e+08, -5.987952711608392e+08, 5.064437616783217e+08,
      -2.780475554312354e+08, 8.937242853146853e+07, -1.276748979020979e+07,
      0.000000000000000e-01, -7.763592643695499e+00, 8.024013730981778e+02,
      -2.693903659000494e+04, 4.360799342555947e+05, -4.038452021764897e+06,
      2.347605516322762e+07, -9.046087127754393e+07, 2.384165701554799e+08,
      -4.360078989902932e+08, 5.526354677146948e+08, -4.761786531169400e+08,
      2.660933883812130e+08, -8.696289827395033e+07, 1.261735673043107e+07,
      0.000000000000000e-01, 6.402657115106499e+00, -6.632652203626004e+02,
      2.237191510124190e+04, -3.647008347515229e+05, 3.408912135873189e+06,
      -2.004238739169268e+07, 7.824760241311718e+07, -2.092325230710293e+08,
      3.885803431690498e+08, -5.004334616015021e+08, 4.381904244262927e+08,
      -2.487967617473505e+08, 8.258400115090981e+07, -1.216316370145453e+07,
      0.000000000000000e-01, -5.303584550798654e+00, 5.502727059685037e+02,
      -1.861988467344103e+04, 3.050015660700474e+05, -2.869271809771707e+06,
      1.700462338220501e+07, -6.701518639186630e+07, 1.811217810430339e+08,
      -3.403535526125257e+08, 4.438883801280693e+08, -3.938531369776322e+08,
      2.266861847568459e+08, -7.628839120592035e+07, 1.139164303704407e+07,
      0.000000000000000e-01, 4.356127432886088e+00, -4.524531153029460e+02,
      1.534317874813675e+04, -2.521565333504139e+05, 2.382646312619400e+06,
      -1.419908547341188e+07, 5.633074020272954e+07, -1.534171364617007e+08,
      2.907942363862238e+08, -3.828802350952767e+08, 3.432339697459343e+08,
      -1.997222564631490e+08, 6.798706212352923e+07, -1.027267982590785e+07,
      0.000000000000000e-01, -3.466327490120249e+00, 3.602867454806401e+02,
      -1.223518159744950e+04, 2.015152850494290e+05, -1.909713800970267e+06,
      1.142278748458884e+07, -4.551909122318173e+07, 1.246206804545239e+08,
      -2.376275441309645e+08, 3.149824199011395e+08, -2.844660497989077e+08,
      1.668669076381706e+08, -5.729784575448167e+07, 8.737812306213077e+06,
      0.000000000000000e-01, 2.512080922932578e+00, -2.612119914965073e+02,
      8.878143206793750e+03, -1.464124202177336e+05, 1.389929291394544e+06,
      -8.332053211967151e+06, 3.329158201137647e+07, -9.143262460433713e+07,
      1.749809182259230e+08, -2.329047114119376e+08, 2.113183971320489e+08,
      -1.245975118891604e+08, 4.302553108480184e+07, -6.601663651220939e+06,
      0.000000000000000e-01, -1.000000000000000e+00, 1.040000000000000e+02,
      -3.536000000000000e+03, 5.834400000000000e+04, -5.542680000000000e+05,
      3.325608000000000e+06, -1.330243200000000e+07, 3.658168800000000e+07,
      -7.011490200000000e+07, 9.348653600000000e+07, -8.498776000000000e+07,
      5.022004000000000e+07, -1.738386000000000e+07, 2.674440000000000e+06 };
 
  static const double fem_coef_gausslob_16[289] =
    { 1.000000000000000e+00, -1.360000000000000e+02, 6.120000000000000e+03,
      -1.356600000000000e+05, 1.763580000000000e+06, -1.481407200000000e+07,
      8.535727200000000e+07, -3.505745100000000e+08, 1.051723530000000e+09,
      -2.337163400000000e+09, 3.866943080000000e+09, -4.745793780000000e+09,
      4.259045700000000e+09, -2.714556600000000e+09, 1.163381400000000e+09,
      -3.005401950000000e+08, 3.535767000000000e+07, 0.000000000000000e-01,
      1.839908474110340e+02, -1.132675577023645e+04, 2.828774748172428e+05,
      -3.903228397113182e+06, 3.393217989215092e+07, -1.997936813387011e+08,
      8.326183533203730e+08, -2.523654441487500e+09, 5.650496513849965e+09,
      -9.402288564814163e+09, 1.159004125992500e+10, -1.043753556167033e+10,
      6.671119113102151e+09, -2.865585732582308e+09, 7.416762022559170e+08,
      -8.739414676532781e+07, 0.000000000000000e-01, -7.365158560983221e+01,
      8.363752486870456e+03, -2.630512922725150e+05, 4.088268892503375e+06,
      -3.814296099037553e+07, 2.350888858038512e+08, -1.010915772862718e+09,
      3.133749897384450e+09, -7.134585101325640e+09, 1.202392366962845e+10,
      -1.496979454364896e+10, 1.358833701849443e+10, -8.740756338653663e+09,
      3.774397412355719e+09, -9.811765345365573e+08, 1.160408606498937e+08,
      0.000000000000000e-01, 4.208515410469884e+01, -5.266887544418802e+03,
      2.004067164817591e+05, -3.534528216742270e+06, 3.586506864599599e+07,
      -2.342572871648561e+08, 1.050195570375994e+09, -3.357626237614668e+09,
      7.826157982970905e+09, -1.343314504492630e+10, 1.696909064108448e+10,
      -1.558480944857237e+10, 1.012167818105129e+10, -4.405611470443590e+09,
      1.152926420144654e+09, -1.371250292488943e+08, 0.000000000000000e-01,
      -2.804154671449467e+01, 3.632134644860107e+03, -1.481030987143817e+05,
      2.845430205439061e+06, -3.103475950537538e+07, 2.145184209516461e+08,
      -1.004950951655481e+09, 3.325504487572721e+09, -7.965634686855131e+09,
      1.397533464514550e+10, -1.797134014690632e+10, 1.674913384367476e+10,
      -1.101142723314074e+10, 4.842306972881947e+09, -1.278281396603255e+09,
      1.531698732399162e+08, 0.000000000000000e-01, 2.036791285538472e+01,
      -2.681212491930484e+03, 1.129589658306899e+05, -2.270501508829428e+06,
      2.601887714173907e+07, -1.882644410377163e+08, 9.175357517098464e+08,
      -3.139134215049928e+09, 7.731773974122787e+09, -1.388518223176127e+10,
      1.820883321323428e+10, -1.725391802961138e+10, 1.150422262078233e+10,
      -5.120395806896533e+09, 1.365808881323657e+09, -1.651383905702324e+08,
      0.000000000000000e-01, -1.562975981216744e+01, 2.075857199589987e+03,
      -8.904128307330524e+04, 1.836683464233000e+06, -2.171339625537650e+07,
      1.623702309776589e+08, -8.168673425840398e+08, 2.877184244039678e+09,
      -7.272633187920932e+09, 1.336161532561319e+10, -1.787465369706492e+10,
      1.723415209556249e+10, -1.166677727073840e+10, 5.262202876034756e+09,
      -1.420107794637259e+09, 1.734782145645626e+08, 0.000000000000000e-01,
      1.245158740600359e+01, -1.662689739514946e+03, 7.210078018619273e+04,
      -1.511262926531350e+06, 1.823010400853032e+07, -1.394732137018548e+08,
      7.186625912313578e+08, -2.591802100670653e+09, 6.699969496687534e+09,
      -1.256822106464210e+10, 1.713562111741701e+10, -1.680796707246793e+10,
      1.155571599892056e+10, -5.285087289024682e+09, 1.444204616664480e+09,
      -1.784123720377501e+08, 0.000000000000000e-01, -1.018430458430458e+01,
      1.364696814296814e+03, -5.959855042735043e+04, 1.262405659052059e+06,
      -1.543602456068376e+07, 1.199989722604507e+08, -6.293065120124320e+08,
      2.311744565308469e+09, -6.087583637383061e+09, 1.162721665412277e+10,
      -1.612769282864336e+10, 1.607722369253147e+10, -1.122097126220979e+10,
      5.203928701314685e+09, -1.440373122685315e+09, 1.800466403356643e+08,
      0.000000000000000e-01, 8.484036081962663e+00, -1.139564172918365e+03,
      5.000628114583172e+04, -1.066865683702600e+06, 1.316848914076596e+07,
      -1.035421266853741e+08, 5.500823983897466e+08, -2.049399257475671e+09,
      5.477078740096760e+09, -1.061962761323504e+10, 1.495184835525608e+10,
      -1.512401891411349e+10, 1.070494963879464e+10, -5.031502683587495e+09,
      1.410393335939522e+09, -1.784123720377501e+08, 0.000000000000000e-01,
      -7.151250283691663e+00, 9.621468006776697e+02, -4.236328367402523e+04,
      9.083924059514159e+05, -1.128778312795541e+07, 8.948673396532559e+07,
      -4.799806642461592e+08, 1.807454740270899e+09, -4.886699456126079e+09,
      9.591077381959552e+09, -1.367409274689175e+10, 1.400781324267719e+10,
      -1.004054471299105e+10, 4.777971704223384e+09, -1.355543638395742e+09,
      1.734782145645626e+08, 0.000000000000000e-01, 6.060147075943005e+00,
      -8.163167553056976e+02, 3.602890134025652e+04, -7.753708269797396e+05,
      9.681484150831325e+06, -7.721340002337801e+07, 4.170906086685726e+08,
      -1.583343835764609e+09, 4.319156776154927e+09, -8.559301071696637e+09,
      1.232825943540154e+10, -1.276387222258454e+10, 9.248884856115279e+09,
      -4.449869455469566e+09, 1.276405367800062e+09, -1.651383905702324e+08,
      0.000000000000000e-01, -5.123522643555085e+00, 6.907394286218810e+02,
      -3.053901287681582e+04, 6.589382296622791e+05, -8.256408016568484e+06,
      6.613528180437879e+07, -3.591109349416936e+08, 1.371451685008549e+09,
      -3.766497172573159e+09, 7.519828793959344e+09, -1.091857986975193e+10,
      1.140164818726860e+10, -8.336452758217789e+09, 4.048470812623064e+09,
      -1.172436575235404e+09, 1.531698732399162e+08, 0.000000000000000e-01,
      4.271888353674923e+00, -5.762713061877649e+02, 2.550919052304026e+04,
      -5.514258520673562e+05, 6.926418073801134e+06, -5.565457178175802e+07,
      3.033329010654617e+08, -1.163492208450654e+09, 3.211252052317023e+09,
      -6.446888262886903e+09, 9.417864122967573e+09, -9.899668972442366e+09,
      7.289624669351123e+09, -3.566718678141100e+09, 1.041074047837655e+09,
      -1.371250292488943e+08, 0.000000000000000e-01, -3.434977405385288e+00,
      4.635617485626016e+02, -2.053688028669370e+04, 4.444943513906060e+05,
      -5.592632684586483e+06, 4.503253959356561e+07, -2.460675245081598e+08,
      9.466718497394927e+08, -2.621823641133609e+09, 5.284002694315866e+09,
      -7.752423037115038e+09, 8.187712309040204e+09, -6.060153271928369e+09,
      2.981652672294607e+09, -8.754772358617425e+08, 1.160408606498937e+08,
      0.000000000000000e-01, 2.505373764729175e+00, -3.381913429670035e+02,
      1.499009100007397e+04, -3.246847962658711e+05, 4.089321087106837e+06,
      -3.296978262323839e+07, 1.804331430493320e+08, -6.954301078105766e+08,
      1.930060872117711e+09, -3.899125665782255e+09, 5.735918309736332e+09,
      -6.075963842786729e+09, 4.511802094762441e+09, -2.227740310582889e+09,
      6.566301459893279e+08, -8.739414676532781e+07, 0.000000000000000e-01,
      -1.000000000000000e+00, 1.350000000000000e+02, -5.985000000000000e+03,
      1.296750000000000e+05, -1.633905000000000e+06, 1.318016700000000e+07,
      -7.217710500000000e+07, 2.783974050000000e+08, -7.733261250000000e+08,
      1.563837275000000e+09, -2.303105805000000e+09, 2.442687975000000e+09,
      -1.816357725000000e+09, 8.981988750000000e+08, -2.651825250000000e+08,
      3.535767000000000e+07 };
 
  static const double fem_coef_gausslob_24[625] =
    { 1.000000000000000e+00, -3.000000000000000e+02, 2.990000000000000e+04,
      -1.480050000000000e+06, 4.351347000000000e+07, -8.412604200000000e+08,
      1.141710570000000e+10, -1.137633032250000e+11, 8.595449577000000e+11,
      -5.042663751840000e+12, 2.337962284944000e+13, -8.678799391080000e+13,
      2.603639817324000e+14, -6.351736697208000e+14, 1.264298066396640e+15,
      -2.054484357894540e+15, 2.719170473683950e+15, -2.914666390092600e+15,
      2.505590405518200e+15, -1.701164012167620e+15, 8.910859111354200e+14,
      -3.471763290138000e+14, 9.468445336740000e+13, -1.612380184155000e+13,
      1.289904147324000e+12, 0.000000000000000e-01, 4.058641271792565e+02,
      -5.527892747506228e+04, 3.080680865091051e+06, -9.608544697986574e+07,
      1.921815551511292e+09, -2.664514371195282e+10, 2.693344195607988e+11,
      -2.055620699316044e+12, 1.214897536784061e+13, -5.664114215003022e+13,
      2.111637034506248e+14, -6.356401824720572e+14, 1.554903769748071e+15,
      -3.101863582054997e+15, 5.049759901096915e+15, -6.693732300831693e+15,
      7.184248666000200e+15, -6.182729556103940e+15, 4.201716500300169e+15,
      -2.202700032171432e+15, 8.588067464630930e+14, -2.343664562972425e+14,
      3.993223187351215e+13, -3.196139551478179e+12, 0.000000000000000e-01,
      -1.624756704002590e+02, 4.076566744849383e+04, -2.856559180759908e+06,
      1.002242332396244e+08, -2.149192199655287e+09, 3.116591524589708e+10,
      -3.248555446646729e+11, 2.534404897909630e+12, -1.522400128496926e+13,
      7.186059820369571e+13, -2.704952845330370e+14, 8.204876425573068e+14,
      -2.019503996546278e+15, 4.049106649016715e+15, -6.619542027745681e+15,
      8.805466617020362e+15, -9.478911036302426e+15, 8.178282706527972e+15,
      -5.570062256804815e+15, 2.925593222099543e+15, -1.142548313369483e+15,
      3.122535907953383e+14, -5.327144417235540e+13, 4.268671053543734e+12,
      0.000000000000000e-01, 9.285790471348058e+01, -2.567292232023452e+04,
      2.172505008713478e+06, -8.632693629032278e+07, 2.009759260312978e+09,
      -3.083881003217180e+10, 3.346965047874316e+11, -2.690206053619648e+12,
      1.652940573703274e+13, -7.940281485029880e+13, 3.030598877175188e+14,
      -9.295754861643471e+14, 2.308921995608790e+15, -4.664329550291535e+15,
      7.673380080847903e+15, -1.026158278560169e+16, 1.109639715400819e+16,
      -9.611028022799504e+15, 6.567868332535721e+15, -3.459765425094874e+15,
      1.354622808613475e+15, -3.710479852977045e+14, 6.342841599973944e+13,
      -5.091588188794019e+12, 0.000000000000000e-01, -6.190263201810896e+01,
      1.771295817306943e+04, -1.605426884400296e+06, 6.937138005497782e+07,
      -1.732266817595608e+09, 2.807090718452186e+10, -3.177491729312126e+11,
      2.638958096874092e+12, -1.663806369766778e+13, 8.158796865135113e+13,
      -3.166342096198069e+14, 9.845663378570931e+14, -2.473337869447428e+15,
      5.044014753850598e+15, -8.364663482728635e+15, 1.126254226203764e+16,
      -1.225027170970453e+16, 1.066427200479037e+16, -7.319785257851210e+15,
      3.870748078365627e+15, -1.520689469004945e+15, 4.177883147633562e+14,
      -7.160927970988626e+13, 5.762006100161049e+12, 0.000000000000000e-01,
      4.501042175430940e+01, -1.308958048352837e+04, 1.225547369430789e+06,
      -5.535760993480313e+07, 1.449945861634959e+09, -2.454369660090688e+10,
      2.883865544213490e+11, -2.470900824125910e+12, 1.598637745987517e+13,
      -8.009303566237079e+13, 3.164491556576764e+14, -9.988976501079246e+14,
      2.541436559580714e+15, -5.239264041161704e+15, 8.769361705836964e+15,
      -1.190220432935039e+16, 1.303612471279278e+16, -1.141726255498977e+16,
      7.878336148801056e+15, -4.185657107040718e+15, 1.651238357905402e+15,
      -4.553312625304459e+14, 7.830179054235796e+13, -6.319165567954503e+12,
      0.000000000000000e-01, -3.460823180120079e+01, 1.015469610544880e+04,
      -9.679531870145017e+05, 4.485134755229398e+07, -1.210735945225935e+09,
      2.114609948086286e+10, -2.559531795223062e+11, 2.252595670603948e+12,
      -1.492191456512238e+13, 7.630937241452699e+13, -3.068986895462340e+14,
      9.837309848042438e+14, -2.536329824251836e+15, 5.289429958022946e+15,
      -8.942835505684019e+15, 1.224496508118832e+16, -1.351567177210176e+16,
      1.191830727462528e+16, -8.273935552638878e+15, 4.419525275428065e+15,
      -1.751884239411906e+15, 4.851652818736986e+14, -8.375521716296401e+13,
      6.782865257514837e+12, 0.000000000000000e-01, 2.766582877253528e+01,
      -8.161941725332768e+03, 7.865526415735044e+05, -3.702889412503837e+07,
      1.019390720188427e+09, -1.819645578676150e+10, 2.252249000861092e+11,
      -2.025483052842698e+12, 1.369084263890056e+13, -7.131369560391163e+13,
      2.915943262851778e+14, -9.485944639479267e+14, 2.478119277999363e+15,
      -5.228788122037598e+15, 8.932615077176731e+15, -1.234455468758280e+16,
      1.373835480205531e+16, -1.220422101884528e+16, 8.528504524807064e+15,
      -4.582579019928659e+15, 1.826238758482371e+15, -5.082003015125067e+14,
      8.811547249928783e+13, -7.164301017225998e+12, 0.000000000000000e-01,
      -2.275670591599455e+01, 6.737590226239552e+03, -6.539504206517175e+05,
      3.111139106156313e+07, -8.679722967903971e+08, 1.573365420642779e+10,
      -1.979909637705803e+11, 1.810880611538906e+12, -1.244463025448231e+13,
      6.585374856883214e+13, -2.732735801573156e+14, 9.011914571652850e+14,
      -2.383831456114384e+15, 5.087291841992724e+15, -8.780953301699755e+15,
      1.224889816733013e+16, -1.374781660040675e+16, 1.230672077620068e+16,
      -8.660226214361099e+15, 4.682883135541219e+15, -1.876981572450838e+15,
      5.250670186065034e+14, -9.147680701891190e+13, 7.470231264323625e+12,
      0.000000000000000e-01, 1.912965144613660e+01, -5.677631395079384e+03,
      5.537935689545739e+05, -2.653927821530689e+07, 7.474036306614312e+08,
      -1.369940651896189e+10, 1.745319507247928e+11, -1.617301631376405e+12,
      1.126327452432594e+13, -6.039298017983792e+13, 2.538313181524926e+14,
      -8.473116290006231e+14, 2.267097817730376e+15, -4.890112526330036e+15,
      8.524655979372967e+15, -1.200076621013731e+16, 1.358349523272153e+16,
      -1.225446423734989e+16, 8.685297402148473e+15, -4.727405822382261e+15,
      1.906317065479628e+15, -5.362492238636365e+14, 9.390516767804315e+13,
      -7.704880889559378e+12, 0.000000000000000e-01, -1.635497697368700e+01,
      4.862656953497025e+03, -4.759804636482565e+05, 2.293041632316653e+07,
      -6.502015538842764e+08, 1.201606379408697e+10, -1.545199230484324e+11,
      1.446437458063616e+12, -1.018096106218994e+13, 5.518468570532575e+13,
      -2.344620385089091e+14, 7.909885622725649e+14, -2.138165417503776e+15,
      4.657339854616682e+15, -8.194527986057532e+15, 1.163730420250179e+16,
      -1.328057957790062e+16, 1.207345046837354e+16, -8.618482475976474e+15,
      4.722435616707614e+15, -1.916176031917556e+15, 5.421466879431477e+14,
      -9.544980641239802e+13, 7.870911362255844e+12, 0.000000000000000e-01,
      1.417066207624127e+01, -4.218698704765500e+03, 4.140273580736256e+05,
      -2.002373113338682e+07, 5.706909517514257e+08, -1.061235738178459e+10,
      1.374488760853391e+11, -1.296867134188179e+12, 9.206001697332345e+12,
      -5.034424091354864e+13, 2.158419782218546e+14, -7.348173021384713e+14,
      2.004252206567683e+15, -4.404149060382797e+15, 7.815178748758451e+15,
      -1.118956431220418e+16, 1.286957131473114e+16, -1.178684055128811e+16,
      8.473168167228462e+15, -4.673705358375533e+15, 1.908297467260950e+15,
      -5.431049066701974e+14, 9.614927445703371e+13, -7.969947411626695e+12,
      0.000000000000000e-01, -1.240846755882427e+01, 3.697723332529632e+03,
      -3.636177333437864e+05, 1.763791694375029e+07, -5.046596469793725e+08,
      9.429433336134134e+09, -1.228099190218494e+11, 1.166008419408402e+12,
      -8.333618884154624e+12, 4.590449180645647e+13, -1.982963080519099e+14,
      6.803133908337801e+14, -1.870091239145240e+15, 4.141349396659428e+15,
      -7.405302748248103e+15, 1.068239700855010e+16, -1.237594459251991e+16,
      1.141465416121965e+16, -8.261228940134621e+15, 4.586380576952005e+15,
      -1.884249466545215e+15, 5.394272826690079e+14, -9.603440259637641e+13,
      8.002866883031367e+12, 0.000000000000000e-01, 1.095558179466470e+01,
      -3.267249008495488e+03, 3.217786804294041e+05, -1.564425754031200e+07,
      4.489762785029423e+08, -8.420409805207488e+09, 1.101506623685581e+11,
      -1.051033145830932e+12, 7.553210200363392e+12, -4.185258869687206e+13,
      1.819278279112502e+14, -6.282336154227903e+14, 1.738507104235632e+15,
      -3.876120340749998e+15, 6.978305450391071e+15, -1.013471839778258e+16,
      1.182005159615111e+16, -1.097352991130819e+16, 7.992846614334175e+15,
      -4.464988119249731e+15, 1.845417602986293e+15, -5.313770797673897e+14,
      9.512946342200696e+13, -7.969947411626695e+12, 0.000000000000000e-01,
      -9.733404845062563e+00, 2.904495367336991e+03, -2.863957452026712e+05,
      1.394908621565543e+07, -4.012835563442479e+08, 7.548227155070059e+09,
      -9.908687724892865e+10, 9.492474279583278e+11, -6.852122094227041e+12,
      3.815223404048325e+13, -1.667054040084672e+14, 5.788252023371522e+14,
      -1.610924210936436e+15, 3.612762299121459e+15, -6.543084510709907e+15,
      9.560030643675167e+15, -1.121725598566102e+16, 1.047660019645038e+16,
      -7.676343347474552e+15, 4.313320840279778e+15, -1.792974677700824e+15,
      5.191726764406062e+14, -9.345206628174223e+13, 7.870911362255844e+12,
      0.000000000000000e-01, 8.685152146887849e+00, -2.592917301003156e+03,
      2.559159080755008e+05, -1.248235381982623e+07, 3.597715754514247e+08,
      -6.783361410773893e+09, 8.929618961811727e+10, -8.582135820103382e+11,
      6.217423158689400e+12, -3.475607425638306e+13, 1.525197220141970e+14,
      -5.320009405626489e+14, 1.487763351500956e+15, -3.353349463718517e+15,
      6.104800041187937e+15, -8.967037300250809e+15, 1.057820092203265e+16,
      -9.933456336414407e+15, 7.318037585534789e+15, -4.134330534453634e+15,
      1.727837357443638e+15, -5.029774927870322e+14, 9.101197367138191e+13,
      -7.704880889559378e+12, 0.000000000000000e-01, -7.768227503689160e+00,
      2.320048262018063e+03, -2.291579825895192e+05, 1.118998178170061e+07,
      -3.230127412726743e+08, 6.101825984880789e+09, -8.050592964465223e+10,
      7.757517929934399e+11, -5.636578411502579e+12, 3.161187883403309e+13,
      -1.392153217132889e+14, 4.874510273126649e+14, -1.368719368735070e+15,
      3.098228256202772e+15, -5.665515780120550e+15, 8.360206053329256e+15,
      -9.909089467205730e+15, 9.350136076665120e+15, -6.922098442148973e+15,
      3.930003596385763e+15, -1.650608740098542e+15, 4.828842861248500e+14,
      -8.780874332485509e+13, 7.470231264323625e+12, 0.000000000000000e-01,
      6.949251795654158e+00, -2.076080736426886e+03, 2.051850668187582e+05,
      -1.002851556374879e+07, 2.898385274463058e+08, -5.483488755332865e+09,
      7.247948130849744e+10, -6.998845716368053e+11, 5.097508994937879e+12,
      -2.866481138828522e+13, 1.266058930533223e+14, -4.447041797385597e+14,
      1.252927498139101e+15, -2.846337193255717e+15, 5.224630116918579e+15,
      -7.740148279482860e+15, 9.211839907560781e+15, -8.729031202403334e+15,
      6.490342376708594e+15, -3.701195553992629e+15, 1.561498591338377e+15,
      -4.588915147832623e+14, 8.382775191413614e+13, -7.164301017225998e+12,
      0.000000000000000e-01, -6.200545901231066e+00, 1.852852310460863e+03,
      -1.832115058838774e+05, 8.961081566680860e+06, -2.592406841374813e+08,
      4.910586592056635e+09, -6.500190149790611e+10, 6.287466876199138e+11,
      -4.588252565931313e+12, 2.585696625207004e+13, -1.144768233740891e+14,
      4.031460020145847e+14, -1.139023606360627e+15, 2.595327972551887e+15,
      -4.779021130665880e+15, 7.103675353349177e+15, -8.483943776806492e+15,
      8.068562477979859e+15, -6.021870692282251e+15, 3.447372991345801e+15,
      -1.460201300789598e+15, 4.308660981996213e+14, -7.903354901739207e+13,
      6.782865257514837e+12, 0.000000000000000e-01, 5.497265260133587e+00,
      -1.643010914375697e+03, 1.625245542407611e+05, -7.953853255844527e+06,
      2.302798044521235e+08, -4.366227075820252e+09, 5.786337032883471e+10,
      -5.604566478207931e+11, 4.096239643506256e+12, -2.312433390636899e+13,
      1.025754064938261e+14, -3.619933583977279e+14, 1.025085118853682e+15,
      -2.341436140046740e+15, 4.322778692729082e+15, -6.443312152718818e+15,
      7.717747123310279e+15, -7.362354286134902e+15, 5.512353176524103e+15,
      -3.166155510128892e+15, 1.345687520087759e+15, -3.984797768116558e+14,
      7.335818308855012e+13, -6.319165567954503e+12, 0.000000000000000e-01,
      -4.814420396783881e+00, 1.439137261510244e+03, -1.424001050225189e+05,
      6.972107765750786e+06, -2.019777804539574e+08, 3.832494908082456e+09,
      -5.083618287186299e+10, 4.929144471099037e+11, -3.606960360420349e+12,
      2.039001482865880e+13, -9.058350359554217e+13, 3.202053356214614e+14,
      -9.083926524953095e+14, 2.078951013598256e+15, -3.846222877659622e+15,
      5.745792065253601e+15, -6.898564020656367e+15, 6.597335811773795e+15,
      -4.952528004193797e+15, 2.852412393699805e+15, -1.215806035913631e+15,
      3.610885650804218e+14, -6.667886669397892e+13, 5.762006100161049e+12,
      0.000000000000000e-01, 4.122502768575858e+00, -1.232445305915745e+03,
      1.219756720552666e+05, -5.974119340666288e+06, 1.731450552812862e+08,
      -3.287266438655767e+09, 4.363384253154989e+10, -4.234185297825935e+11,
      3.101259925467429e+12, -1.754945237689139e+13, 7.805398522969187e+13,
      -2.762644893141324e+14, 7.848217578377621e+14, -1.798840649081661e+15,
      3.333370625337169e+15, -4.988259107765304e+15, 6.000070847701146e+15,
      -5.749271353120764e+15, 4.324788417805135e+15, -2.496262406568326e+15,
      1.066418114121039e+15, -3.174727574108465e+14, 5.876970053131701e+13,
      -5.091588188794019e+12, 0.000000000000000e-01, -3.378098091794216e+00,
      1.009981094027902e+03, -9.997415292076045e+04, 4.897701324351262e+06,
      -1.419932385393241e+08, 2.696914741487987e+09, -3.581511558046563e+10,
      3.477438352489987e+11, -2.548653261914288e+12, 1.443296944374584e+13,
      -6.424560812216257e+13, 2.275969917931459e+14, -6.472060159591627e+14,
      1.485016620793322e+15, -2.755030677045359e+15, 4.127938042773006e+15,
      -4.971860898078912e+15, 4.770796079822418e+15, -3.594142516031414e+15,
      2.077829100777859e+15, -8.891455208945626e+14, 2.651635856092349e+14,
      -4.917666111269423e+13, 4.268671053543734e+12, 0.000000000000000e-01,
      2.493031698759675e+00, -7.454011644125376e+02, 7.379166798110049e+04,
      -3.615566630393531e+06, 1.048426846839652e+08, -1.991802210178859e+09,
      2.645916950749144e+10, -2.569938254028299e+11, 1.884300408926143e+12,
      -1.067564580288449e+13, 4.754491961430585e+13, -1.685282542848975e+14,
      4.795322158961934e+14, -1.101030336950954e+15, 2.044141709763629e+15,
      -3.065196721720791e+15, 3.694953407319221e+15, -3.548705940334544e+15,
      2.676012339710790e+15, -1.548605987126603e+15, 6.633870802695019e+14,
      -1.980596394144405e+14, 3.677511736196415e+13, -3.196139551478179e+12,
      0.000000000000000e-01, -1.000000000000000e+00, 2.990000000000000e+02,
      -2.960100000000000e+04, 1.450449000000000e+06, -4.206302100000000e+07,
      7.991973990000000e+08, -1.061790830100000e+10, 1.031453949240000e+11,
      -7.563995627760000e+11, 4.286264189064000e+12, -1.909335866037600e+13,
      6.769463525042400e+13, -1.926693464819760e+14, 4.425043232388240e+14,
      -8.217937431578160e+14, 1.232690614736724e+15, -1.486479858947226e+15,
      1.428186531145374e+15, -1.077403874372826e+15, 6.237601377947940e+14,
      -2.673257733406260e+14, 7.985055567317400e+13, -1.483389769422600e+13,
      1.289904147324000e+12 };
 
  static const double fem_coef_gausslob_32[1089] =
    { 1.000000000000000e+00, -5.280000000000000e+02, 9.275200000000000e+04,
      -8.115800000000000e+06, 4.236447600000000e+08, -1.462986571200000e+10,
      3.573867195360000e+11, -6.471252385884000e+12, 8.987850535950000e+13,
      -9.826716585972000e+14, 8.629643838226320e+15, -6.184578084062196e+16,
      3.663173172867608e+17, -1.811459261308158e+18, 7.539120925634905e+18,
      -2.657540126286304e+19, 7.972620378858912e+19, -2.042658293145551e+20,
      4.479513800757788e+20, -8.416770667739633e+20, 1.354699278902855e+21,
      -1.864910695632502e+21, 2.189242990525111e+21, -2.181310950704368e+21,
      1.832301198591669e+21, -1.285429763935079e+21, 7.434251911077520e+20,
      -3.481117958361696e+20, 1.286127324517868e+20, -3.607069737728273e+19,
      7.214139475456547e+18, -9.163120704712953e+17, 5.553406487704820e+16,
      0.000000000000000e-01, 7.143214682034540e+02, -1.714133648848133e+05,
      1.688198752795928e+07, -9.347155797218437e+08, 3.338933321339879e+10,
      -8.331872976904955e+11, 1.530331866800279e+13, -2.146891160964000e+14,
      2.364531150945150e+15,  -2.087974965454055e+16, 1.502766434956253e+17,
      -8.930913114593062e+17, 4.428285547618542e+18, -1.847053591535166e+19,
      6.522650112096549e+19, -1.959752299775458e+20, 5.027465769775919e+20,
      -1.103711061615762e+21, 2.075747280986809e+21, -3.343657099618051e+21,
      4.606186697816854e+21, -5.410587449851726e+21, 5.393904740095235e+21,
      -4.533054367123241e+21, 3.181470251931304e+21, -1.840698151594781e+21,
      8.622096114618880e+20, -3.186487205259235e+20, 8.939310241964093e+19,
      -1.788315002192243e+19, 2.271972224754086e+18, -1.377242653770284e+17,
      0.000000000000000e-01, -2.859599139140576e+02, 1.263498115660723e+05,
      -1.563762114667972e+07, 9.735262028286150e+08, -3.727077001705718e+10,
      9.724728419064618e+11, -1.841437932011095e+13, 2.640184847320257e+14,
      -2.955001754223333e+15, 2.641501188215526e+16, -1.919333187312935e+17,
      1.149300695265413e+18, -5.733477927220639e+18, 2.403400433666191e+19,
      -8.522440623029794e+19, 2.569470483161848e+20, -6.610937005531425e+20,
      1.454972285683492e+21, -2.742254919882459e+21, 4.425523119288710e+21,
      -6.106476088763322e+21, 7.183124069502222e+21, -7.170000472640595e+21,
      6.032425263760434e+21, -4.238001893187230e+21, 2.454158919196046e+21,
      -1.150483843370953e+21, 4.254938488755735e+20, -1.194452209653931e+20,
      2.390927098916194e+19, -3.039204212011252e+18, 1.843238572103207e+17,
      0.000000000000000e-01, 1.634368826450262e+02, -7.956978254181237e+04,
      1.188506225092293e+07, -8.373897346830345e+08, 3.478333876388285e+10,
      -9.598393772867426e+11, 1.891593013667364e+13, -2.793128536050818e+14,
      3.196655214719607e+15, -2.907291848273606e+16, 2.141468780998689e+17,
      -1.296440193900301e+18, 6.525499955570070e+18, -2.755634128784674e+19,
      9.831738858873329e+19, -2.979627701541287e+20, 7.700118175018154e+20,
      -1.701110443974214e+21, 3.216663602425906e+21, -5.205922417479511e+21,
      7.201199888459469e+21, -8.489441777094031e+21, 8.490373918346112e+21,
      -7.155631502801035e+21, 5.034836798969603e+21, -2.919618180375494e+21,
      1.370386674687711e+21, -5.073912120857659e+20, 1.425804015450344e+20,
      -2.856660788223791e+19, 3.634277715764878e+18, -2.205841595819251e+17,
      0.000000000000000e-01, -1.089624682152292e+02, 5.490286557349962e+04,
      -8.781653856365252e+06, 6.724119974793967e+08, -2.993574836604601e+10,
      8.717420300323828e+11, -1.790617757376325e+13, 2.730388108651176e+14,
      -3.204823752640508e+15, 2.974036495407882e+16, -2.226577421064578e+17,
      1.366028725355318e+18, -6.951897867287077e+18, 2.962837694696566e+19,
      -1.065339992825695e+20, 3.250038142137898e+20, -8.446629720680499e+20,
      1.875178351120271e+21, -3.560917379168056e+21, 5.784532874816462e+21,
      -8.027770029015426e+21, 9.491255762425382e+21, -9.516676697643140e+21,
      8.038962462344214e+21, -5.667965587919967e+21, 3.292818343972966e+21,
      -1.548126823651629e+21, 5.740631512413311e+20, -1.615361972833789e+20,
      3.240477969405117e+19, -4.127262310667621e+18, 2.507669351857205e+17,
      0.000000000000000e-01, 7.924220268285228e+01, -4.057928795465258e+04,
      6.704332297502558e+06, -5.364604956918865e+08, 2.503646198522573e+10,
      -7.610195551213778e+11, 1.621371476972958e+13, -2.548664517300070e+14,
      3.067722734616486e+15, -2.906734241270257e+16, 2.214249975700081e+17,
      -1.378338812356282e+18, 7.101001568235638e+18, -3.058038376939241e+19,
      1.109399051188740e+20, -3.410471802537065e+20, 8.922579803599067e+20,
      -1.992320553154499e+21, 3.802564715982721e+21, -6.204660474711801e+21,
      8.644825045503473e+21, -1.025665248235139e+22, 1.031630159337718e+22,
      -8.738843624083011e+21, 6.176944575139345e+21, -3.596664164242111e+21,
      1.694457654020145e+21, -6.294971162057861e+20, 1.774358410868178e+20,
      -3.564953335849152e+19, 4.546981308698057e+18, -2.766285114923478e+17,
      0.000000000000000e-01, -6.094810716099949e+01, 3.149084615481936e+04,
      -5.296674502372112e+06, 4.346997745238327e+08, -2.090081537881765e+10,
      6.551264846374850e+11, -1.436792716105805e+13, 2.318076448212311e+14,
      -2.854539793592080e+15, 2.758693070051896e+16, -2.137570364709653e+17,
      1.350278421513405e+18, -7.045142142783022e+18, 3.067459655658656e+19,
      -1.123483766643411e+20, 3.482651662900422e+20, -9.178174833360337e+20,
      2.062604456625362e+21, -3.959135109240898e+21, 6.492786887393485e+21,
      -9.086987198956428e+21, 1.082464243206250e+22, -1.092690369776145e+22,
      9.286162259274268e+21, -6.583075803335110e+21, 3.843334757170004e+21,
      -1.815042549975247e+21, 6.757786079688690e+20, -1.908640091400089e+20,
      3.841802724253007e+19, -4.908370532055551e+18, 2.990786503802649e+17,
      0.000000000000000e-01, 4.874762271398016e+01, -2.532459924379460e+04,
      4.306289111489640e+06, -3.590409832532830e+08, 1.760136771205027e+10,
      -5.636351004936871e+11, 1.263327395525212e+13, -2.081295681322519e+14,
      2.613155513828582e+15, -2.570230208401814e+16, 2.023153786109967e+17,
      -1.296022489490302e+18, 6.846470251956364e+18, -3.013872322115386e+19,
      1.114644422978128e+20, -3.485183237291328e+20, 9.255531281005032e+20,
      -2.094244828319674e+21, 4.044473148145946e+21, -6.669094801072411e+21,
      9.379689895649970e+21, -1.122286017421572e+22, 1.137425276024759e+22,
      -9.701397530693913e+21, 6.900090712824613e+21, -4.040491989417283e+21,
      1.913372195656766e+21, -7.141717651945982e+20, 2.021704584417841e+20,
      -4.077964478229480e+19, 5.220210907153068e+18, -3.186495777814819e+17,
      0.000000000000000e-01, -4.013085745911482e+01, 2.092265556924549e+04,
      -3.583307397906160e+06, 3.019036845311046e+08, -1.499682567334736e+10,
      4.875419883971985e+11, -1.110534210962378e+13, 1.859662144990967e+14,
      -2.372232833495577e+15, 2.368570555021916e+16, -1.890606460337203e+17,
      1.226710982466906e+18, -6.556236864370271e+18, 2.916718348444145e+19,
      -1.089043462231537e+20, 3.434548724554294e+20, -9.192120758302752e+20,
      2.094521354663323e+21, -4.070706954543361e+21, 6.750946202167051e+21,
      -9.544297656724117e+21, 1.147387384551836e+22, -1.167874642179135e+22,
      1.000023491239336e+22, -7.138163819152158e+21, 4.193633752898508e+21,
      -1.991877420085289e+21, 7.455334138647318e+20, -2.115867266165027e+20,
      4.277941431548355e+19, -5.488110031116298e+18, 3.356769581362759e+17,
      0.000000000000000e-01, 3.377658303742900e+01, -1.765321539233928e+04,
      3.038340440514934e+06, -2.578584606807322e+08, 1.292884612296167e+10,
      -4.249332478575872e+11, 9.796453340445514e+12, -1.661322016030880e+14,
      2.146412285165976e+15, -2.170063067679739e+16, 1.753071524663780e+17,
      -1.150445213494342e+18, 6.214123504683469e+18, -2.791805131931260e+19,
      1.051885109311114e+20, -3.345072958360405e+20, 9.021185531229185e+20,
      -2.069976464960887e+21, 4.048800325514122e+21, -6.754020775002269e+21,
      9.599959543792231e+21, -1.159763408249418e+22, 1.185806715330365e+22,
      -1.019594067425241e+22, 7.305655578610319e+21, -4.307135697275206e+21,
      2.052428738060058e+21, -7.705007825196310e+20, 2.192794049399991e+20,
      -4.444874922276286e+19, 5.715873383290923e+18, -3.503832522669480e+17,
      0.000000000000000e-01, -2.892964123166297e+01, 1.514678464159201e+04,
      -2.616230196259168e+06, 2.232056370284014e+08, -1.126780271399590e+10,
      3.733563821201946e+11, -8.686293054983950e+12, 1.487584698214981e+14,
      -1.941627928264148e+15, 1.983312707085937e+16, -1.618550810194144e+17,
      1.072675099399285e+18, -5.848903282733471e+18, 2.651290247911275e+19,
      -1.007365700808825e+20, 3.228755039109242e+20, -8.771430385090572e+20,
      2.026394505724555e+21, -3.988616003559215e+21, 6.692583651440911e+21,
      -9.564190176741721e+21, 1.161237649908247e+22, -1.192826855978758e+22,
      1.030040441263462e+22, -7.409917359938946e+21, 4.384750316004612e+21,
      -2.096582122064279e+21, 7.895856427194460e+20, -2.253772189916786e+20,
      4.581095810605600e+19, -5.906211978906061e+18, 3.629208802827826e+17,
      0.000000000000000e-01, 2.513021588285159e+01, -1.317482888832123e+04,
      2.281636380391524e+06, -1.954241069277852e+08, 9.915879553852877e+09,
      -3.305907224268159e+11, 7.745610540950591e+12, -1.336744677758923e+14,
      1.759053047347007e+15, -1.812022604295268e+16, 1.491398073535739e+17,
      -9.967823580408796e+17, 5.480122870501913e+18, -2.504020332812559e+19,
      9.587106380784488e+19, -3.095239755572209e+20, 8.466795723536765e+20,
      -1.968748626171725e+21, 3.898845137794306e+21, -6.579450581556874e+21,
      9.452948974562387e+21, -1.153486704390072e+22, 1.190416296509360e+22,
      -1.032457212621363e+22, 7.457660044656362e+21, -4.429851231444799e+21,
      2.125705105699995e+21, -8.032242691112874e+20, 2.299857495593084e+20,
      -4.688429048988357e+19, 6.061137852550284e+18, -3.733965028540281e+17,
      0.000000000000000e-01, -2.208386882313948e+01,  1.158936145601140e+04,
      -2.011104324155879e+06, 1.727696977814699e+08, -8.800873732349352e+09,
      2.948204525687799e+11, -6.945679729271358e+12, 1.206045724118162e+14,
      -1.597549338531331e+15, 1.657073913780072e+16, -1.373597907391489e+17,
      9.246697496818755e+17, -5.120171153308929e+18, 2.356084529994870e+19,
      -9.082843603687474e+19, 2.951965214544448e+20, -8.126535950009281e+20,
      1.901181843699431e+21, -3.786945462658442e+21, 6.425892713139728e+21,
      -9.280555773871743e+21, 1.138038464790677e+22, -1.179939706531010e+22,
      1.027858871223671e+22, -7.455106011504378e+21, 4.445547883990988e+21,
      -2.141041178542259e+21, 8.118044630932303e+20, -2.331957228005324e+20,
      4.768370096691469e+19, -6.182197530650295e+18, 3.818855272205113e+17,
      0.000000000000000e-01, 1.959414243804128e+01, -1.029081802559516e+04,
      1.788568172854703e+06, -1.540118150914350e+08, 7.869521604579143e+09,
      -2.646147373681910e+11, 6.261419532268355e+12, -1.092584911616051e+14,
      1.455025808801897e+15, -1.517863644621452e+16, 1.265704706962783e+17,
      -8.572524657708660e+17, 4.776247526570586e+18, -2.211426166036566e+19,
      8.577380339701724e+19, -2.804435586124510e+20, 7.765584710933440e+20,
      -1.827036131219099e+21, 3.659136530638906e+21, -6.241580538370836e+21,
      9.059595877078572e+21, -1.116262879266932e+22, 1.162640401193756e+22,
      -1.017180635719499e+22, 7.408032173818105e+21, -4.434731590925861e+21,
      2.143740476291448e+21, -8.156798285468939e+20, 2.350876961959084e+20,
      -4.822189338581046e+19, 6.270614615621034e+18, -3.884411477495515e+17,
      0.000000000000000e-01, -1.752536768073206e+01, 9.210004044516183e+03,
      -1.602710362254713e+06, 1.382643168640165e+08, -7.082209191374808e+09,
      2.388593247336584e+11, -5.671955040601066e+12, 9.936819687224059e+13,
      -1.329133616843959e+15, 1.393095353694366e+16, -1.167468019716584e+17,
      7.948230830253194e+17, -4.451987054929648e+18, 2.072406047410437e+19,
      -8.081630829587058e+19, 2.656550395781670e+20, -7.395103797274482e+20,
      1.748920815324121e+21, -3.520456085792643e+21, 6.034597075619834e+21,
      -8.800877008017673e+21, 1.089363982096666e+22, -1.139632655550776e+22,
      1.001274028996965e+22, -7.321760958016250e+21, 4.400085249649023e+21,
      -2.134871514399776e+21, 8.151766211905410e+20, -2.357346041272461e+20,
      4.850993390791870e+19, -6.327377892472496e+18, 3.931001229299567e+17,
      0.000000000000000e-01, 1.578106132320405e+01, -8.297465341460668e+03,
      1.445356637015484e+06, -1.248763055461196e+08, 6.409121288934440e+09,
      -2.166867323836831e+11, 5.160255428324136e+12, -9.069980923376013e+13,
      1.217593153323901e+15, -1.281217702947176e+16, 1.078222149705347e+17,
      -7.373025946691738e+17, 4.148685852384210e+18, -1.940267191328230e+19,
      7.602301795140850e+19, -2.510934651183994e+20, 7.023105241061891e+20,
      -1.668804442314510e+21, 3.374862747204858e+21, -5.811516441898685e+21,
      8.513453677248631e+21, -1.058376696971405e+22, 1.111895644585412e+22,
      -9.809014479986894e+21, 7.201130357372178e+21, -4.344074707924148e+21,
      2.115422234515678e+21, -8.105961395642809e+20, 2.352029716821803e+20,
      -4.855759087275033e+19, 6.353294711027987e+18, -3.958864781209368e+17,
      0.000000000000000e-01, -1.429082487951404e+01, 7.516973886624383e+03,
      -1.310468641451437e+06, 1.133605392742571e+08, -5.827511997735239e+09,
      1.974178249055285e+11, -4.712515373341917e+12, 8.305450385112180e+13,
      -1.118328972822835e+15, 1.180653064142852e+16, -9.971166758181266e+16,
      6.844038883665075e+17, -3.866168854945464e+18, 1.815490936983781e+19,
      -7.143043815405257e+19, 2.369235292422868e+20, -6.655061581666021e+20,
      1.588114879269808e+21, -3.225364968659972e+21, 5.577529629049808e+21,
      -8.204710901105034e+21, 1.024169036500672e+22, -1.080270746628453e+22,
      9.567317871713345e+21, -7.050459812170524e+21, 4.268932026970228e+21,
      -2.086295163199683e+21, 8.022143863884770e+20, -2.335532639673231e+20,
      4.837349156604751e+19, -6.349020768043736e+18, 3.968137980027335e+17,
      0.000000000000000e-01, 1.300209931372656e+01, -6.841393840416191e+03,
      1.193492661387020e+06, -1.033456200764548e+08, 5.319778623982719e+09,
      -1.805161947596514e+11, 4.317533185036128e+12, -7.626509475155619e+13,
      1.029508945956098e+15, -1.089906772387182e+16, 9.232461935811958e+16,
      -6.357336264452389e+17, 3.603376239255576e+18, -1.698055636889742e+19,
      6.705347306221760e+19, -2.232368248073979e+20, 6.294452022548759e+20,
      -1.507836188614350e+21, 3.074157987189278e+21, -5.336595912526426e+21,
      7.880489249366001e+21, -9.874488341898865e+21, 1.045462363245188e+22,
      -9.293378659359744e+21, 6.873520006757566e+21, -4.176636208226636e+21,
      2.048299449271481e+21, -7.902800175855315e+20, 2.308396453521622e+20,
      -4.796514797886738e+19, 6.315072588841990e+18, -3.958864781209368e+17,
      0.000000000000000e-01, -1.187480374787823e+01, 6.249975469245804e+03,
      -1.090926975816432e+06, 9.454341715249161e+07, -4.872094426264295e+09,
      1.655534657183789e+11, -3.966168303111664e+12, 7.019129535830873e+13,
      -9.495382387492220e+14, 1.007610862599221e+16, -8.557186874129430e+16,
      5.908530248892806e+17, -3.358744210862016e+18, 1.587616779782261e+19,
      -6.289207145157216e+19, 2.100713184820985e+20, -5.943219992524428e+20,
      1.428595119033459e+21, -2.922754970902040e+21, 5.091600520853282e+21,
      -7.545231007708337e+21, 9.487734903251288e+21, -1.008041543577854e+22,
      8.991956191593799e+21, -6.673509171379833e+21, 4.068893946683280e+21,
      -2.002141237919232e+21, 7.750111453424084e+20, -2.271093028431890e+20,
      4.733888021452981e+19, -6.251826041286116e+18, 3.931001229299567e+17,
      0.000000000000000e-01, 1.087774446529347e+01, -5.726532522119743e+03,
      1.000026921141720e+06, -8.672640379209894e+07, 4.473426627332310e+09,
      -1.521830785057648e+11, 3.650893458666403e+12, -6.471493237339299e+13,
      8.770338012270925e+14, -9.325329559951301e+15, 7.936874744339234e+16,
      -5.493120646406107e+17, 3.130441935246294e+18, -1.483627518494573e+19,
      5.893595494590624e+19, -1.974260043524307e+20, 5.602136541583943e+20,
      -1.350733620606971e+21, 2.772103374812838e+21, -4.844503527518292e+21,
      7.202129004924270e+21, -9.085610385118606e+21, 9.684512731454793e+21,
      -8.666845324741530e+21, 6.453034816508831e+21, -3.947120866745059e+21,
      1.948412902571092e+21, -7.565912375504262e+20, 2.224014019524002e+20,
      -4.649964958533596e+19, 6.159502112364614e+18, -3.884411477495515e+17,
      0.000000000000000e-01, -9.986140223631332e+00, 5.258180249016786e+03,
      -9.185985927746997e+05, 7.971153565650261e+07, -4.114819554974130e+09,
      1.401203840137855e+11, -3.365432249133715e+12, 5.973558126785002e+13,
      -8.107918476058353e+14, 8.635671857055418e+15, -7.363618322103649e+16,
      5.106667921634938e+17,  -2.916510065095704e+18, 1.385415474244305e+19,
      -5.516783142970931e+19, 1.852714215024819e+20, -5.271074464383260e+20,
      1.274366202537370e+21, -2.622681385943975e+21, 4.596469304424025e+21,
      -6.853262795168416e+21, 8.671008970964122e+21, -9.270120999117532e+21,
      8.320885075782802e+21, -6.214097196642535e+21, 3.812422050430118e+21,
      -1.887580884371838e+21, 7.351640810621922e+20, -2.167456694682573e+20,
      4.545079901812915e+19, -6.038139340406067e+18, 3.818855272205113e+17,
      0.000000000000000e-01, 9.179865311132163e+00, -4.834435688158695e+03,
      8.448504512588955e+05, -7.334848338052211e+07, 3.788859743400152e+09,
      -1.291272970868570e+11, 3.104465491388026e+12, -5.516672358352230e+13,
      7.497538905042167e+14, -7.997160527339879e+15, 6.830050930752028e+16,
      -4.744858033268795e+17, 2.714931990926294e+18, -1.292227727069609e+19,
      5.156543365822560e+19, -1.735567332289304e+20, 4.949202035461436e+20,
      -1.199422235955881e+21, 2.474571823687724e+21, -4.347969142753134e+21,
      6.499709691677232e+21, -8.245627749154246e+21, 8.839266680637617e+21,
      -7.955961175047700e+21, 5.958068545674167e+21, -3.665569136784129e+21,
      1.819971125502232e+21, -7.108274904080239e+20, 2.101605178572964e+20,
      -4.419368247642274e+19, 5.887550238778617e+18, -3.733965028540281e+17,
      0.000000000000000e-01, -8.442157387023658e+00, 4.446553379007233e+03,
      -7.772828492878541e+05, 6.751075382325961e+07, -3.489264209336539e+09,
      1.190001385182876e+11, -2.863388547750731e+12, 5.093235749821151e+13,
      -6.929732092529343e+14, 7.400674318693824e+15, -6.329249553366740e+16,
      4.403495020338399e+17, -2.523657530207839e+18, 1.203252080967825e+19,
      -4.810263202989546e+19, 1.622139284826088e+20, -4.635104707077193e+20,
      1.125673633860567e+21, -2.327511860731306e+21, 4.098851177315484e+21,
      -6.141619421494262e+21, 7.810022020551529e+21, -8.392815674436741e+21,
      7.572989471395764e+21, -5.685659534933809e+21, 3.506969466398297e+21,
      -1.745750145592470e+21, 6.836250778812437e+20, -2.026505202012846e+20,
      4.272714338022826e+19, -5.707256190142981e+18, 3.629208802827826e+17,
      0.000000000000000e-01, 7.758620432176431e+00, -4.087010780643976e+03,
      7.146017483117621e+05, -6.208866298624787e+07, 3.210548205967185e+09,
      -1.095595506626254e+11, 2.638102073062174e+12, -4.696390598303079e+13,
      6.395814992099562e+14, -6.837680401038084e+15, 5.854580721511378e+16,
      -4.078439182746133e+17, 2.340589674143266e+18, -1.117619206099512e+19,
      4.474976818167516e+19, -1.511594767319364e+20, 4.326839217130804e+20,
      -1.052747833818875e+21, 2.180916376866737e+21, -3.848371023940319e+21,
      5.778239784350781e+21, -7.363608932672221e+21, 7.930445122237345e+21,
      -7.171862635842173e+21, 5.396860454297140e+21, -3.336620071079666e+21,
      1.664898414527208e+21, -6.535348447882521e+20, 1.942028797566696e+20,
      -4.104676746515045e+19, 5.496390689251412e+18, -3.503832522669480e+17,
      0.000000000000000e-01,  -7.116397718865410e+00, 3.749079648317452e+03,
      -6.556462179536947e+05, 5.698334876317732e+07, -2.947736407946278e+09,
      1.006414850064029e+11, -2.424817814258759e+12, 4.319719539696606e+13,
      -5.887538442457942e+14, 6.299924892323114e+15, -5.399488828717867e+16,
      3.765493816021260e+17, -2.163536905559044e+18, 1.034386803998833e+19,
      -4.147323866260055e+19, 1.402934443266815e+20, -4.021917192563760e+20,
      9.801245774469403e+20, -2.033870287254597e+21, 3.595172619362371e+21,
      -5.407874926278053e+21, 6.904593808205439e+21, -7.450539640286401e+21,
      6.751333782133479e+21, -5.090837485270753e+21, 3.154034661016916e+21,
      -1.577170275266693e+21, 6.204523939342194e+20, -1.847822507348537e+20,
      3.914377458647115e+19, -5.253552629244530e+18, 3.356769581362759e+17,
      0.000000000000000e-01, 6.503405049947615e+00, -3.426426844095358e+03,
      5.993202798185401e+05, -5.210105929407407e+07, 2.696081626301935e+09,
      -9.208817752140222e+10, 2.219856964226088e+12, -3.956916806058024e+13,
      5.396682472814272e+14, -5.779046608150560e+15, 4.957204191401071e+16,
      -3.460227309661112e+17, 1.990124376198793e+18, -9.525027063011665e+18,
      3.823419916527245e+19, -1.294955870987350e+20, 3.717202435921960e+20,
      -9.071121012468153e+20, 1.885079625280803e+21, -3.337199379905083e+21,
      5.027744029565647e+21, -6.429775852763725e+21, 6.949963688436526e+21,
      -6.308792489007701e+21, 4.765750352591050e+21, -2.958122806424166e+21,
      1.482030100296912e+21, -5.841647400501416e+20, 1.743227189970331e+20,
      -3.700329724016468e+19, 4.976575581854351e+18, -3.186495777814819e+17,
      0.000000000000000e-01, -5.907497514106532e+00, 3.112678595829013e+03,
      -5.445178292388210e+05, 4.734677219183480e+07, -2.450744429063559e+09,
      8.373760838189771e+10, -2.019407140707212e+12, 3.601376581781814e+13,
      -4.914525865076907e+14, 5.266042927583444e+15, -4.520313655093684e+16,
      3.157692701941383e+17, -1.817642904902012e+18, 8.707370090912989e+18,
      -3.498599158499494e+19, 1.186170493957938e+20, -3.408681448268105e+20,
      8.327925173944575e+20, -1.732759335262471e+21, 3.071495239478500e+21,
      -4.633677555133575e+21, 5.934150774303194e+21, -6.423619232920034e+21,
      5.839849717449368e+21, -4.418427829351906e+21, 2.746981777317360e+21,
      -1.378544874829898e+21, 5.443069194938188e+20, -1.627146166161762e+20,
      3.460155133741941e+19, -4.662146280112927e+18, 2.990786503802649e+17,
      0.000000000000000e-01, 5.315369270978286e+00, -2.800843064095075e+03,
      4.900224139921877e+05, -4.261558203433840e+07, 2.206354210150574e+09,
      -7.540878769771313e+10, 1.819175365822707e+12, -3.245589496559513e+13,
      4.431044900769581e+14, -4.750435903707279e+15, 4.080066038394767e+16,
      -2.851956961225284e+17, 1.642785899790599e+18, -7.875595013484737e+18,
      3.166934141851351e+19, -1.074644294009643e+20, 3.091013384752239e+20,
      -7.559132271023600e+20, 1.574408999942342e+21, -2.793807988826481e+21,
      4.219517259473626e+21, -5.410137053845813e+21, 5.863599367756335e+21,
      -5.337558212358710e+21, 4.043769184497504e+21, -2.517518733791446e+21,
      1.265191806068099e+21, -5.002850262989436e+20, 1.497812681253764e+20,
      -3.190085448905626e+19, 4.305131059057072e+18, -2.766285114923478e+17,
      0.000000000000000e-01, -4.710772351834030e+00, 2.482373368683492e+03,
      -4.343438634071739e+05, 3.777856440947868e+07, -1.956282178028853e+09,
      6.687710881826573e+10, -1.613799071688983e+12, 2.880102840436904e+13,
      -3.933509953446171e+14, 4.218785748675981e+15, -3.625111116519231e+16,
      2.535230396849270e+17, -1.461153893751222e+18, 7.009048271913450e+18,
      -2.820301208876073e+19, 9.576835312471221e+19, -2.756632919785344e+20,
      6.746687832309763e+20, -1.406360250100074e+21, 2.497787658823857e+21,
      -3.775905442111079e+21, 4.846020652107188e+21, -5.257496778824818e+21,
      4.790865278530114e+21, -3.633565158982738e+21, 2.264711978580421e+21,
      -1.139484606704540e+21, 4.511274997631573e+20, -1.352342175988863e+20,
      2.884004791547230e+19, -3.897279615275436e+18, 2.507669351857205e+17,
      0.000000000000000e-01, 4.070989770987400e+00, -2.145310206509185e+03,
      3.753936962817004e+05, -3.265459454227024e+07, 1.691185508550510e+09,
      -5.782472303217436e+10, 1.395652621524724e+12, -2.491398584416296e+13,
      3.403599167085817e+14, -3.651608452505086e+15, 3.138862675548854e+16,
      -2.196030665561983e+17, 1.266200972268635e+18, -6.076691812326445e+18,
      2.446363025550984e+19, -8.311520079137181e+19, 2.393790730533568e+20,
      -5.862224225217512e+20, 1.222781062900224e+21, -2.173219858365351e+21,
      3.287615107047680e+21, -4.222527251629046e+21, 4.584681175928391e+21,
      -4.181215238485294e+21, 3.173915831158933e+21, -1.979997675938663e+21,
      9.971596317435381e+20, -3.951620422561584e+20, 1.185761286177830e+20,
      -2.531374184616153e+19, 3.424415390856724e+18, -2.205841595819251e+17,
      0.000000000000000e-01, -3.358090451409268e+00, 1.769674233094550e+03,
      -3.096791496404924e+05, 2.694027470522958e+07, -1.395380783042074e+09,
      4.771664530498290e+10, -1.151859937920620e+12, 2.056566436202355e+13,
      -2.810129791321012e+14, 3.015587336970416e+15, -2.592812452685370e+16,
      1.814511218788989e+17, -1.046544742872117e+18, 5.024209499496014e+18,
      -2.023384009065594e+19, 6.877115067771995e+19, -1.981490581745830e+20,
      4.854671646250063e+20, -1.013093139325803e+21, 1.801437605201838e+21,
      -2.726610427990654e+21, 3.503909183048419e+21, -3.806619619529808e+21,
      3.473717880327268e+21, -2.638522642188920e+21, 1.647082019674793e+21,
      -8.300649678444784e+20, 3.291782934571713e+20, -9.884928188742345e+19,
      2.111857359313218e+19, -2.859159218719010e+18, 1.843238572103207e+17,
      0.000000000000000e-01, 2.488636218117664e+00, -1.311502616747683e+03,
      2.295099147207366e+05, -1.996696431094197e+07, 1.034261165808075e+09,
      -3.537054923192207e+10, 8.539117568460436e+11, -1.524770635014292e+13,
      2.083740755848528e+14, -2.236412246676083e+15, 1.923184048547382e+16,
      -1.346128075294067e+17, 7.765487558375990e+17, -3.728808938617731e+18,
      1.502032959138465e+19, -5.106384693104742e+19, 1.471677691807719e+20,
      -3.606628515829967e+20, 7.528686576360846e+20, -1.339136443296313e+21,
      2.027551807534738e+21, -2.606468672347013e+21, 2.832680005398635e+21,
      -2.585940609411422e+21, 1.964981315222245e+21, -1.227139920687813e+21,
      6.186996825719858e+20, -2.454684425809199e+20, 7.374666999744300e+19,
      -1.576324668155186e+19, 2.135204267310823e+18, -1.377242653770284e+17,
      0.000000000000000e-01, -1.000000000000000e+00, 5.270000000000000e+02,
      -9.222500000000000e+04, 8.023575000000000e+06, -4.156211850000000e+08,
      1.421424452700000e+10, -3.431724750090000e+11, 6.128079910875000e+12,
      -8.375042544862500e+13, 8.989212331485750e+14, -7.730722605077745e+15,
      5.411505823554421e+16, -3.122022590512166e+17, 1.499257002256941e+18,
      -6.039863923377964e+18, 2.053553733948508e+19, -5.919066644910405e+19,
      1.450751628654511e+20, -3.028762172103277e+20, 5.388008495636356e+20,
      -8.158984293392197e+20, 1.049012266293282e+21, -1.140230724231829e+21,
      1.041080226472539e+21, -7.912209721191298e+20, 4.942087918159488e+20,
      -2.492163992918032e+20, 9.889539654436636e+19, -2.971733590742043e+19,
      6.353361469862300e+18, -8.607780055942471e+17, 5.553406487704820e+16 };
 
  static const double *fem_coeff_gausslob[] =
    { 0, fem_coef_gausslob_1, fem_coef_gausslob_2, fem_coef_gausslob_3,
      fem_coef_gausslob_4, fem_coef_gausslob_5, fem_coef_gausslob_6, fem_coef_gausslob_7, fem_coef_gausslob_8,
      fem_coef_gausslob_9, fem_coef_gausslob_10, fem_coef_gausslob_11, fem_coef_gausslob_12, fem_coef_gausslob_13,
      fem_coef_gausslob_14, 0, fem_coef_gausslob_16, 0, 0, 0, 0, 0, 0, 0, fem_coef_gausslob_24, 0, 0, 0, 0, 0, 0, 0,
      fem_coef_gausslob_32, };
 
  const unsigned fem_coeff_gausslob_max_k = 33;
 
 
  class PK_GL_fem_ : public fem<base_poly> {
  public :
    PK_GL_fem_(unsigned k);
  };
 
  PK_GL_fem_::PK_GL_fem_(unsigned k) {
    cvr = bgeot::simplex_of_reference(1);
    dim_ = cvr->structure()->dim();
    is_standard_fem = is_equiv = is_pol = is_lag = true;
    es_degree = short_type(k);
    GMM_ASSERT1(k < fem_coeff_gausslob_max_k && fem_coeff_gausslob[k],
                "try another degree");
    init_cvs_node();
    std::stringstream sstr; sstr << "IM_GAUSSLOBATTO1D(" << k*2-1 << ")";
    pintegration_method gl_im = int_method_descriptor(sstr.str());
    std::vector<base_node> points(k+1);
    for (size_type i = 0; i < k+1; ++i) {
      points[i] = gl_im->approx_method()->point(i);
    }
    std::sort(points.begin(),points.end());
    for (size_type i = 0; i < k+1; ++i)
      add_node(lagrange_dof(1), points[i]);
    base_.resize(k+1);
    const double *coefs = fem_coeff_gausslob[k];
    for (size_type r = 0; r < k+1; r++) {
      base_[r] = base_poly(1,short_type(k));
      std::copy(coefs, coefs+k+1, base_[r].begin());
      coefs += k+1;
    }
  }
 
  static pfem PK_GL_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    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));
    pfem p = std::make_shared<PK_GL_fem_>(k);
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*        Hermite element on the segment                                    */
  /* ******************************************************************** */
 
  struct hermite_segment__ : public fem<base_poly> {
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    hermite_segment__();
  };
 
  void hermite_segment__::mat_trans(base_matrix &M,
                                    const base_matrix &G,
                                    bgeot::pgeometric_trans pgt) const {
    THREAD_SAFE_STATIC bgeot::pgeotrans_precomp pgp;
    THREAD_SAFE_STATIC bgeot::pgeometric_trans pgt_stored = nullptr;
    THREAD_SAFE_STATIC base_matrix K(1, 1);
    THREAD_SAFE_STATIC base_vector r(1);
    dim_type N = dim_type(G.nrows());
 
    if (pgt != pgt_stored) {
      gmm::resize(r, N);
      for (size_type i = 0; i < N; ++i) r[i] = ::exp(double(i));
      pgt_stored = pgt;
      pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
      gmm::resize(K, N, 1);
    }
    gmm::copy(gmm::identity_matrix(), M);
    // gradient at point 0
    gmm::mult(G, pgp->grad(1), K);
    if (N == 1) M(1, 1) = K(0,0);
    else  M(1, 1) = gmm::mat_euclidean_norm(K)
      * gmm::sgn(gmm::vect_sp(gmm::mat_col(K, 0), r));
    // gradient at point 1
    if (!(pgt->is_linear())) gmm::mult(G, pgp->grad(3), K);
    if (N == 1) M(3, 3) = K(0,0);
    else M(3, 3) = gmm::mat_euclidean_norm(K)
      * gmm::sgn(gmm::vect_sp(gmm::mat_col(K, 0), r));
  }
 
  // Hermite element on the segment. when the real element lies in
  // a 2 or 3 dimensional domain, the element should still work if
  // the tangent coincides.
  hermite_segment__::hermite_segment__() {
    base_node pt(1);
    cvr = bgeot::simplex_of_reference(1);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = 3;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
    base_.resize(4);
 
    pt[0] = 0.0; add_node(lagrange_dof(1), pt);
    read_poly(base_[0], 1, "(1 - x)^2*(2*x + 1)");
 
    pt[0] = 0.0; add_node(derivative_dof(1, 0), pt);
    read_poly(base_[1], 1, "x*(x - 1)*(x - 1)");
 
    pt[0] = 1.0; add_node(lagrange_dof(1), pt);
    read_poly(base_[2], 1, "x*x*(3  - 2*x)");
 
    pt[0] = 1.0; add_node(derivative_dof(1, 0), pt);
    read_poly(base_[3], 1, "x*x*(x - 1)");
  }
 
  /* ******************************************************************** */
  /*    Hermite element on the triangle                                   */
  /* ******************************************************************** */
 
  struct hermite_triangle__ : public fem<base_poly> {
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    hermite_triangle__();
  };
 
  void hermite_triangle__::mat_trans(base_matrix &M,
                                    const base_matrix &G,
                                    bgeot::pgeometric_trans pgt) const {
 
    THREAD_SAFE_STATIC bgeot::pgeotrans_precomp pgp;
    THREAD_SAFE_STATIC bgeot::pgeometric_trans pgt_stored = nullptr;
    THREAD_SAFE_STATIC base_matrix K(2, 2);
    dim_type N = dim_type(G.nrows());
 
    GMM_ASSERT1(N == 2, "Sorry, this version of hermite "
                "element works only in dimension two.")
    if (pgt != pgt_stored)
      { pgt_stored = pgt; pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0); }
    gmm::copy(gmm::identity_matrix(), M);
 
    gmm::mult(G, pgp->grad(0), K);
    for (size_type i = 0; i < 3; ++i) {
      if (i && !(pgt->is_linear())) gmm::mult(G, pgp->grad(i*3), K);
      gmm::copy(K, gmm::sub_matrix(M, gmm::sub_interval(1+3*i, 2)));
    }
  }
 
  hermite_triangle__::hermite_triangle__() {
    cvr = bgeot::simplex_of_reference(2);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = 3;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
    base_.resize(10);
 
    add_node(lagrange_dof(2), base_node(0.0, 0.0));
    read_poly(base_[0], 2, "(1 - x - y)*(1 + x + y - 2*x*x - 11*x*y - 2*y*y)");
 
    add_node(derivative_dof(2, 0), base_node(0.0, 0.0));
    read_poly(base_[1], 2, "x*(1 - x - y)*(1 - x - 2*y)");
 
    add_node(derivative_dof(2, 1), base_node(0.0, 0.0));
    read_poly(base_[2], 2, "y*(1 - x - y)*(1 - 2*x - y)");
 
    add_node(lagrange_dof(2), base_node(1.0, 0.0));
    read_poly(base_[3], 2, "-2*x*x*x + 7*x*x*y + 7*x*y*y + 3*x*x - 7*x*y");
 
    add_node(derivative_dof(2, 0), base_node(1.0, 0.0));
    read_poly(base_[4], 2, "x*x*x - 2*x*x*y - 2*x*y*y - x*x + 2*x*y");
 
    add_node(derivative_dof(2, 1), base_node(1.0, 0.0));
    read_poly(base_[5], 2, "x*y*(2*x + y - 1)");
 
    add_node(lagrange_dof(2), base_node(0.0, 1.0));
    read_poly(base_[6], 2, "7*x*x*y + 7*x*y*y - 2*y*y*y + 3*y*y - 7*x*y");
 
    add_node(derivative_dof(2, 0), base_node(0.0, 1.0));
    read_poly(base_[7], 2, "x*y*(x + 2*y - 1)");
 
    add_node(derivative_dof(2, 1), base_node(0.0, 1.0));
    read_poly(base_[8], 2, "y*y*y - 2*y*y*x - 2*y*x*x - y*y + 2*x*y");
 
    add_node(lagrange_dof(2), base_node(1.0/3.0, 1.0/3.0));
    read_poly(base_[9], 2, "27*x*y*(1 - x - y)");
 
 }
 
  /* ******************************************************************** */
  /*    Hermite element on the tetrahedron                                */
  /* ******************************************************************** */
 
  struct hermite_tetrahedron__ : public fem<base_poly> {
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    hermite_tetrahedron__();
  };
 
  void hermite_tetrahedron__::mat_trans(base_matrix &M,
                                    const base_matrix &G,
                                    bgeot::pgeometric_trans pgt) const {
    THREAD_SAFE_STATIC bgeot::pgeotrans_precomp pgp;
    THREAD_SAFE_STATIC bgeot::pgeometric_trans pgt_stored = nullptr;
    THREAD_SAFE_STATIC base_matrix K(3, 3);
    dim_type N = dim_type(G.nrows());
    GMM_ASSERT1(N == 3, "Sorry, this version of hermite "
                "element works only on dimension three.")
    if (pgt != pgt_stored)
      { pgt_stored = pgt; pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0); }
    gmm::copy(gmm::identity_matrix(), M);
 
    gmm::mult(G, pgp->grad(0), K);
    for (size_type k = 0; k < 4; ++k) {
      if (k && !(pgt->is_linear())) gmm::mult(G, pgp->grad(k*4), K);
      gmm::copy(K, gmm::sub_matrix(M, gmm::sub_interval(1+4*k, 3)));
    }
  }
 
  hermite_tetrahedron__::hermite_tetrahedron__() {
    cvr = bgeot::simplex_of_reference(3);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = 3;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
    base_.resize(20);
    std::stringstream s
      ( "1 - 3*x*x - 13*x*y - 13*x*z - 3*y*y - 13*y*z - 3*z*z + 2*x*x*x"
        "+ 13*x*x*y + 13*x*x*z + 13*x*y*y + 33*x*y*z + 13*x*z*z + 2*y*y*y"
        "+ 13*y*y*z + 13*y*z*z + 2*z*z*z;"
        "x - 2*x*x - 3*x*y - 3*x*z + x*x*x + 3*x*x*y + 3*x*x*z + 2*x*y*y"
        "+ 4*x*y*z + 2*x*z*z;"
        "y - 3*x*y - 2*y*y - 3*y*z + 2*x*x*y + 3*x*y*y + 4*x*y*z"
        "+ y*y*y + 3*y*y*z + 2*y*z*z;"
        "z - 3*x*z - 3*y*z - 2*z*z + 2*x*x*z + 4*x*y*z + 3*x*z*z"
        "+ 2*y*y*z + 3*y*z*z + z*z*z;"
        "3*x*x - 7*x*y - 7*x*z - 2*x*x*x + 7*x*x*y + 7*x*x*z + 7*x*y*y"
        "+ 7*x*y*z + 7*x*z*z;"
        "-x*x + 2*x*y + 2*x*z + x*x*x - 2*x*x*y - 2*x*x*z - 2*x*y*y"
        "- 2*x*y*z - 2*x*z*z;"
        "-x*y + 2*x*x*y + x*y*y;"
        "-x*z + 2*x*x*z + x*z*z;"
        "-7*x*y + 3*y*y - 7*y*z + 7*x*x*y + 7*x*y*y + 7*x*y*z - 2*y*y*y"
        "+ 7*y*y*z + 7*y*z*z;"
        "-x*y + x*x*y + 2*x*y*y;"
        "2*x*y - y*y + 2*y*z - 2*x*x*y - 2*x*y*y - 2*x*y*z + y*y*y"
        "- 2*y*y*z - 2*y*z*z;"
        "-y*z + 2*y*y*z + y*z*z;"
        "-7*x*z - 7*y*z + 3*z*z + 7*x*x*z + 7*x*y*z + 7*x*z*z + 7*y*y*z"
        "+ 7*y*z*z - 2*z*z*z;"
        "-x*z + x*x*z + 2*x*z*z;"
        "-y*z + y*y*z + 2*y*z*z;"
        "2*x*z + 2*y*z - z*z - 2*x*x*z - 2*x*y*z - 2*x*z*z - 2*y*y*z"
        "- 2*y*z*z + z*z*z;"
        "27*x*y*z;"
        "27*y*z - 27*x*y*z - 27*y*y*z - 27*y*z*z;"
        "27*x*z - 27*x*x*z - 27*x*y*z - 27*x*z*z;"
        "27*x*y - 27*x*x*y - 27*x*y*y - 27*x*y*z;");
 
    base_node pt(3);
    for (unsigned k = 0; k < 5; ++k) {
      for (unsigned i = 0; i < 4; ++i) {
        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;
        if (k == 4 || i == 0)  add_node(lagrange_dof(3), pt);
        else add_node(derivative_dof(3, dim_type(i-1)), pt);
      }
    }
  }
 
  static pfem Hermite_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    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 d = int(::floor(params[0].num() + 0.01));
    pfem p;
    switch(d) {
    case 1 : p = std::make_shared<hermite_segment__>(); break;
    case 2 : p = std::make_shared<hermite_triangle__>(); break;
    case 3 : p = std::make_shared<hermite_tetrahedron__>(); break;
    default : GMM_ASSERT1(false, "Sorry, Hermite element in dimension "
                          << d << " not available");
    }
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return pfem(p);
  }
 
  /* ******************************************************************** */
  /*    Argyris element on the triangle                                   */
  /* ******************************************************************** */
 
  struct argyris_triangle__ : public fem<base_poly> {
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    argyris_triangle__();
  };
 
  void argyris_triangle__::mat_trans(base_matrix &M,
                                    const base_matrix &G,
                                    bgeot::pgeometric_trans pgt) const {
 
    THREAD_SAFE_STATIC bgeot::pgeotrans_precomp pgp;
    THREAD_SAFE_STATIC pfem_precomp pfp;
    THREAD_SAFE_STATIC bgeot::pgeometric_trans pgt_stored = nullptr;
    THREAD_SAFE_STATIC base_matrix K(2, 2);
    dim_type N = dim_type(G.nrows());
    GMM_ASSERT1(N == 2, "Sorry, this version of argyris "
                "element works only on dimension two.")
    if (pgt != pgt_stored) {
      pgt_stored = pgt;
      pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
      pfp = fem_precomp(std::make_shared<argyris_triangle__>(), node_tab(0), 0);
    }
    gmm::copy(gmm::identity_matrix(), M);
 
    gmm::mult(G, pgp->grad(0), K); // gradient at point (0, 0)
    for (unsigned k = 0; k < 3; ++k) {
      if (k && !(pgt->is_linear())) gmm::mult(G, pgp->grad(6*k), K);
      M(1+6*k, 1+6*k) = K(0,0); M(1+6*k, 2+6*k) = K(0,1);
      M(2+6*k, 1+6*k) = K(1,0); M(2+6*k, 2+6*k) = K(1,1);
      if (!(pgt->is_linear())) {
        base_matrix XX[2], H(2,4), B(2,2), X(2,2);
        XX[0] = XX[1] = base_matrix(2,2);
        gmm::copy(gmm::transposed(K), B); gmm::lu_inverse(B);
        gmm::mult(G, pgp->hessian(6*k), H);
        for (unsigned j = 0; j < 2; ++j) {
          XX[j](0,0) = B(0, j)*H(0, 0) + B(1, j)*H(1, 0);
          XX[j](0,1) = XX[j](1,0) = B(0, j)*H(0, 1) + B(1, j)*H(1, 1);
          XX[j](1,1) = B(0, j)*H(0, 3) + B(1, j)*H(1, 3);
        }
        for (unsigned j = 0; j < 2; ++j) {
          gmm::copy(gmm::scaled(XX[0], K(j,0)), X);
          gmm::add(gmm::scaled(XX[1], K(j,1)), X);
          M(1+j+6*k, 3+6*k) = X(0,0); M(1+j+6*k, 4+6*k) = X(1, 0);
          M(1+j+6*k, 5+6*k) = X(1, 1);
        }
      }
      scalar_type a = K(0,0), b = K(0,1), c = K(1,0), d = K(1,1);
      M(3+6*k, 3+6*k) = a*a;     M(3+6*k, 4+6*k) = a*b;       M(3+6*k, 5+6*k) = b*b;
      M(4+6*k, 3+6*k) = 2.0*a*c; M(4+6*k, 4+6*k) = b*c + a*d; M(4+6*k, 5+6*k) = 2.0*b*d;
      M(5+6*k, 3+6*k) = c*c;     M(5+6*k, 4+6*k) = c*d;       M(5+6*k, 5+6*k) = d*d;
    }
 
    THREAD_SAFE_STATIC base_matrix W(3, 21);
    base_small_vector norient(M_PI, M_PI * M_PI);
    if (pgt->is_linear()) gmm::lu_inverse(K);
    for (unsigned i = 18; i < 21; ++i) {
      if (!(pgt->is_linear()))
        { gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
      bgeot::base_small_vector n(2), v(2);
      gmm::mult(gmm::transposed(K), cvr->normals()[i-18], n);
      n /= gmm::vect_norm2(n);
 
      scalar_type ps = gmm::vect_sp(n, norient);
      if (ps < 0) n *= scalar_type(-1);
      if (gmm::abs(ps) < 1E-8)
        GMM_WARNING2("Argyris : The normal orientation may be incorrect");
      gmm::mult(K, n, v);
      const bgeot::base_tensor &t = pfp->grad(i);
      for (unsigned j = 0; j < 21; ++j)
        W(i-18, j) = t(j, 0, 0) * v[0] + t(j, 0, 1) * v[1];
    }
    THREAD_SAFE_STATIC base_matrix A(3,3);
    THREAD_SAFE_STATIC bgeot::base_vector w(3);
    THREAD_SAFE_STATIC bgeot::base_vector coeff(3);
    THREAD_SAFE_STATIC gmm::sub_interval SUBI(18, 3);
    THREAD_SAFE_STATIC gmm::sub_interval SUBJ(0, 3);
    gmm::copy(gmm::sub_matrix(W, SUBJ, SUBI), A);
    gmm::lu_inverse(A);
    gmm::copy(gmm::transposed(A), gmm::sub_matrix(M, SUBI));
 
    for (unsigned j = 0; j < 18; ++j) {
      gmm::mult(W, gmm::mat_row(M, j), w);
      gmm::mult(A, gmm::scaled(w, -1.0), coeff);
      gmm::copy(coeff, gmm::sub_vector(gmm::mat_row(M, j), SUBI));
    }
  }
 
  argyris_triangle__::argyris_triangle__() {
    cvr = bgeot::simplex_of_reference(2);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = 5;
    is_pol = true;
    is_lag = false;
    is_standard_fem = is_equiv = false;
    base_.resize(21);
 
    std::stringstream s
      ("1 - 10*x^3 - 10*y^3 + 15*x^4 - 30*x*x*y*y"
       "+ 15*y*y*y*y - 6*x^5 + 30*x*x*x*y*y + 30*x*x*y*y*y - 6*y^5;"
       "x - 6*x*x*x - 11*x*y*y + 8*x*x*x*x + 10*x*x*y*y"
       "+ 18*x*y*y*y - 3*x*x*x*x*x + x*x*x*y*y - 10*x*x*y*y*y - 8*x*y*y*y*y;"
       "y - 11*x*x*y - 6*y*y*y + 18*x*x*x*y + 10*x*x*y*y"
       "+ 8*y*y*y*y - 8*x*x*x*x*y - 10*x*x*x*y*y + x*x*y*y*y - 3*y*y*y*y*y;"
       "0.5*x*x - 1.5*x*x*x + 1.5*x*x*x*x - 1.5*x*x*y*y"
       "- 0.5*x*x*x*x*x + 1.5*x*x*x*y*y + x*x*y*y*y;"
       "x*y - 4*x*x*y - 4*x*y*y + 5*x*x*x*y + 10*x*x*y*y"
       "+ 5*x*y*y*y - 2*x*x*x*x*y - 6*x*x*x*y*y - 6*x*x*y*y*y - 2*x*y*y*y*y;"
       "0.5*y*y - 1.5*y*y*y - 1.5*x*x*y*y + 1.5*y*y*y*y + x*x*x*y*y"
       "+ 1.5*x*x*y*y*y - 0.5*y*y*y*y*y;"
       "10*x^3 - 15*x^4 + 15*x*x*y*y + 6*x^5 - 15*x*x*x*y*y - 15*x*x*y*y*y;"
       "-4*x*x*x + 7*x*x*x*x - 3.5*x*x*y*y - 3*x*x*x*x*x + 3.5*x*x*x*y*y"
       "+ 3.5*x*x*y*y*y;"
       "-5*x*x*y + 14*x*x*x*y + 18.5*x*x*y*y - 8*x*x*x*x*y"
       "- 18.5*x*x*x*y*y - 13.5*x*x*y*y*y;"
       "0.5*x*x*x - x*x*x*x + 0.25*x*x*y*y + 0.5*x*x*x*x*x"
       "- 0.25*x*x*x*y*y - 0.25*x*x*y*y*y;"
       "x*x*y - 3*x*x*x*y - 3.5*x*x*y*y + 2*x*x*x*x*y + 3.5*x*x*x*y*y"
       "+ 2.5*x*x*y*y*y;"
       "1.25*x*x*y*y - 0.75*x*x*x*y*y - 1.25*x*x*y*y*y;"
       "10*y*y*y + 15*x*x*y*y - 15*y^4 - 15*x*x*x*y*y - 15*x*x*y*y*y + 6*y^5;"
       "-5*x*y*y + 18.5*x*x*y*y + 14*x*y*y*y - 13.5*x*x*x*y*y"
       "- 18.5*x*x*y*y*y - 8*x*y*y*y*y;"
       "-4*y*y*y - 3.5*x*x*y*y + 7*y*y*y*y + 3.5*x*x*x*y*y"
       "+ 3.5*x*x*y*y*y - 3*y*y*y*y*y;"
       "1.25*x*x*y*y - 1.25*x*x*x*y*y - 0.75*x*x*y*y*y;"
       "x*y*y - 3.5*x*x*y*y - 3*x*y*y*y + 2.5*x*x*x*y*y + 3.5*x*x*y*y*y"
       "+ 2*x*y*y*y*y;"
       "0.5*y*y*y + 0.25*x*x*y*y - y*y*y*y - 0.25*x*x*x*y*y"
       "- 0.25*x*x*y*y*y + 0.5*y*y*y*y*y;"
       "sqrt(2) * (-8*x*x*y*y + 8*x*x*x*y*y + 8*x*x*y*y*y);"
       "-16*x*y*y + 32*x*x*y*y + 32*x*y*y*y - 16*x*x*x*y*y"
       "- 32*x*x*y*y*y - 16*x*y*y*y*y;"
       "-16*x*x*y + 32*x*x*x*y + 32*x*x*y*y - 16*x*x*x*x*y"
       "- 32*x*x*x*y*y - 16*x*x*y*y*y;");
 
    base_node pt(2);
    for (unsigned k = 0; k < 7; ++k) {
      for (unsigned i = 0; i < 3; ++i) {
        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);
        }
        else {
          pt[0] = pt[1] = 0.0; if (k/2) pt[k/2-1] = 1.0;
          if (k & 1)
            add_node(second_derivative_dof(2, dim_type((i) ? 1:0),
                                           dim_type((i == 2) ? 1:0)), pt);
          else {
            if (i) add_node(derivative_dof(2, dim_type(i-1)), pt);
            else add_node(lagrange_dof(2), pt);
          }
        }
      }
    }
  }
 
  static pfem triangle_Argyris_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
    pfem p = std::make_shared<argyris_triangle__>();
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*    Morley element on the triangle and tetrahedron                    */
  /* ******************************************************************** */
 
  struct morley_element__ : public fem<base_poly> {
    virtual void mat_trans(base_matrix &M, const base_matrix &G,
                           bgeot::pgeometric_trans pgt) const;
    morley_element__(dim_type);
  };
 
  void morley_element__::mat_trans(base_matrix &M,
                                    const base_matrix &G,
                                    bgeot::pgeometric_trans pgt) const {
 
    THREAD_SAFE_STATIC bgeot::pgeotrans_precomp pgp;
    THREAD_SAFE_STATIC pfem_precomp pfp;
    THREAD_SAFE_STATIC bgeot::pgeometric_trans pgt_stored = nullptr;
    dim_type n = dim();
    dim_type N = dim_type(G.nrows());
    size_type nddl = base_.size();
    GMM_ASSERT1(N == n, "Sorry, this version of morley "
                "element works only for indentical dimensions.")
 
    if (pgt != pgt_stored) {
      pgt_stored = pgt;
      pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
      pfp = fem_precomp(std::make_shared<morley_element__>(n), node_tab(0), 0);
    }
    gmm::copy(gmm::identity_matrix(), M);
 
    THREAD_SAFE_STATIC base_matrix W(n+1, nddl);
    dim_type nfd = (n == 2) ? 3 : 6;
    
    base_matrix K(n, n);
    base_small_vector norient(n);
    norient[0] = M_PI;
    for (dim_type i = 1; i < n; ++i)
      norient[i] = M_PI * norient[i-1];
    if (pgt->is_linear())
      { gmm::mult(G, pgp->grad(0), K); gmm::lu_inverse(K); }
    for (unsigned i = nfd; i < nddl; ++i) {
      if (!(pgt->is_linear()))
        { gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
      bgeot::base_small_vector nn(n), v(n);
      gmm::mult(gmm::transposed(K), cvr->normals()[i-nfd], nn);
      nn /= gmm::vect_norm2(nn);
 
      scalar_type ps = gmm::vect_sp(nn, norient);
      if (ps < 0) nn *= scalar_type(-1);
      if (gmm::abs(ps) < 1E-8)
        GMM_WARNING2("Morley : The normal orientation may be incorrect");
      gmm::mult(K, nn, v);
      const bgeot::base_tensor &t = pfp->grad(i);
      for (unsigned j = 0; j < nddl; ++j)
        if (n == 2)
          W(i-nfd, j) = t(j, 0, 0) * v[0] + t(j, 0, 1) * v[1];
        else
          W(i-nfd, j) = t(j, 0, 0) * v[0] + t(j, 0, 1) * v[1] + t(j, 0, 2)*v[2];
    }
 
    base_matrix A(n+1, n+1);
    base_vector w(n+1);
    base_vector coeff(n+1);
    gmm::sub_interval SUBI(nfd, n+1);
    gmm::sub_interval SUBJ(0, n+1);
    gmm::copy(gmm::sub_matrix(W, SUBJ, SUBI), A);
    gmm::lu_inverse(A);
    gmm::copy(gmm::transposed(A), gmm::sub_matrix(M, SUBI));
 
    for (unsigned j = 0; j < nfd; ++j) {
      gmm::mult(W, gmm::mat_row(M, j), w);
      gmm::mult(A, gmm::scaled(w, -1.0), coeff);
      gmm::copy(coeff, gmm::sub_vector(gmm::mat_row(M, j), SUBI));
    }
  }
 
  morley_element__::morley_element__(dim_type n) {
    cvr = bgeot::simplex_of_reference(n);
    dim_ = cvr->structure()->dim();
    init_cvs_node();
    es_degree = 2;
    is_pol = true;
    is_standard_fem = is_lag = is_equiv = false;
 
    if (n == 2) {
      base_.resize(6);
 
      std::stringstream s("1 - x - y + 2*x*y;  (x + y + x^2 - 2*x*y - y^2)/2;"
                          "(x + y - x^2 - 2*x*y + y^2)/2;"
                          "((x+y)^2 - x - y)*sqrt(2)/2;  x*(x-1);  y*(y-1);");
 
      for (unsigned k = 0; k < 6; ++k)
        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));
      add_node(lagrange_dof(2), base_node(0.0, 1.0));
      add_node(normal_derivative_dof(2), base_node(0.5, 0.5));
      add_node(normal_derivative_dof(2), base_node(0.0, 0.5));
      add_node(normal_derivative_dof(2), base_node(0.5, 0.0));
    } else {
      base_.resize(10);
 
      std::stringstream s
        ("5/3 - 16/9*x - 28/9*y - 28/9*z + 8/9*x^2 + 8/3*x*y + 8/3*x*z"
         " - 4/9*y^2 + 20/3*y*z - 4/9*z^2;"
         "5/3 - 28/9*x - 16/9*y - 28/9*z - 4/9*x^2 + 8/3*x*y + 20/3*x*z"
         " + 8/9*y^2 + 8/3*y*z - 4/9*z^2;"
         "5/3 - 28/9*x - 28/9*y - 16/9*z - 4/9*x^2 + 20/3*x*y + 8/3*x*z"
         " - 4/9*y^2 + 8/3*y*z + 8/9*z^2;"
         "-4/3 + 20/9*x + 20/9*y + 32/9*z + 8/9*x^2 - 4/3*x*y - 16/3*x*z"
         " + 8/9*y^2 - 16/3*y*z - 16/9*z^2;"
         "-4/3 + 20/9*x + 32/9*y + 20/9*z + 8/9*x^2 - 16/3*x*y - 4/3*x*z"
         " - 16/9*y^2 - 16/3*y*z + 8/9*z^2;"
         "-4/3 + 32/9*x + 20/9*y + 20/9*z - 16/9*x^2 - 16/3*x*y - 16/3*x*z"
         " + 8/9*y^2 - 4/3*y*z + 8/9*z^2;"
         "sqrt(3)*(3/8 - x - y - z + 1/2*x^2 + 3/2*x*y + 3/2*x*z + 1/2*y^2"
         " + 3/2*y*z + 1/2*z^2);"
         "-1/8 - 2/3*x + 1/3*y + 1/3*z + 11/6*x^2 - 1/2*x*y - 1/2*x*z"
         " - 1/6*y^2 - 1/2*y*z - 1/6*z^2;"
         "-1/8 + 1/3*x - 2/3*y + 1/3*z - 1/6*x^2 - 1/2*x*y - 1/2*x*z"
         " + 11/6*y^2 - 1/2*y*z - 1/6*z^2;"
         "-1/8 + 1/3*x + 1/3*y - 2/3*z - 1/6*x^2 - 1/2*x*y - 1/2*x*z"
         " - 1/6*y^2 - 1/2*y*z + 11/6*z^2;");
      
      for (unsigned k = 0; k < 10; ++k)
        base_[k] = read_base_poly(3, s);
 
      add_node(lagrange_dof(3), base_node(0.5, 0.0, 0.0));
      add_node(lagrange_dof(3), base_node(0.0, 0.5, 0.0));
      add_node(lagrange_dof(3), base_node(0.0, 0.0, 0.5));
      add_node(lagrange_dof(3), base_node(0.5, 0.5, 0.0));
      add_node(lagrange_dof(3), base_node(0.5, 0.0, 0.5));
      add_node(lagrange_dof(3), base_node(0.0, 0.5, 0.5));
      add_node(normal_derivative_dof(3), base_node(1./3., 1./3., 1./3.));
      add_node(normal_derivative_dof(3), base_node(0.0, 1./3., 1./3.));
      add_node(normal_derivative_dof(3), base_node(1./3., 0.0, 1./3.));
      add_node(normal_derivative_dof(3), base_node(1./3., 1./3., 0.0));
    }
  }
 
  static pfem element_Morley_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    GMM_ASSERT1(params.size() == 1, "Bad number of parameters");
    int n = dim_type(::floor(params[0].num() + 0.01));
    GMM_ASSERT1(n > 0 && n < 4, "Bad parameters");
    pfem p = std::make_shared<morley_element__>(dim_type(n));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*    DISCONTINUOUS PK                                                  */
  /* ******************************************************************** */
 
  struct PK_discont_ : public PK_fem_ {
  public :
 
    PK_discont_(dim_type nc, short_type k, scalar_type alpha=scalar_type(0))
      : PK_fem_(nc, k) {
 
      std::fill(dof_types_.begin(),
                dof_types_.end(), lagrange_nonconforming_dof(nc));
     
      if (alpha != scalar_type(0)) {
        base_node G =
          gmm::mean_value(cv_node.points().begin(), cv_node.points().end());
        for (size_type i=0; i < cv_node.nb_points(); ++i)
          cv_node.points()[i] = (1-alpha)*cv_node.points()[i] + alpha*G;
        for (size_type d = 0; d < nc; ++d) {
          base_poly S(1,2);
          S[0] = -alpha * G[d] / (1-alpha);
          S[1] = 1. / (1-alpha);
          for (size_type j=0; j < nb_base(0); ++j)
            base_[j] = bgeot::poly_substitute_var(base_[j],S,d);
        }
      }
    }
  };
 
  static pfem PK_discontinuous_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies) {
    GMM_ASSERT1(params.size() == 2 || params.size() == 3,
                "Bad number of parameters : "
                << params.size() << " should be 2 or 3.");
    GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0 &&
                (params.size() != 3 || params[2].type() == 0),
                "Bad type of parameters");
    int n = int(::floor(params[0].num() + 0.01));
    int k = int(::floor(params[1].num() + 0.01));
    scalar_type alpha = 0.;
    if (params.size() == 3) alpha = params[2].num();
    GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 && alpha >= 0 &&
                alpha < 1 && double(n) == params[0].num()
                && double(k) == params[1].num(), "Bad parameters");
    pfem p = std::make_shared<PK_discont_>(dim_type(n), short_type(k), alpha);
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*    Connectable interior PK                                           */
  /* ******************************************************************** */
 
  struct PK_conn_int_ : public PK_fem_ {
  public :
 
    PK_conn_int_(dim_type nc, short_type k)
      : PK_fem_(nc, k) {
      scalar_type alpha = 0.1;
 
      std::fill(dof_types_.begin(), dof_types_.end(), lagrange_dof(nc));
      if (alpha != scalar_type(0)) {
        base_node G =
          gmm::mean_value(cv_node.points().begin(), cv_node.points().end());
        for (size_type i=0; i < cv_node.nb_points(); ++i)
          cv_node.points()[i] = (1-alpha)*cv_node.points()[i] + alpha*G;
        for (size_type d = 0; d < nc; ++d) {
          base_poly S(1,2);
          S[0] = -alpha * G[d] / (1-alpha);
          S[1] = 1. / (1-alpha);
          for (size_type j=0; j < nb_base(0); ++j)
            base_[j] = bgeot::poly_substitute_var(base_[j],S,d);
        }
      }
    }
  };
 
  static pfem conn_int_PK_fem(fem_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() == 0 && params[1].type() == 0,
                "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");
    pfem p = std::make_shared<PK_conn_int_>(dim_type(n), short_type(k));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*    Interior PK                                                       */
  /* ******************************************************************** */
  /* Interior PK element for arbitrary shape                              */
  /* (equivalent to DISCONTINUOUS_PK(n,k,0.1) for simplicies              */
 
  
  struct PK_int_ : public PK_fem_ {
  public :
 
    PK_int_(dim_type nc, short_type k, bgeot::pconvex_ref cvr_)
      : PK_fem_(nc, k) {
      cvr = cvr_;
      
      scalar_type alpha = 0.1;
      std::fill(dof_types_.begin(),
                  dof_types_.end(), lagrange_nonconforming_dof(nc));
 
      if (alpha != scalar_type(0)) {
        base_node G =
          gmm::mean_value(cv_node.points().begin(), cv_node.points().end());
        for (size_type i=0; i < cv_node.nb_points(); ++i)
          cv_node.points()[i] = (1-alpha)*cv_node.points()[i] + alpha*G;
        for (size_type d = 0; d < nc; ++d) {
          base_poly S(1,2);
          S[0] = -alpha * G[d] / (1-alpha);
          S[1] = 1. / (1-alpha);
          for (size_type j=0; j < nb_base(0); ++j)
            base_[j] = bgeot::poly_substitute_var(base_[j],S,d);
        }
      }
    }
  };
 
  static pfem simplex_IPK_fem(fem_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() == 0 && params[1].type() == 0,
                "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");
    pfem p = std::make_shared<PK_int_>(dim_type(n), short_type(k),
                                     bgeot::simplex_of_reference(dim_type(n)));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  static pfem parallelepiped_IPK_fem(fem_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() == 0 && params[1].type() == 0,
                "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");
    pfem p = std::make_shared<PK_int_>(dim_type(n), short_type(k),
                              bgeot::parallelepiped_of_reference(dim_type(n)));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  static pfem prism_IPK_fem(fem_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() == 0 && params[1].type() == 0,
                "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");
    pfem p = std::make_shared<PK_int_>(dim_type(n), short_type(k),
                                       bgeot::prism_of_reference(dim_type(n)));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
 
  /* ******************************************************************** */
  /*    PK element with a bubble base function                            */
  /* ******************************************************************** */
 
  struct PK_with_cubic_bubble_ : public PK_fem_ {
    PK_with_cubic_bubble_(dim_type nc, short_type k);
  };
 
  PK_with_cubic_bubble_::PK_with_cubic_bubble_(dim_type nc, short_type k)
    : PK_fem_(nc, k) {
    unfreeze_cvs_node();
    is_lag = false; es_degree = short_type(nc+1);
    base_node pt(nc);
    size_type j;
    PK_fem_ P1(nc, 1);
 
    pt.fill(1./(nc+1)); /* barycenter of the convex */
 
    add_node(bubble1_dof(nc), pt);
    base_.resize(nb_dof(0));
 
    j = nb_dof(0) - 1;
    base_[j] = base_poly(nc, 0);
    base_[j].one();
    for (size_type i = 0; i < P1.nb_dof(0); i++) base_[j] *= P1.base()[i];
  }
 
  static pfem PK_with_cubic_bubble(fem_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() == 0 && params[1].type() == 0,
                "Bad type of parameters");
    int n = int(::floor(params[0].num() + 0.01));
    int k = int(::floor(params[1].num() + 0.01));
    GMM_ASSERT1(k < n+1, "dimensions mismatch");
    GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
                double(n) == params[0].num() && double(k) == params[1].num(),
                "Bad parameters");
    pfem p = std::make_shared<PK_with_cubic_bubble_>(dim_type(n),
                                                     short_type(k));
    dependencies.push_back(p->ref_convex(0));
    dependencies.push_back(p->node_tab(0));
    return p;
  }
 
  /* ******************************************************************** */
  /*    classical fem                                                     */
  /* ******************************************************************** */
 
  static pfem classical_fem_(bgeot::pgeometric_trans pgt,
                             short_type k, bool complete=false,
                             bool discont=false, scalar_type alpha=0) {
    THREAD_SAFE_STATIC bgeot::pgeometric_trans pgt_last = nullptr;
    THREAD_SAFE_STATIC short_type k_last = short_type(-1);
    THREAD_SAFE_STATIC pfem fem_last = nullptr;
    THREAD_SAFE_STATIC bool complete_last = 0;
    THREAD_SAFE_STATIC bool discont_last = 0;
    THREAD_SAFE_STATIC scalar_type alpha_last = 0;
 
    bool found = false;
    if (pgt_last == pgt && k_last == k && complete_last == complete &&
        discont_last == discont && alpha_last == alpha)
      return fem_last;
 
    dim_type n = pgt->structure()->dim();
    dim_type nbp = dim_type(pgt->basic_structure()->nb_points());
    std::stringstream name;
    std::string suffix(discont ? "_DISCONTINUOUS" : "");
    GMM_ASSERT2(discont || alpha == scalar_type(0),
                "Cannot use an alpha parameter in continuous elements.");
    std::stringstream arg_;
    if (discont && alpha != scalar_type(0))
      arg_ << "," << alpha;
    std::string arg(arg_.str());
 
    // Identifying if it is a torus structure
    if (bgeot::is_torus_geom_trans(pgt) && n == 3) n = 2;
 
    /* Identifying P1-simplexes.                                          */
    if (nbp == n+1)
      if (pgt->basic_structure() == bgeot::simplex_structure(n)) {
        name << "FEM_PK" << suffix << "(" << n << "," << k << arg << ")";
        found = true;
      }
 
    /* Identifying Q1-parallelepiped.                                     */
    if (!found && nbp == (size_type(1) << n))
      if (pgt->basic_structure() == bgeot::parallelepiped_structure(n)) {
        if (!complete && k == 2 && (n == 2 || n == 3) &&
            pgt->structure() == bgeot::Q2_incomplete_structure(n)) {
          GMM_ASSERT2(alpha == scalar_type(0),
                      "Cannot use an alpha parameter in incomplete Q2"
                      " elements.");
          name << "FEM_Q2_INCOMPLETE" << suffix << "(" << n << ")";
        } else
          name << "FEM_QK" << suffix << "(" << n << "," << k << arg << ")";
        found = true;
      }
 
    /* Identifying Q1-prisms.                                             */
    if (!found && nbp == 2 * n)
      if (pgt->basic_structure() == bgeot::prism_P1_structure(n)) {
        if (!complete && k == 2 && n == 3 &&
            pgt->structure() == bgeot::prism_incomplete_P2_structure()) {
          GMM_ASSERT2(alpha == scalar_type(0),
                      "Cannot use an alpha parameter in incomplete prism"
                      " elements.");
          name << "FEM_PRISM_INCOMPLETE_P2" << suffix;
        } else
          name << "FEM_PRISM_PK" << suffix << "(" << n << "," << k << arg << ")";
        found = true;
      }
 
    /* Identifying pyramids.                                              */
    if (!found && nbp == 5)
      if (pgt->basic_structure() == bgeot::pyramid_QK_structure(1)) {
        if (!complete && k == 2 &&
            pgt->structure() == bgeot::pyramid_Q2_incomplete_structure()) {
          GMM_ASSERT2(alpha == scalar_type(0),
                      "Cannot use an alpha parameter in incomplete pyramid"
                      " elements.");
          name << "FEM_PYRAMID_Q2_INCOMPLETE" << suffix;
        } else
          name << "FEM_PYRAMID_QK" << suffix << "(" << k << arg << ")";
        found = true;;
      }
 
    // To be completed
 
    GMM_ASSERT1(found, "This element is not taken into account. Contact us");
    fem_last = fem_descriptor(name.str());
    pgt_last = pgt;
    k_last = k;
    complete_last = complete;
    discont_last = discont;
    alpha_last = alpha;
    return fem_last;
  }
 
  pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k,
                     bool complete) {
    return classical_fem_(pgt, k, complete);
  }
 
  pfem classical_discontinuous_fem(bgeot::pgeometric_trans pgt, short_type k,
                                   scalar_type alpha, bool complete) {
    return classical_fem_(pgt, k, complete, true, alpha);
  }
 
  /* ******************************************************************** */
  /*    Naming system                                                     */
  /* ******************************************************************** */
 
  pfem structured_composite_fem_method(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  pfem PK_composite_hierarch_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  pfem PK_composite_full_hierarch_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  pfem HCT_triangle_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  pfem reduced_HCT_triangle_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  pfem reduced_quadc1p3_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  pfem quadc1p3_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  pfem P1bubbletriangle_fem(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  pfem hho_method(fem_param_list &params,
        std::vector<dal::pstatic_stored_object> &dependencies);
  
  struct fem_naming_system : public dal::naming_system<virtual_fem> {
    fem_naming_system() : dal::naming_system<virtual_fem>("FEM") {
      add_suffix("HERMITE", Hermite_fem);
      add_suffix("ARGYRIS", triangle_Argyris_fem);
      add_suffix("MORLEY", element_Morley_fem);
      add_suffix("PK", PK_fem);
      add_suffix("QK", QK_fem);
      add_suffix("QK_DISCONTINUOUS", QK_discontinuous_fem);
      add_suffix("PRISM_PK", prism_PK_fem);
      add_suffix("PK_PRISM", prism_PK_fem); // for backwards compatibility
      add_suffix("PK_DISCONTINUOUS", PK_discontinuous_fem);
      add_suffix("PRISM_PK_DISCONTINUOUS", prism_PK_discontinuous_fem);
      add_suffix("PK_PRISM_DISCONTINUOUS", prism_PK_discontinuous_fem); // for backwards compatibility
      add_suffix("SIMPLEX_IPK", simplex_IPK_fem);
      add_suffix("PRISM_IPK", prism_IPK_fem);
      add_suffix("QUAD_IPK", parallelepiped_IPK_fem);
      add_suffix("SIMPLEX_CIPK", conn_int_PK_fem);
      add_suffix("PK_WITH_CUBIC_BUBBLE", PK_with_cubic_bubble);
      add_suffix("PRODUCT", product_fem);
      add_suffix("P1_NONCONFORMING", P1_nonconforming_fem);
      add_suffix("P1_BUBBLE_FACE", P1_with_bubble_on_a_face);
      add_suffix("P1_BUBBLE_FACE_LAG", P1_with_bubble_on_a_face_lagrange);
      add_suffix("P1_PIECEWISE_LINEAR_BUBBLE", P1bubbletriangle_fem);
      add_suffix("GEN_HIERARCHICAL", gen_hierarchical_fem);
      add_suffix("PK_HIERARCHICAL", PK_hierarch_fem);
      add_suffix("QK_HIERARCHICAL", QK_hierarch_fem);
      add_suffix("PRISM_PK_HIERARCHICAL", prism_PK_hierarch_fem);
      add_suffix("PK_PRISM_HIERARCHICAL", prism_PK_hierarch_fem); // for backwards compatibility
      add_suffix("STRUCTURED_COMPOSITE", structured_composite_fem_method);
      add_suffix("PK_HIERARCHICAL_COMPOSITE", PK_composite_hierarch_fem);
      add_suffix("PK_FULL_HIERARCHICAL_COMPOSITE",
                 PK_composite_full_hierarch_fem);
      add_suffix("PK_GAUSSLOBATTO1D", PK_GL_fem);
      add_suffix("QUAD_CIPK", CIPK_SQUARE);
      add_suffix("Q2_INCOMPLETE", Q2_incomplete_fem);
      add_suffix("Q2_INCOMPLETE_DISCONTINUOUS", Q2_incomplete_discontinuous_fem);
      add_suffix("HCT_TRIANGLE", HCT_triangle_fem);
      add_suffix("REDUCED_HCT_TRIANGLE", reduced_HCT_triangle_fem);
      add_suffix("QUADC1_COMPOSITE", quadc1p3_fem);
      add_suffix("REDUCED_QUADC1_COMPOSITE", reduced_quadc1p3_fem);
      add_suffix("HHO", hho_method);
      add_suffix("RT0", P1_RT0);
      add_suffix("RTK", RTk);
      add_suffix("BDMK", BDMk);
      add_suffix("RT0Q", P1_RT0Q);
      add_suffix("NEDELEC", P1_nedelec);
      add_suffix("PYRAMID_QK", pyramid_QK_fem);
      add_suffix("PYRAMID_QK_DISCONTINUOUS", pyramid_QK_disc_fem);
      add_suffix("PYRAMID_LAGRANGE", pyramid_QK_fem); // for backwards compatibility
      add_suffix("PYRAMID_DISCONTINUOUS_LAGRANGE", pyramid_QK_disc_fem); // for backwards compatibility
      add_suffix("PYRAMID_Q2_INCOMPLETE", pyramid_Q2_incomplete_fem);
      add_suffix("PYRAMID_Q2_INCOMPLETE_DISCONTINUOUS",
                 pyramid_Q2_incomplete_disc_fem);
      add_suffix("PRISM_INCOMPLETE_P2", prism_incomplete_P2_fem);
      add_suffix("PRISM_INCOMPLETE_P2_DISCONTINUOUS",
                 prism_incomplete_P2_disc_fem);
    }
  };
 
  // get a fem descriptor from a string name of a fem.
  pfem fem_descriptor(const std::string &name) {
    size_type i = 0;
    pfem  pf = dal::singleton<fem_naming_system>::instance().method(name, i);
    const_cast<virtual_fem &>(*pf).debug_name()
      = dal::singleton<fem_naming_system>::instance().shorter_name_of_method(pf);
    return pf;
  }
 
  // get the string name of a fem descriptor.
  std::string name_of_fem(pfem p) {
 
    auto &instance = dal::singleton<fem_naming_system>::instance();
    auto *p_torus = dynamic_cast<const torus_fem *>(p.get());
    if (p_torus) return instance.shorter_name_of_method(p_torus->get_original_pfem());
    return instance.shorter_name_of_method(p);
 
    return dal::singleton<fem_naming_system>::instance().
      shorter_name_of_method(p);
  }
 
  // allows the add of a fem.
  void add_fem_name(std::string name,
                    dal::naming_system<virtual_fem>::pfunction f) {
    dal::singleton<fem_naming_system>::instance().add_suffix(name, f);
  }
 
  /* ******************************************************************** */
  /*    Aliases functions                                                 */
  /* ******************************************************************** */
 
  pfem PK_fem(size_type n, short_type k) {
    THREAD_SAFE_STATIC pfem pf = nullptr;
    THREAD_SAFE_STATIC size_type d = size_type(-2);
    THREAD_SAFE_STATIC short_type r = short_type(-2);
    if (d != n || r != k) {
      std::stringstream name;
      name << "FEM_PK(" << n << "," << k << ")";
      pf = fem_descriptor(name.str());
      d = n; r = k;
    }
    return pf;
  }
 
  pfem QK_fem(size_type n, short_type k) {
    THREAD_SAFE_STATIC pfem pf = nullptr;
    THREAD_SAFE_STATIC size_type d = size_type(-2);
    THREAD_SAFE_STATIC short_type r = short_type(-2);
    if (d != n || r != k) {
      std::stringstream name;
      name << "FEM_QK(" << n << "," << k << ")";
      pf = fem_descriptor(name.str());
      d = n; r = k;
    }
    return pf;
  }
 
  pfem prism_PK_fem(size_type n, short_type k) {
    THREAD_SAFE_STATIC pfem pf = nullptr;
    THREAD_SAFE_STATIC size_type d = size_type(-2);
    THREAD_SAFE_STATIC short_type r = short_type(-2);
    if (d != n || r != k) {
      std::stringstream name;
      name << "FEM_PRISM_PK(" << n << "," << k << ")";
      pf = fem_descriptor(name.str());
      d = n; r = k;
    }
    return pf;
  }
 
 
  /* ********************************************************************* */
  /*       Precomputation on fem.                                          */
  /* ********************************************************************* */
 
  DAL_DOUBLE_KEY(pre_fem_key_, pfem, bgeot::pstored_point_tab);
 
  fem_precomp_::fem_precomp_(const pfem pff, const bgeot::pstored_point_tab ps) :
    pf(pff), pspt(ps) {
    DAL_STORED_OBJECT_DEBUG_CREATED(this, "Fem_precomp");
    for (const auto &p : *pspt)
      GMM_ASSERT1(p.size() == pf->dim(), "dimensions mismatch");
  }
 
  void fem_precomp_::init_val() const {
    c.resize(pspt->size());
    for (size_type i = 0; i < pspt->size(); ++i)
      pf->base_value((*pspt)[i], c[i]);
  }
 
  void fem_precomp_::init_grad() const {
    pc.resize(pspt->size());
    for (size_type i = 0; i < pspt->size(); ++i)
      pf->grad_base_value((*pspt)[i], pc[i]);
  }
 
  void fem_precomp_::init_hess() const {
    hpc.resize(pspt->size());
    for (size_type i = 0; i < pspt->size(); ++i)
      pf->hess_base_value((*pspt)[i], hpc[i]);
  }
 
  pfem_precomp fem_precomp(pfem pf, bgeot::pstored_point_tab pspt,
                           dal::pstatic_stored_object dep) {
    dal::pstatic_stored_object_key pk = std::make_shared<pre_fem_key_>(pf,pspt);
    dal::pstatic_stored_object o = dal::search_stored_object(pk);
    if (o) return std::dynamic_pointer_cast<const fem_precomp_>(o);
    pfem_precomp p = std::make_shared<fem_precomp_>(pf, pspt);
    dal::add_stored_object(pk, p, pspt, dal::AUTODELETE_STATIC_OBJECT);
    if (dal::exists_stored_object(pf)) dal::add_dependency(p, pf);
    if (dep) dal::add_dependency(p, dep);
    return p;
  }
 
  void fem_precomp_pool::clear() {
    for (const pfem_precomp &p : precomps)
      dal::del_stored_object(p, true);
    precomps.clear();
  }
 
 
}  /* end of namespace getfem.                                            */