diff --git a/src/fe_engine/fe_engine_template_tmpl.hh b/src/fe_engine/fe_engine_template_tmpl.hh
index 8d5495c1e..7d520b481 100644
--- a/src/fe_engine/fe_engine_template_tmpl.hh
+++ b/src/fe_engine/fe_engine_template_tmpl.hh
@@ -1,1402 +1,1403 @@
 /**
  * @file   fe_engine_template_tmpl.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Mauro Corrado <mauro.corrado@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Tue Feb 15 2011
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Template implementation of FEEngineTemplate
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 #include "dof_manager.hh"
 #include "fe_engine_template.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::FEEngineTemplate(
     Mesh & mesh, UInt spatial_dimension, ID id, MemoryID memory_id)
     : FEEngine(mesh, spatial_dimension, id, memory_id),
-      integrator(mesh, id, memory_id), shape_functions(mesh, id, memory_id) {}
+      integrator(mesh, spatial_dimension, id, memory_id),
+      shape_functions(mesh, spatial_dimension, id, memory_id) {}
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::~FEEngineTemplate() =
     default;
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GradientOnIntegrationPointsHelper {
       template <class S>
       static void call(const S &, Mesh &, const Array<Real> &, Array<Real> &,
                        const UInt, const ElementType &, const GhostType &,
                        const Array<UInt> &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_GRADIENT(type)                                                 \
   if (element_dimension == ElementClass<type>::getSpatialDimension())          \
     shape_functions.template gradientOnIntegrationPoints<type>(                \
         u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER(kind)          \
   template <> struct GradientOnIntegrationPointsHelper<kind> {                 \
     template <class S>                                                         \
     static void call(const S & shape_functions, Mesh & mesh,                   \
                      const Array<Real> & u, Array<Real> & nablauq,             \
                      const UInt nb_degree_of_freedom,                          \
                      const ElementType & type, const GhostType & ghost_type,   \
                      const Array<UInt> & filter_elements) {                    \
       UInt element_dimension = mesh.getSpatialDimension(type);                 \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_GRADIENT, kind);                \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER,
         AKANTU_FE_ENGINE_LIST_GRADIENT_ON_INTEGRATION_POINTS)
 
 #undef AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER
 #undef COMPUTE_GRADIENT
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     gradientOnIntegrationPoints(const Array<Real> & u, Array<Real> & nablauq,
                                 const UInt nb_degree_of_freedom,
                                 const ElementType & type,
                                 const GhostType & ghost_type,
                                 const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter)
     nb_element = filter_elements.size();
   UInt nb_points =
       shape_functions.getIntegrationPoints(type, ghost_type).cols();
 
 #ifndef AKANTU_NDEBUG
 
   UInt element_dimension = mesh.getSpatialDimension(type);
 
   AKANTU_DEBUG_ASSERT(u.size() == mesh.getNbNodes(),
                       "The vector u(" << u.getID()
                                       << ") has not the good size.");
   AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom,
                       "The vector u("
                           << u.getID()
                           << ") has not the good number of component.");
 
   AKANTU_DEBUG_ASSERT(
       nablauq.getNbComponent() == nb_degree_of_freedom * element_dimension,
       "The vector nablauq(" << nablauq.getID()
                             << ") has not the good number of component.");
 
 // AKANTU_DEBUG_ASSERT(nablauq.size() == nb_element * nb_points,
 //                  "The vector nablauq(" << nablauq.getID()
 //                  << ") has not the good size.");
 #endif
 
   nablauq.resize(nb_element * nb_points);
 
   fe_engine::details::GradientOnIntegrationPointsHelper<kind>::call(
       shape_functions, mesh, u, nablauq, nb_degree_of_freedom, type, ghost_type,
       filter_elements);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::initShapeFunctions(
     const GhostType & ghost_type) {
   initShapeFunctions(mesh.getNodes(), ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::initShapeFunctions(
     const Array<Real> & nodes, const GhostType & ghost_type) {
   AKANTU_DEBUG_IN();
 
   for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) {
     integrator.initIntegrator(nodes, type, ghost_type);
     const auto & control_points = getIntegrationPoints(type, ghost_type);
     shape_functions.initShapeFunctions(nodes, control_points, type, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct IntegrateHelper {};
 
 #define INTEGRATE(type)                                                        \
   integrator.template integrate<type>(f, intf, nb_degree_of_freedom,           \
                                       ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_INTEGRATE_HELPER(kind)                               \
   template <> struct IntegrateHelper<kind> {                                   \
     template <class I>                                                         \
     static void call(const I & integrator, const Array<Real> & f,              \
                      Array<Real> & intf, UInt nb_degree_of_freedom,            \
                      const ElementType & type, const GhostType & ghost_type,   \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_HELPER)
 
 #undef AKANTU_SPECIALIZE_INTEGRATE_HELPER
 #undef INTEGRATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & f, Array<Real> & intf, UInt nb_degree_of_freedom,
     const ElementType & type, const GhostType & ghost_type,
     const Array<UInt> & filter_elements) const {
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter)
     nb_element = filter_elements.size();
 #ifndef AKANTU_NDEBUG
 
   UInt nb_quadrature_points = getNbIntegrationPoints(type);
 
   AKANTU_DEBUG_ASSERT(f.size() == nb_element * nb_quadrature_points,
                       "The vector f(" << f.getID() << " size " << f.size()
                                       << ") has not the good size ("
                                       << nb_element << ").");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
                       "The vector intf("
                           << intf.getID()
                           << ") has not the good number of component.");
 #endif
 
   intf.resize(nb_element);
 
   fe_engine::details::IntegrateHelper<kind>::call(integrator, f, intf,
                                                   nb_degree_of_freedom, type,
                                                   ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct IntegrateScalarHelper {};
 
 #define INTEGRATE(type)                                                        \
   integral =                                                                   \
       integrator.template integrate<type>(f, ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER(kind)                        \
   template <> struct IntegrateScalarHelper<kind> {                             \
     template <class I>                                                         \
     static Real call(const I & integrator, const Array<Real> & f,              \
                      const ElementType & type, const GhostType & ghost_type,   \
                      const Array<UInt> & filter_elements) {                    \
       Real integral = 0.;                                                      \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
       return integral;                                                         \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER)
 
 #undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER
 #undef INTEGRATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 Real FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & f, const ElementType & type,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
 #ifndef AKANTU_NDEBUG
   //   std::stringstream sstr; sstr << ghost_type;
   //   AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
   //                  "The vector " << nablauq.getID() << " is not taged " <<
   //                  ghost_type);
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter)
     nb_element = filter_elements.size();
 
   UInt nb_quadrature_points = getNbIntegrationPoints(type, ghost_type);
 
   AKANTU_DEBUG_ASSERT(
       f.size() == nb_element * nb_quadrature_points,
       "The vector f(" << f.getID() << ") has not the good size. (" << f.size()
                       << "!=" << nb_quadrature_points * nb_element << ")");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == 1,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
 #endif
 
   Real integral = fe_engine::details::IntegrateScalarHelper<kind>::call(
       integrator, f, type, ghost_type, filter_elements);
   AKANTU_DEBUG_OUT();
   return integral;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct IntegrateScalarOnOneElementHelper {};
 
 #define INTEGRATE(type)                                                        \
   res = integrator.template integrate<type>(f, index, ghost_type);
 
 #define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER(kind)         \
   template <> struct IntegrateScalarOnOneElementHelper<kind> {                 \
     template <class I>                                                         \
     static Real call(const I & integrator, const Vector<Real> & f,             \
                      const ElementType & type, UInt index,                     \
                      const GhostType & ghost_type) {                           \
       Real res = 0.;                                                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
       return res;                                                              \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(
         AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER)
 
 #undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER
 #undef INTEGRATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 Real FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
     const Vector<Real> & f, const ElementType & type, UInt index,
     const GhostType & ghost_type) const {
 
   Real res = fe_engine::details::IntegrateScalarOnOneElementHelper<kind>::call(
       integrator, f, type, index, ghost_type);
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct IntegrateOnIntegrationPointsHelper {};
 
 #define INTEGRATE(type)                                                        \
   integrator.template integrateOnIntegrationPoints<type>(                      \
       f, intf, nb_degree_of_freedom, ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER(kind)         \
   template <> struct IntegrateOnIntegrationPointsHelper<kind> {                \
     template <class I>                                                         \
     static void call(const I & integrator, const Array<Real> & f,              \
                      Array<Real> & intf, UInt nb_degree_of_freedom,            \
                      const ElementType & type, const GhostType & ghost_type,   \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(
         AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER)
 
 #undef AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER
 #undef INTEGRATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     integrateOnIntegrationPoints(const Array<Real> & f, Array<Real> & intf,
                                  UInt nb_degree_of_freedom,
                                  const ElementType & type,
                                  const GhostType & ghost_type,
                                  const Array<UInt> & filter_elements) const {
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter)
     nb_element = filter_elements.size();
   UInt nb_quadrature_points = getNbIntegrationPoints(type);
 #ifndef AKANTU_NDEBUG
   //   std::stringstream sstr; sstr << ghost_type;
   //   AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
   //                  "The vector " << nablauq.getID() << " is not taged " <<
   //                  ghost_type);
 
   AKANTU_DEBUG_ASSERT(f.size() == nb_element * nb_quadrature_points,
                       "The vector f(" << f.getID() << " size " << f.size()
                                       << ") has not the good size ("
                                       << nb_element << ").");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
                       "The vector intf("
                           << intf.getID()
                           << ") has not the good number of component.");
 #endif
 
   intf.resize(nb_element * nb_quadrature_points);
   fe_engine::details::IntegrateOnIntegrationPointsHelper<kind>::call(
       integrator, f, intf, nb_degree_of_freedom, type, ghost_type,
       filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct InterpolateOnIntegrationPointsHelper {
       template <class S>
       static void call(const S &, const Array<Real> &, Array<Real> &,
                        const UInt, const ElementType &, const GhostType &,
                        const Array<UInt> &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define INTERPOLATE(type)                                                      \
   shape_functions.template interpolateOnIntegrationPoints<type>(               \
       u, uq, nb_degree_of_freedom, ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER(kind)       \
   template <> struct InterpolateOnIntegrationPointsHelper<kind> {              \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & u,         \
                      Array<Real> & uq, const UInt nb_degree_of_freedom,        \
                      const ElementType & type, const GhostType & ghost_type,   \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind);                     \
     }                                                                          \
   };
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER,
         AKANTU_FE_ENGINE_LIST_INTERPOLATE_ON_INTEGRATION_POINTS)
 
 #undef AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER
 #undef INTERPOLATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateOnIntegrationPoints(const Array<Real> & u, Array<Real> & uq,
                                    const UInt nb_degree_of_freedom,
                                    const ElementType & type,
                                    const GhostType & ghost_type,
                                    const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_points =
       shape_functions.getIntegrationPoints(type, ghost_type).cols();
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter)
     nb_element = filter_elements.size();
 #ifndef AKANTU_NDEBUG
 
   AKANTU_DEBUG_ASSERT(u.size() == mesh.getNbNodes(),
                       "The vector u(" << u.getID()
                                       << ") has not the good size.");
   AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom,
                       "The vector u("
                           << u.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(uq.getNbComponent() == nb_degree_of_freedom,
                       "The vector uq("
                           << uq.getID()
                           << ") has not the good number of component.");
 #endif
 
   uq.resize(nb_element * nb_points);
 
   fe_engine::details::InterpolateOnIntegrationPointsHelper<kind>::call(
       shape_functions, u, uq, nb_degree_of_freedom, type, ghost_type,
       filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateOnIntegrationPoints(
         const Array<Real> & u, ElementTypeMapArray<Real> & uq,
         const ElementTypeMapArray<UInt> * filter_elements) const {
   AKANTU_DEBUG_IN();
 
   const Array<UInt> * filter = nullptr;
 
   for (auto ghost_type : ghost_types) {
     for (auto & type : uq.elementTypes(_all_dimensions, ghost_type, kind)) {
       UInt nb_quad_per_element = getNbIntegrationPoints(type, ghost_type);
 
       UInt nb_element = 0;
 
       if (filter_elements) {
         filter = &((*filter_elements)(type, ghost_type));
         nb_element = filter->size();
       } else {
         filter = &empty_filter;
         nb_element = mesh.getNbElement(type, ghost_type);
       }
 
       UInt nb_tot_quad = nb_quad_per_element * nb_element;
 
       Array<Real> & quad = uq(type, ghost_type);
       quad.resize(nb_tot_quad);
 
       interpolateOnIntegrationPoints(u, quad, quad.getNbComponent(), type,
                                      ghost_type, *filter);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeBtDHelper {};
 
 #define COMPUTE_BTD(type)                                                      \
   shape_functions.template computeBtD<type>(Ds, BtDs, ghost_type,              \
                                             filter_elements);
 
 #define AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER(kind)                             \
   template <> struct ComputeBtDHelper<kind> {                                  \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & Ds,        \
                      Array<Real> & BtDs, const ElementType & type,             \
                      const GhostType & ghost_type,                             \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_BTD, kind);                     \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER
 #undef COMPUTE_BTD
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeBtD(
     const Array<Real> & Ds, Array<Real> & BtDs, const ElementType & type,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
   fe_engine::details::ComputeBtDHelper<kind>::call(
       shape_functions, Ds, BtDs, type, ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeBtDBHelper {};
 
 #define COMPUTE_BTDB(type)                                                     \
   shape_functions.template computeBtDB<type>(Ds, BtDBs, order_d, ghost_type,   \
                                              filter_elements);
 
 #define AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER(kind)                            \
   template <> struct ComputeBtDBHelper<kind> {                                 \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & Ds,        \
                      Array<Real> & BtDBs, UInt order_d,                        \
                      const ElementType & type, const GhostType & ghost_type,   \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_BTDB, kind);                    \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER
 #undef COMPUTE_BTDB
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeBtDB(
     const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
     const ElementType & type, const GhostType & ghost_type,
     const Array<UInt> & filter_elements) const {
   fe_engine::details::ComputeBtDBHelper<kind>::call(
       shape_functions, Ds, BtDBs, order_d, type, ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeNtbHelper {};
 
 #define COMPUTE_Ntb(type)                                                      \
   shape_functions.template computeNtb<type>(bs, Ntbs, ghost_type,              \
                                             filter_elements);
 
 #define AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER(kind)                             \
   template <> struct ComputeNtbHelper<kind> {                                  \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & bs,        \
                      Array<Real> & Ntbs, const ElementType & type,             \
                      const GhostType & ghost_type,                             \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_Ntb, kind);                     \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER
 #undef COMPUTE_Ntb
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeNtb(
     const Array<Real> & bs, Array<Real> & Ntbs, const ElementType & type,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
   fe_engine::details::ComputeNtbHelper<kind>::call(
       shape_functions, bs, Ntbs, type, ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeIntegrationPointsCoordinates(
         ElementTypeMapArray<Real> & quadrature_points_coordinates,
         const ElementTypeMapArray<UInt> * filter_elements) const {
 
   const Array<Real> & nodes_coordinates = mesh.getNodes();
 
   interpolateOnIntegrationPoints(
       nodes_coordinates, quadrature_points_coordinates, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeIntegrationPointsCoordinates(
         Array<Real> & quadrature_points_coordinates, const ElementType & type,
         const GhostType & ghost_type,
         const Array<UInt> & filter_elements) const {
 
   const Array<Real> & nodes_coordinates = mesh.getNodes();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
 
   interpolateOnIntegrationPoints(
       nodes_coordinates, quadrature_points_coordinates, spatial_dimension, type,
       ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     initElementalFieldInterpolationFromIntegrationPoints(
         const ElementTypeMapArray<Real> & interpolation_points_coordinates,
         ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
         ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
         const ElementTypeMapArray<UInt> * element_filter) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = this->mesh.getSpatialDimension();
 
   ElementTypeMapArray<Real> quadrature_points_coordinates(
       "quadrature_points_coordinates_for_interpolation", getID(),
       getMemoryID());
 
   quadrature_points_coordinates.initialize(*this,
                                            _nb_component = spatial_dimension);
 
   computeIntegrationPointsCoordinates(quadrature_points_coordinates,
                                       element_filter);
   shape_functions.initElementalFieldInterpolationFromIntegrationPoints(
       interpolation_points_coordinates,
       interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, quadrature_points_coordinates,
       element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateElementalFieldFromIntegrationPoints(
         const ElementTypeMapArray<Real> & field,
         const ElementTypeMapArray<Real> & interpolation_points_coordinates,
         ElementTypeMapArray<Real> & result, const GhostType ghost_type,
         const ElementTypeMapArray<UInt> * element_filter) const {
 
   ElementTypeMapArray<Real> interpolation_points_coordinates_matrices(
       "interpolation_points_coordinates_matrices", id, memory_id);
   ElementTypeMapArray<Real> quad_points_coordinates_inv_matrices(
       "quad_points_coordinates_inv_matrices", id, memory_id);
 
   initElementalFieldInterpolationFromIntegrationPoints(
       interpolation_points_coordinates,
       interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, element_filter);
 
   interpolateElementalFieldFromIntegrationPoints(
       field, interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, result, ghost_type, element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateElementalFieldFromIntegrationPoints(
         const ElementTypeMapArray<Real> & field,
         const ElementTypeMapArray<Real> &
             interpolation_points_coordinates_matrices,
         const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
         ElementTypeMapArray<Real> & result, const GhostType ghost_type,
         const ElementTypeMapArray<UInt> * element_filter) const {
 
   shape_functions.interpolateElementalFieldFromIntegrationPoints(
       field, interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, result, ghost_type, element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct InterpolateHelper {
       template <class S>
       static void call(const S &, const Vector<Real> &, UInt,
                        const Matrix<Real> &, Vector<Real> &,
                        const ElementType &, const GhostType &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define INTERPOLATE(type)                                                      \
   shape_functions.template interpolate<type>(                                  \
       real_coords, element, nodal_values, interpolated, ghost_type);
 
 #define AKANTU_SPECIALIZE_INTERPOLATE_HELPER(kind)                             \
   template <> struct InterpolateHelper<kind> {                                 \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, UInt element,           \
                      const Matrix<Real> & nodal_values,                        \
                      Vector<Real> & interpolated, const ElementType & type,    \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind);                     \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INTERPOLATE_HELPER,
                                AKANTU_FE_ENGINE_LIST_INTERPOLATE)
 
 #undef AKANTU_SPECIALIZE_INTERPOLATE_HELPER
 #undef INTERPOLATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::interpolate(
     const Vector<Real> & real_coords, const Matrix<Real> & nodal_values,
     Vector<Real> & interpolated, const Element & element) const {
 
   AKANTU_DEBUG_IN();
 
   fe_engine::details::InterpolateHelper<kind>::call(
       shape_functions, real_coords, element.element, nodal_values, interpolated,
       element.type, element.ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const GhostType & ghost_type) {
   AKANTU_DEBUG_IN();
 
   computeNormalsOnIntegrationPoints(mesh.getNodes(), ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                       const GhostType & ghost_type) {
   AKANTU_DEBUG_IN();
 
   //  Real * coord = mesh.getNodes().storage();
   UInt spatial_dimension = mesh.getSpatialDimension();
 
   // allocate the normal arrays
   normals_on_integration_points.initialize(
       *this, _nb_component = spatial_dimension,
       _spatial_dimension = element_dimension, _ghost_type = ghost_type,
       _element_kind = kind);
 
   // loop over the type to build the normals
   for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) {
     auto & normals_on_quad = normals_on_integration_points(type, ghost_type);
     computeNormalsOnIntegrationPoints(field, normals_on_quad, type, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeNormalsOnIntegrationPoints {
       template <template <ElementKind, class> class I,
                 template <ElementKind> class S, ElementKind k, class IOF>
       static void call(const FEEngineTemplate<I, S, k, IOF> &,
                        const Array<Real> &, Array<Real> &, const ElementType &,
                        const GhostType &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_NORMALS_ON_INTEGRATION_POINTS(type)                            \
   fem.template computeNormalsOnIntegrationPoints<type>(field, normal,          \
                                                        ghost_type);
 
 #define AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS(kind)          \
   template <> struct ComputeNormalsOnIntegrationPoints<kind> {                 \
     template <template <ElementKind, class> class I,                           \
               template <ElementKind> class S, ElementKind k, class IOF>        \
     static void call(const FEEngineTemplate<I, S, k, IOF> & fem,               \
                      const Array<Real> & field, Array<Real> & normal,          \
                      const ElementType & type, const GhostType & ghost_type) { \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_INTEGRATION_POINTS,  \
                                        kind);                                  \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS,
         AKANTU_FE_ENGINE_LIST_COMPUTE_NORMALS_ON_INTEGRATION_POINTS)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS
 #undef COMPUTE_NORMALS_ON_INTEGRATION_POINTS
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                       Array<Real> & normal,
                                       const ElementType & type,
                                       const GhostType & ghost_type) const {
   fe_engine::details::ComputeNormalsOnIntegrationPoints<kind>::call(
       *this, field, normal, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                       Array<Real> & normal,
                                       const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   if (type == _point_1) {
     computeNormalsOnIntegrationPointsPoint1(field, normal, ghost_type);
     return;
   }
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_points = getNbIntegrationPoints(type, ghost_type);
 
   UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
   normal.resize(nb_element * nb_points);
   Array<Real>::matrix_iterator normals_on_quad =
       normal.begin_reinterpret(spatial_dimension, nb_points, nb_element);
   Array<Real> f_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, field, f_el, type, ghost_type);
 
   const Matrix<Real> & quads =
       integrator.template getIntegrationPoints<type>(ghost_type);
 
   Array<Real>::matrix_iterator f_it =
       f_el.begin(spatial_dimension, nb_nodes_per_element);
 
   for (UInt elem = 0; elem < nb_element; ++elem) {
     ElementClass<type>::computeNormalsOnNaturalCoordinates(quads, *f_it,
                                                            *normals_on_quad);
     ++normals_on_quad;
     ++f_it;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 template <ElementKind kind> struct InverseMapHelper {
   template <class S>
   static void
   call(const S & /*shape_functions*/, const Vector<Real> & /*real_coords*/,
        UInt /*element*/, const ElementType & /*type*/,
        Vector<Real> & /*natural_coords*/, const GhostType & /*ghost_type*/) {
     AKANTU_TO_IMPLEMENT();
   }
 };
 
 #define INVERSE_MAP(type)                                                      \
   shape_functions.template inverseMap<type>(real_coords, element,              \
                                             natural_coords, ghost_type);
 
 #define AKANTU_SPECIALIZE_INVERSE_MAP_HELPER(kind)                             \
   template <> struct InverseMapHelper<kind> {                                  \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, UInt element,           \
                      const ElementType & type, Vector<Real> & natural_coords,  \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INVERSE_MAP, kind);                     \
     }                                                                          \
   };
 
 AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INVERSE_MAP_HELPER,
                            AKANTU_FE_ENGINE_LIST_INVERSE_MAP)
 
 #undef AKANTU_SPECIALIZE_INVERSE_MAP_HELPER
 #undef INVERSE_MAP
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::inverseMap(
     const Vector<Real> & real_coords, UInt element, const ElementType & type,
     Vector<Real> & natural_coords, const GhostType & ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   InverseMapHelper<kind>::call(shape_functions, real_coords, element, type,
                                natural_coords, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ContainsHelper {
       template <class S>
       static void call(const S &, const Vector<Real> &, UInt,
                        const ElementType &, const GhostType &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define CONTAINS(type)                                                         \
   contain = shape_functions.template contains<type>(real_coords, element,      \
                                                     ghost_type);
 
 #define AKANTU_SPECIALIZE_CONTAINS_HELPER(kind)                                \
   template <> struct ContainsHelper<kind> {                                    \
     template <template <ElementKind> class S, ElementKind k>                   \
     static bool call(const S<k> & shape_functions,                             \
                      const Vector<Real> & real_coords, UInt element,           \
                      const ElementType & type, const GhostType & ghost_type) { \
       bool contain = false;                                                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(CONTAINS, kind);                        \
       return contain;                                                          \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_CONTAINS_HELPER,
                                AKANTU_FE_ENGINE_LIST_CONTAINS)
 
 #undef AKANTU_SPECIALIZE_CONTAINS_HELPER
 #undef CONTAINS
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline bool FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::contains(
     const Vector<Real> & real_coords, UInt element, const ElementType & type,
     const GhostType & ghost_type) const {
   return fe_engine::details::ContainsHelper<kind>::call(
       shape_functions, real_coords, element, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeShapesHelper {
       template <class S>
       static void call(const S &, const Vector<Real> &, UInt, const ElementType,
                        Vector<Real> &, const GhostType &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_SHAPES(type)                                                   \
   shape_functions.template computeShapes<type>(real_coords, element, shapes,   \
                                                ghost_type);
 
 #define AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER(kind)                          \
   template <> struct ComputeShapesHelper<kind> {                               \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, UInt element,           \
                      const ElementType type, Vector<Real> & shapes,            \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPES, kind);                  \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER,
                                AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER
 #undef COMPUTE_SHAPES
   } // namespace details
 } // namespace fe_engine
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapes(
     const Vector<Real> & real_coords, UInt element, const ElementType & type,
     Vector<Real> & shapes, const GhostType & ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   fe_engine::details::ComputeShapesHelper<kind>::call(
       shape_functions, real_coords, element, type, shapes, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeShapeDerivativesHelper {
       template <class S>
       static void call(__attribute__((unused)) const S & shape_functions,
                        __attribute__((unused)) const Vector<Real> & real_coords,
                        __attribute__((unused)) UInt element,
                        __attribute__((unused)) const ElementType type,
                        __attribute__((unused)) Matrix<Real> & shape_derivatives,
                        __attribute__((unused)) const GhostType & ghost_type) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_SHAPE_DERIVATIVES(type)                                        \
   Matrix<Real> coords_mat(real_coords.storage(), shape_derivatives.rows(), 1); \
   Tensor3<Real> shapesd_tensor(shape_derivatives.storage(),                    \
                                shape_derivatives.rows(),                       \
                                shape_derivatives.cols(), 1);                   \
   shape_functions.template computeShapeDerivatives<type>(                      \
       coords_mat, element, shapesd_tensor, ghost_type);
 
 #define AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER(kind)               \
   template <> struct ComputeShapeDerivativesHelper<kind> {                     \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, UInt element,           \
                      const ElementType type, Matrix<Real> & shape_derivatives, \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPE_DERIVATIVES, kind);       \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER,
         AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES_DERIVATIVES)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER
 #undef COMPUTE_SHAPE_DERIVATIVES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapeDerivatives(
     const Vector<Real> & real_coords, UInt element, const ElementType & type,
     Matrix<Real> & shape_derivatives, const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   fe_engine::details::ComputeShapeDerivativesHelper<kind>::call(
       shape_functions, real_coords, element, type, shape_derivatives,
       ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetNbIntegrationPointsHelper {};
 
 #define GET_NB_INTEGRATION_POINTS(type)                                        \
   nb_quad_points = integrator.template getNbIntegrationPoints<type>(ghost_type);
 
 #define AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER(kind)               \
   template <> struct GetNbIntegrationPointsHelper<kind> {                      \
     template <template <ElementKind, class> class I, ElementKind k, class IOF> \
     static UInt call(const I<k, IOF> & integrator, const ElementType type,     \
                      const GhostType & ghost_type) {                           \
       UInt nb_quad_points = 0;                                                 \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_NB_INTEGRATION_POINTS, kind);       \
       return nb_quad_points;                                                   \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER
 #undef GET_NB_INTEGRATION
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline UInt
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getNbIntegrationPoints(
     const ElementType & type, const GhostType & ghost_type) const {
   return fe_engine::details::GetNbIntegrationPointsHelper<kind>::call(
       integrator, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetShapesHelper {};
 
 #define GET_SHAPES(type) ret = &(shape_functions.getShapes(type, ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_SHAPES_HELPER(kind)                              \
   template <> struct GetShapesHelper<kind> {                                   \
     template <class S>                                                         \
     static const Array<Real> & call(const S & shape_functions,                 \
                                     const ElementType type,                    \
                                     const GhostType & ghost_type) {            \
       const Array<Real> * ret = NULL;                                          \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES, kind);                      \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_SHAPES_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_SHAPES_HELPER
 #undef GET_SHAPES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Array<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getShapes(
     const ElementType & type, const GhostType & ghost_type,
     __attribute__((unused)) UInt id) const {
   return fe_engine::details::GetShapesHelper<kind>::call(shape_functions, type,
                                                          ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetShapesDerivativesHelper {
       template <template <ElementKind> class S, ElementKind k>
       static const Array<Real> & call(const S<k> &, const ElementType &,
                                       const GhostType &, UInt) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define GET_SHAPES_DERIVATIVES(type)                                           \
   ret = &(shape_functions.getShapesDerivatives(type, ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER(kind)                  \
   template <> struct GetShapesDerivativesHelper<kind> {                        \
     template <template <ElementKind> class S, ElementKind k>                   \
     static const Array<Real> &                                                 \
     call(const S<k> & shape_functions, const ElementType type,                 \
          const GhostType & ghost_type, __attribute__((unused)) UInt id) {      \
       const Array<Real> * ret = NULL;                                          \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES_DERIVATIVES, kind);          \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER,
                                AKANTU_FE_ENGINE_LIST_GET_SHAPES_DERIVATIVES)
 
 #undef AKANTU_SPECIALIZE_GET_SHAPE_DERIVATIVES_HELPER
 #undef GET_SHAPES_DERIVATIVES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Array<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getShapesDerivatives(
     const ElementType & type, const GhostType & ghost_type,
     __attribute__((unused)) UInt id) const {
   return fe_engine::details::GetShapesDerivativesHelper<kind>::call(
       shape_functions, type, ghost_type, id);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetIntegrationPointsHelper {};
 
 #define GET_INTEGRATION_POINTS(type)                                           \
   ret = &(integrator.template getIntegrationPoints<type>(ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER(kind)                  \
   template <> struct GetIntegrationPointsHelper<kind> {                        \
     template <template <ElementKind, class> class I, ElementKind k, class IOF> \
     static const Matrix<Real> & call(const I<k, IOF> & integrator,             \
                                      const ElementType type,                   \
                                      const GhostType & ghost_type) {           \
       const Matrix<Real> * ret = NULL;                                         \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_INTEGRATION_POINTS, kind);          \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER
 #undef GET_INTEGRATION_POINTS
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Matrix<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getIntegrationPoints(
     const ElementType & type, const GhostType & ghost_type) const {
   return fe_engine::details::GetIntegrationPointsHelper<kind>::call(
       integrator, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::printself(
     std::ostream & stream, int indent) const {
   std::string space(indent, AKANTU_INDENT);
 
   stream << space << "FEEngineTemplate [" << std::endl;
   stream << space << " + parent [" << std::endl;
   FEEngine::printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << " + shape functions [" << std::endl;
   shape_functions.printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << " + integrator [" << std::endl;
   integrator.printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsAdded(
     const Array<Element> & new_elements, const NewElementsEvent &) {
   integrator.onElementsAdded(new_elements);
   shape_functions.onElementsAdded(new_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsRemoved(
     const Array<Element> &, const ElementTypeMapArray<UInt> &,
     const RemovedElementsEvent &) {}
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsChanged(
     const Array<Element> &, const Array<Element> &,
     const ElementTypeMapArray<UInt> &, const ChangedElementsEvent &) {}
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPointsPoint1(
         const Array<Real> &, Array<Real> & normal,
         const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == 1,
                       "Mesh dimension must be 1 to compute normals on points!");
   const auto type = _point_1;
   auto spatial_dimension = mesh.getSpatialDimension();
   // UInt nb_nodes_per_element  = Mesh::getNbNodesPerElement(type);
   auto nb_points = getNbIntegrationPoints(type, ghost_type);
   const auto & connectivity = mesh.getConnectivity(type, ghost_type);
   auto nb_element = connectivity.size();
 
   normal.resize(nb_element * nb_points);
   auto normals_on_quad =
       normal.begin_reinterpret(spatial_dimension, nb_points, nb_element);
   const auto & segments = mesh.getElementToSubelement(type, ghost_type);
   const auto & coords = mesh.getNodes();
 
   const Mesh * mesh_segment;
   if (mesh.isMeshFacets())
     mesh_segment = &(mesh.getMeshParent());
   else
     mesh_segment = &mesh;
 
   for (UInt elem = 0; elem < nb_element; ++elem) {
     UInt nb_segment = segments(elem).size();
     AKANTU_DEBUG_ASSERT(
         nb_segment > 0,
         "Impossible to compute a normal on a point connected to 0 segments");
 
     Real normal_value = 1;
     if (nb_segment == 1) {
       auto point = connectivity(elem);
       const auto segment = segments(elem)[0];
       const auto & segment_connectivity =
           mesh_segment->getConnectivity(segment.type, segment.ghost_type);
       Vector<UInt> segment_points = segment_connectivity.begin(
           Mesh::getNbNodesPerElement(segment.type))[segment.element];
       Real difference;
       if (segment_points(0) == point) {
         difference = coords(elem) - coords(segment_points(1));
       } else {
         difference = coords(elem) - coords(segment_points(0));
       }
 
       normal_value = difference / std::abs(difference);
     }
 
     for (UInt n(0); n < nb_points; ++n) {
       (*normals_on_quad)(0, n) = normal_value;
     }
     ++normals_on_quad;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
diff --git a/src/fe_engine/integrator.hh b/src/fe_engine/integrator.hh
index c751cfff4..1cc92d1f0 100644
--- a/src/fe_engine/integrator.hh
+++ b/src/fe_engine/integrator.hh
@@ -1,135 +1,139 @@
 /**
  * @file   integrator.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Sun Dec 03 2017
  *
  * @brief  interface for integrator classes
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 #ifndef __AKANTU_INTEGRATOR_HH__
 #define __AKANTU_INTEGRATOR_HH__
 
 /* -------------------------------------------------------------------------- */
 #include "aka_memory.hh"
 #include "mesh.hh"
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 
 class Integrator : protected Memory {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
-  Integrator(const Mesh & mesh, const ID & id = "integrator",
-             const MemoryID & memory_id = 0)
+  Integrator(const Mesh & mesh, UInt spatial_dimension,
+             const ID & id = "integrator", const MemoryID & memory_id = 0)
       : Memory(id, memory_id), mesh(mesh),
+        _spatial_dimension(spatial_dimension),
         jacobians("jacobians", id, memory_id) {
     AKANTU_DEBUG_IN();
 
     AKANTU_DEBUG_OUT();
   };
 
   ~Integrator() override = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// empty method
   template <ElementType type>
   inline void precomputeJacobiansOnQuadraturePoints(__attribute__((unused))
                                                     GhostType ghost_type) {}
 
   /// empty method
   void integrateOnElement(const Array<Real> & /*f*/, Real * /*intf*/,
                           UInt /*nb_degree_of_freedom*/,
                           const Element & /*elem*/,
                           GhostType /*ghost_type*/) const {};
 
   /// function to print the contain of the class
   virtual void printself(std::ostream & stream, int indent = 0) const {
     std::string space;
     for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
       ;
     stream << space << "Integrator [" << std::endl;
     jacobians.printself(stream, indent + 1);
     stream << space << "]" << std::endl;
   };
 
   /* ------------------------------------------------------------------------ */
 public:
   virtual void onElementsAdded(const Array<Element> &) {}
   virtual void
   onElementsRemoved(const Array<Element> &,
                     const ElementTypeMapArray<UInt> & new_numbering) {
     jacobians.onElementsRemoved(new_numbering);
   }
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// access to the jacobians
   Array<Real> & getJacobians(const ElementType & type,
                              const GhostType & ghost_type = _not_ghost) {
     return jacobians(type, ghost_type);
   };
 
   /// access to the jacobians const
   const Array<Real> &
   getJacobians(const ElementType & type,
                const GhostType & ghost_type = _not_ghost) const {
     return jacobians(type, ghost_type);
   };
 
   AKANTU_GET_MACRO(Jacobians, jacobians, const ElementTypeMapArray<Real> &);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// mesh associated to the integrator
   const Mesh & mesh;
 
+  // spatial dimension of the elements to consider
+  UInt _spatial_dimension;
+
   /// jacobians for all elements
   ElementTypeMapArray<Real> jacobians;
 };
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 
 //#include "integrator_inline_impl.cc"
 
 /// standard output stream operator
 inline std::ostream & operator<<(std::ostream & stream,
                                  const Integrator & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #endif /* __AKANTU_INTEGRATOR_HH__ */
diff --git a/src/fe_engine/integrator_gauss.hh b/src/fe_engine/integrator_gauss.hh
index da9211479..f05a6471f 100644
--- a/src/fe_engine/integrator_gauss.hh
+++ b/src/fe_engine/integrator_gauss.hh
@@ -1,205 +1,206 @@
 /**
  * @file   integrator_gauss.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Wed Nov 08 2017
  *
  * @brief  Gauss integration facilities
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "integrator.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_INTEGRATOR_GAUSS_HH__
 #define __AKANTU_INTEGRATOR_GAUSS_HH__
 
 namespace akantu {
 namespace integrator {
   namespace details {
     template <ElementKind> struct GaussIntegratorComputeJacobiansHelper;
   } // namespace details
 } // namespace integrator
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 class IntegratorGauss : public Integrator {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
-  IntegratorGauss(const Mesh & mesh, const ID & id = "integrator_gauss",
+  IntegratorGauss(const Mesh & mesh, UInt spatial_dimension,
+                  const ID & id = "integrator_gauss",
                   const MemoryID & memory_id = 0);
 
   ~IntegratorGauss() override = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   void initIntegrator(const Array<Real> & nodes, const ElementType & type,
                       const GhostType & ghost_type);
 
   template <ElementType type>
   inline void initIntegrator(const Array<Real> & nodes,
                              const GhostType & ghost_type);
 
   /// integrate f on the element "elem" of type "type"
   template <ElementType type>
   inline void integrateOnElement(const Array<Real> & f, Real * intf,
                                  UInt nb_degree_of_freedom, const UInt elem,
                                  const GhostType & ghost_type) const;
 
   /// integrate f for all elements of type "type"
   template <ElementType type>
   void integrate(const Array<Real> & in_f, Array<Real> & intf,
                  UInt nb_degree_of_freedom, const GhostType & ghost_type,
                  const Array<UInt> & filter_elements) const;
 
   /// integrate scalar field in_f
   template <ElementType type, UInt polynomial_degree>
   Real integrate(const Array<Real> & in_f,
                  const GhostType & ghost_type = _not_ghost) const;
 
   /// integrate partially around a quadrature point (@f$ intf_q = f_q * J_q *
   /// w_q @f$)
   template <ElementType type>
   Real integrate(const Vector<Real> & in_f, UInt index,
                  const GhostType & ghost_type) const;
 
   /// integrate scalar field in_f
   template <ElementType type>
   Real integrate(const Array<Real> & in_f, const GhostType & ghost_type,
                  const Array<UInt> & filter_elements) const;
 
   /// integrate a field without using the pre-computed values
   template <ElementType type, UInt polynomial_degree>
   void integrate(const Array<Real> & in_f, Array<Real> & intf,
                  UInt nb_degree_of_freedom, const GhostType & ghost_type) const;
 
   /// integrate partially around a quadrature point (@f$ intf_q = f_q * J_q *
   /// w_q @f$)
   template <ElementType type>
   void integrateOnIntegrationPoints(const Array<Real> & in_f,
                                     Array<Real> & intf,
                                     UInt nb_degree_of_freedom,
                                     const GhostType & ghost_type,
                                     const Array<UInt> & filter_elements) const;
 
   /// return a matrix with quadrature points natural coordinates
   template <ElementType type>
   const Matrix<Real> & getIntegrationPoints(const GhostType & ghost_type) const;
 
   /// return number of quadrature points
   template <ElementType type>
   UInt getNbIntegrationPoints(const GhostType & ghost_type) const;
 
   template <ElementType type, UInt n> Matrix<Real> getIntegrationPoints() const;
   template <ElementType type, UInt n>
   Vector<Real> getIntegrationWeights() const;
 
 protected:
   friend struct integrator::details::GaussIntegratorComputeJacobiansHelper<
       kind>;
 
   template <ElementType type>
   void computeJacobiansOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & quad_points,
       Array<Real> & jacobians, const GhostType & ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   void computeJacobiansOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & quad_points,
       Array<Real> & jacobians, const ElementType & type,
       const GhostType & ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// precompute jacobians on elements of type "type"
   template <ElementType type>
   void precomputeJacobiansOnQuadraturePoints(const Array<Real> & nodes,
                                              const GhostType & ghost_type);
 
   // multiply the jacobians by the integration weights and stores the results in
   // jacobians
   template <ElementType type, UInt polynomial_degree>
   void multiplyJacobiansByWeights(
       Array<Real> & jacobians,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// compute the vector of quadrature points natural coordinates
   template <ElementType type>
   void computeQuadraturePoints(const GhostType & ghost_type);
 
   /// check that the jacobians are not negative
   template <ElementType type>
   void checkJacobians(const GhostType & ghost_type) const;
 
   /// internal integrate partially around a quadrature point (@f$ intf_q = f_q *
   /// J_q *
   /// w_q @f$)
   void integrateOnIntegrationPoints(const Array<Real> & in_f,
                                     Array<Real> & intf,
                                     UInt nb_degree_of_freedom,
                                     const Array<Real> & jacobians,
                                     UInt nb_element) const;
 
   void integrate(const Array<Real> & in_f, Array<Real> & intf,
                  UInt nb_degree_of_freedom, const Array<Real> & jacobians,
                  UInt nb_element) const;
 
 public:
   /// compute the jacobians on quad points for a given element
   template <ElementType type>
   void
   computeJacobianOnQuadPointsByElement(const Matrix<Real> & node_coords,
                                        const Matrix<Real> & integration_points,
                                        Vector<Real> & jacobians) const;
 
 public:
   void onElementsAdded(const Array<Element> & elements) override;
 
   template <ElementType type>
   void onElementsAddedByType(const Array<UInt> & new_elements,
                              const GhostType & ghost_type);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// integrate the field f with the jacobian jac -> inte
   inline void integrate(Real * f, Real * jac, Real * inte,
                         UInt nb_degree_of_freedom,
                         UInt nb_quadrature_points) const;
 
 private:
   /// ElementTypeMap of the quadrature points
   ElementTypeMap<Matrix<Real>> quadrature_points;
 };
 
 } // namespace akantu
 
 #include "integrator_gauss_inline_impl.cc"
 
 #endif /* __AKANTU_INTEGRATOR_GAUSS_HH__ */
diff --git a/src/fe_engine/integrator_gauss_inline_impl.cc b/src/fe_engine/integrator_gauss_inline_impl.cc
index b13c68154..c24de5cb5 100644
--- a/src/fe_engine/integrator_gauss_inline_impl.cc
+++ b/src/fe_engine/integrator_gauss_inline_impl.cc
@@ -1,756 +1,760 @@
 /**
  * @file   integrator_gauss_inline_impl.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Feb 15 2011
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  inline function of gauss integrator
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "fe_engine.hh"
 #include "mesh_iterators.hh"
 #if defined(AKANTU_DEBUG_TOOLS)
 #include "aka_debug_tools.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 namespace debug {
   struct IntegratorGaussException : public Exception {};
 } // namespace debug
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::integrateOnElement(
     const Array<Real> & f, Real * intf, UInt nb_degree_of_freedom,
     const UInt elem, const GhostType & ghost_type) const {
   Array<Real> & jac_loc = jacobians(type, ghost_type);
 
   UInt nb_quadrature_points = ElementClass<type>::getNbQuadraturePoints();
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
                       "The vector f do not have the good number of component.");
 
   Real * f_val = f.storage() + elem * f.getNbComponent();
   Real * jac_val = jac_loc.storage() + elem * nb_quadrature_points;
 
   integrate(f_val, jac_val, intf, nb_degree_of_freedom, nb_quadrature_points);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline Real IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Vector<Real> & in_f, UInt index, const GhostType & ghost_type) const {
   const Array<Real> & jac_loc = jacobians(type, ghost_type);
 
   UInt nb_quadrature_points =
       GaussIntegrationElement<type>::getNbQuadraturePoints();
   AKANTU_DEBUG_ASSERT(in_f.size() == nb_quadrature_points,
                       "The vector f do not have nb_quadrature_points entries.");
 
   Real * jac_val = jac_loc.storage() + index * nb_quadrature_points;
   Real intf;
 
   integrate(in_f.storage(), jac_val, &intf, 1, nb_quadrature_points);
 
   return intf;
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom,
     UInt nb_quadrature_points) const {
   memset(inte, 0, nb_degree_of_freedom * sizeof(Real));
 
   Real * cjac = jac;
   for (UInt q = 0; q < nb_quadrature_points; ++q) {
     for (UInt dof = 0; dof < nb_degree_of_freedom; ++dof) {
       inte[dof] += *f * *cjac;
       ++f;
     }
     ++cjac;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline const Matrix<Real> &
 IntegratorGauss<kind, IntegrationOrderFunctor>::getIntegrationPoints(
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_ASSERT(
       quadrature_points.exists(type, ghost_type),
       "Quadrature points for type "
           << quadrature_points.printType(type, ghost_type)
           << " have not been initialized."
           << " Did you use 'computeQuadraturePoints' function ?");
   return quadrature_points(type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline UInt
 IntegratorGauss<kind, IntegrationOrderFunctor>::getNbIntegrationPoints(
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_ASSERT(
       quadrature_points.exists(type, ghost_type),
       "Quadrature points for type "
           << quadrature_points.printType(type, ghost_type)
           << " have not been initialized."
           << " Did you use 'computeQuadraturePoints' function ?");
   return quadrature_points(type, ghost_type).cols();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 inline Matrix<Real>
 IntegratorGauss<kind, IntegrationOrderFunctor>::getIntegrationPoints() const {
   return GaussIntegrationElement<type,
                                  polynomial_degree>::getQuadraturePoints();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 inline Vector<Real>
 IntegratorGauss<kind, IntegrationOrderFunctor>::getIntegrationWeights() const {
   return GaussIntegrationElement<type, polynomial_degree>::getWeights();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void
 IntegratorGauss<kind, IntegrationOrderFunctor>::computeQuadraturePoints(
     const GhostType & ghost_type) {
   Matrix<Real> & quads = quadrature_points(type, ghost_type);
   const UInt polynomial_degree =
       IntegrationOrderFunctor::template getOrder<type>();
   quads =
       GaussIntegrationElement<type, polynomial_degree>::getQuadraturePoints();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::
     computeJacobianOnQuadPointsByElement(const Matrix<Real> & node_coords,
                                          const Matrix<Real> & quad,
                                          Vector<Real> & jacobians) const {
   // jacobian
   ElementClass<type>::computeJacobian(quad, node_coords, jacobians);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 IntegratorGauss<kind, IntegrationOrderFunctor>::IntegratorGauss(
-    const Mesh & mesh, const ID & id, const MemoryID & memory_id)
-    : Integrator(mesh, id, memory_id) {
+    const Mesh & mesh, UInt spatial_dimension, const ID & id,
+    const MemoryID & memory_id)
+    : Integrator(mesh, spatial_dimension, id, memory_id) {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::checkJacobians(
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_quadrature_points = this->quadrature_points(type, ghost_type).cols();
 
   UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
 
   Real * jacobians_val = jacobians(type, ghost_type).storage();
 
   for (UInt i = 0; i < nb_element * nb_quadrature_points;
        ++i, ++jacobians_val) {
     if (*jacobians_val < 0)
       AKANTU_CUSTOM_EXCEPTION_INFO(debug::IntegratorGaussException{},
                                    "Negative jacobian computed,"
                                        << " possible problem in the element "
                                           "node ordering (Quadrature Point "
                                        << i % nb_quadrature_points << ":"
                                        << i / nb_quadrature_points << ":"
                                        << type << ":" << ghost_type << ")");
   }
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     computeJacobiansOnIntegrationPoints(
         const Array<Real> & nodes, const Matrix<Real> & quad_points,
         Array<Real> & jacobians, const GhostType & ghost_type,
         const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_quadrature_points = quad_points.cols();
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   jacobians.resize(nb_element * nb_quadrature_points);
 
   auto jacobians_it =
       jacobians.begin_reinterpret(nb_quadrature_points, nb_element);
   auto jacobians_begin = jacobians_it;
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type,
                                        filter_elements);
 
   auto x_it = x_el.begin(spatial_dimension, nb_nodes_per_element);
 
   nb_element = x_el.size();
 
   //  Matrix<Real> local_coord(spatial_dimension, nb_nodes_per_element);
   for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) {
     const Matrix<Real> & x = *x_it;
     if (filter_elements != empty_filter) {
       jacobians_it = jacobians_begin + filter_elements(elem);
     }
 
     Vector<Real> & J = *jacobians_it;
     computeJacobianOnQuadPointsByElement<type>(x, quad_points, J);
 
     if (filter_elements == empty_filter) {
       ++jacobians_it;
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 template <>
 template <ElementType type>
 void IntegratorGauss<_ek_structural, DefaultIntegrationOrderFunctor>::
     computeJacobiansOnIntegrationPoints(
         const Array<Real> & nodes, const Matrix<Real> & quad_points,
         Array<Real> & jacobians, const GhostType & ghost_type,
         const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   const UInt spatial_dimension = mesh.getSpatialDimension();
   const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   const UInt nb_quadrature_points = quad_points.cols();
   const UInt nb_dofs = ElementClass<type>::getNbDegreeOfFreedom();
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   jacobians.resize(nb_element * nb_quadrature_points);
 
   auto jacobians_it =
       jacobians.begin_reinterpret(nb_quadrature_points, nb_element);
   auto jacobians_begin = jacobians_it;
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type,
                                        filter_elements);
 
   auto x_it = x_el.begin(spatial_dimension, nb_nodes_per_element);
 
   nb_element = x_el.size();
 
   const bool has_extra_normal =
       mesh.hasData<Real>("extra_normal", type, ghost_type);
   Array<Real>::const_vector_iterator extra_normal, extra_normal_begin;
   if (has_extra_normal) {
     extra_normal = mesh.getData<Real>("extra_normal", type, ghost_type)
                        .begin(spatial_dimension);
     extra_normal_begin = extra_normal;
   }
 
   //  Matrix<Real> local_coord(spatial_dimension, nb_nodes_per_element);
   for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) {
     if (filter_elements != empty_filter) {
       jacobians_it = jacobians_begin + filter_elements(elem);
       extra_normal = extra_normal_begin + filter_elements(elem);
     }
 
     const Matrix<Real> & X = *x_it;
     Vector<Real> & J = *jacobians_it;
     Matrix<Real> R(nb_dofs, nb_dofs);
 
     if (has_extra_normal)
       ElementClass<type>::computeRotationMatrix(R, X, *extra_normal);
     else
       ElementClass<type>::computeRotationMatrix(R, X, Vector<Real>(X.rows()));
     // Extracting relevant lines
     auto x = (R.block(0, 0, spatial_dimension, spatial_dimension) * X)
                  .block(0, 0, ElementClass<type>::getNaturalSpaceDimension(),
                         ElementClass<type>::getNbNodesPerElement());
 
     computeJacobianOnQuadPointsByElement<type>(x, quad_points, J);
 
     if (filter_elements == empty_filter) {
       ++jacobians_it;
       ++extra_normal;
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 #endif
 
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_COHESIVE_ELEMENT)
 template <>
 template <ElementType type>
 void IntegratorGauss<_ek_cohesive, DefaultIntegrationOrderFunctor>::
     computeJacobiansOnIntegrationPoints(
         const Array<Real> & nodes, const Matrix<Real> & quad_points,
         Array<Real> & jacobians, const GhostType & ghost_type,
         const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_quadrature_points = quad_points.cols();
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   jacobians.resize(nb_element * nb_quadrature_points);
 
   auto jacobians_begin =
       jacobians.begin_reinterpret(nb_quadrature_points, nb_element);
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type,
                                        filter_elements);
 
   auto x_it = x_el.begin(spatial_dimension, nb_nodes_per_element);
 
   UInt nb_nodes_per_subelement = nb_nodes_per_element / 2;
   Matrix<Real> x(spatial_dimension, nb_nodes_per_subelement);
 
   nb_element = x_el.size();
   UInt l_el = 0;
   auto compute = [&](const auto & el) {
     Vector<Real> J(jacobians_begin[el]);
     Matrix<Real> X(x_it[l_el]);
     ++l_el;
 
     for (UInt n = 0; n < nb_nodes_per_subelement; ++n)
       Vector<Real>(x(n)) =
           (Vector<Real>(X(n)) + Vector<Real>(X(n + nb_nodes_per_subelement))) /
           2.;
 
     if (type == _cohesive_1d_2)
       J(0) = 1;
     else
       this->computeJacobianOnQuadPointsByElement<type>(x, quad_points, J);
   };
 
   for_each_element(nb_element, filter_elements, compute);
 
   AKANTU_DEBUG_OUT();
 }
 #endif
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     precomputeJacobiansOnQuadraturePoints(const Array<Real> & nodes,
                                           const GhostType & ghost_type) {
   AKANTU_DEBUG_IN();
 
   Array<Real> & jacobians_tmp = jacobians.alloc(0, 1, type, ghost_type);
 
   this->computeJacobiansOnIntegrationPoints<type>(
       nodes, quadrature_points(type, ghost_type), jacobians_tmp, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::multiplyJacobiansByWeights(
     Array<Real> & jacobians, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_quadrature_points =
       GaussIntegrationElement<type, polynomial_degree>::getNbQuadraturePoints();
 
   Vector<Real> weights =
       GaussIntegrationElement<type, polynomial_degree>::getWeights();
 
   auto && view = make_view(jacobians, nb_quadrature_points);
 
   if (filter_elements != empty_filter) {
     auto J_it = view.begin();
     for (auto el : filter_elements) {
       Vector<Real> J(J_it[el]);
       J *= weights;
     }
   } else {
     for (auto & J : make_view(jacobians, nb_quadrature_points)) {
       J *= weights;
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom,
     const Array<Real> & jacobians, UInt nb_element) const {
   AKANTU_DEBUG_IN();
 
   intf.resize(nb_element);
   if (nb_element == 0)
     return;
 
   UInt nb_points = jacobians.size() / nb_element;
 
   Array<Real>::const_matrix_iterator J_it;
   Array<Real>::matrix_iterator inte_it;
   Array<Real>::const_matrix_iterator f_it;
 
   f_it = in_f.begin_reinterpret(nb_degree_of_freedom, nb_points, nb_element);
   inte_it = intf.begin_reinterpret(nb_degree_of_freedom, 1, nb_element);
   J_it = jacobians.begin_reinterpret(nb_points, 1, nb_element);
 
   for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) {
     const Matrix<Real> & f = *f_it;
     const Matrix<Real> & J = *J_it;
     Matrix<Real> & inte_f = *inte_it;
 
     inte_f.mul<false, false>(f, J);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type),
                       "No jacobians for the type "
                           << jacobians.printType(type, ghost_type));
 
   const Array<Real> & jac_loc = jacobians(type, ghost_type);
   if (filter_elements != empty_filter) {
     UInt nb_element = filter_elements.size();
     Array<Real> * filtered_J = new Array<Real>(0, jac_loc.getNbComponent());
     FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type,
                                   filter_elements);
     this->integrate(in_f, intf, nb_degree_of_freedom, *filtered_J, nb_element);
     delete filtered_J;
   } else {
     UInt nb_element = mesh.getNbElement(type, ghost_type);
     this->integrate(in_f, intf, nb_degree_of_freedom, jac_loc, nb_element);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom,
     const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   Matrix<Real> quads = this->getIntegrationPoints<type, polynomial_degree>();
 
   Array<Real> jacobians;
   this->computeJacobiansOnIntegrationPoints<type>(mesh.getNodes(), quads,
                                                   jacobians, ghost_type);
   this->multiplyJacobiansByWeights<type, polynomial_degree>(jacobians);
 
   this->integrate(in_f, intf, nb_degree_of_freedom, jacobians,
                   mesh.getNbElement(type, ghost_type));
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type, UInt polynomial_degree>
 Real IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, const GhostType & ghost_type) const {
   AKANTU_DEBUG_IN();
 
   Array<Real> intfv(0, 1);
   integrate<type, polynomial_degree>(in_f, intfv, 1, ghost_type);
 
   Real res = Math::reduce(intfv);
 
   AKANTU_DEBUG_OUT();
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 Real IntegratorGauss<kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & in_f, const GhostType & ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type),
                       "No jacobians for the type "
                           << jacobians.printType(type, ghost_type));
 
   Array<Real> intfv(0, 1);
   integrate<type>(in_f, intfv, 1, ghost_type, filter_elements);
 
   Real res = Math::reduce(intfv);
 
   AKANTU_DEBUG_OUT();
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     integrateOnIntegrationPoints(const Array<Real> & in_f, Array<Real> & intf,
                                  UInt nb_degree_of_freedom,
                                  const Array<Real> & jacobians,
                                  UInt nb_element) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_points = jacobians.size() / nb_element;
 
   intf.resize(nb_element * nb_points);
 
   auto J_it = jacobians.begin();
   auto f_it = in_f.begin(nb_degree_of_freedom);
   auto inte_it = intf.begin(nb_degree_of_freedom);
 
   for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) {
     const Real & J = *J_it;
     const Vector<Real> & f = *f_it;
     Vector<Real> & inte_f = *inte_it;
 
     inte_f = f;
     inte_f *= J;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     integrateOnIntegrationPoints(const Array<Real> & in_f, Array<Real> & intf,
                                  UInt nb_degree_of_freedom,
                                  const GhostType & ghost_type,
                                  const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type),
                       "No jacobians for the type "
                           << jacobians.printType(type, ghost_type));
 
   const Array<Real> & jac_loc = this->jacobians(type, ghost_type);
 
   if (filter_elements != empty_filter) {
 
     UInt nb_element = filter_elements.size();
     Array<Real> * filtered_J = new Array<Real>(0, jac_loc.getNbComponent());
     FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type,
                                   filter_elements);
 
     this->integrateOnIntegrationPoints(in_f, intf, nb_degree_of_freedom,
                                        *filtered_J, nb_element);
   } else {
     UInt nb_element = mesh.getNbElement(type, ghost_type);
     this->integrateOnIntegrationPoints(in_f, intf, nb_degree_of_freedom,
                                        jac_loc, nb_element);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void
 IntegratorGauss<kind, IntegrationOrderFunctor>::onElementsAddedByType(
     const Array<UInt> & elements, const GhostType & ghost_type) {
   const auto & nodes = mesh.getNodes();
 
   if (not quadrature_points.exists(type, ghost_type)) {
     computeQuadraturePoints<type>(ghost_type);
   }
 
   if (not jacobians.exists(type, ghost_type))
     jacobians.alloc(0, 1, type, ghost_type);
 
   this->computeJacobiansOnIntegrationPoints(
       nodes, quadrature_points(type, ghost_type), jacobians(type, ghost_type),
       type, ghost_type, elements);
 
   constexpr UInt polynomial_degree =
       IntegrationOrderFunctor::template getOrder<type>();
 
   multiplyJacobiansByWeights<type, polynomial_degree>(
       this->jacobians(type, ghost_type), elements);
 }
 
 /* -------------------------------------------------------------------------- */
 namespace integrator {
   namespace details {
     template <ElementKind kind> struct IntegratorOnElementsAddedHelper {};
 
 #define ON_ELEMENT_ADDED(type)                                                 \
   integrator.template onElementsAddedByType<type>(elements, ghost_type);
 
 #define AKANTU_SPECIALIZE_ON_ELEMENT_ADDED_HELPER(kind)                        \
   template <> struct IntegratorOnElementsAddedHelper<kind> {                   \
     template <class I>                                                         \
     static void call(I & integrator, const Array<UInt> & elements,             \
                      const ElementType & type, const GhostType & ghost_type) { \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(ON_ELEMENT_ADDED, kind);                \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_ON_ELEMENT_ADDED_HELPER)
 
 #undef AKANTU_SPECIALIZE_ON_ELEMENT_ADDED_HELPER
 #undef ON_ELEMENT_ADDED
   } // namespace details
 } // namespace integrator
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::onElementsAdded(
     const Array<Element> & new_elements) {
   for (auto elements_range : MeshElementsByTypes(new_elements)) {
     auto type = elements_range.getType();
     auto ghost_type = elements_range.getGhostType();
 
+    if (mesh.getSpatialDimension(type) != _spatial_dimension)
+      continue;
+    
     if (mesh.getKind(type) != kind)
       continue;
 
     integrator::details::IntegratorOnElementsAddedHelper<kind>::call(
         *this, elements_range.getElements(), type, ghost_type);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::initIntegrator(
     const Array<Real> & nodes, const GhostType & ghost_type) {
   computeQuadraturePoints<type>(ghost_type);
   precomputeJacobiansOnQuadraturePoints<type>(nodes, ghost_type);
   checkJacobians<type>(ghost_type);
   constexpr UInt polynomial_degree =
       IntegrationOrderFunctor::template getOrder<type>();
   multiplyJacobiansByWeights<type, polynomial_degree>(
       this->jacobians(type, ghost_type));
 }
 
 namespace integrator {
   namespace details {
     template <ElementKind kind> struct GaussIntegratorInitHelper {};
 
 #define INIT_INTEGRATOR(type)                                                  \
   _int.template initIntegrator<type>(nodes, ghost_type)
 
 #define AKANTU_GAUSS_INTERGRATOR_INIT_HELPER(kind)                             \
   template <> struct GaussIntegratorInitHelper<kind> {                         \
     template <ElementKind k, class IOF>                                        \
     static void call(IntegratorGauss<k, IOF> & _int,                           \
                      const Array<Real> & nodes, const ElementType & type,      \
                      const GhostType & ghost_type) {                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INIT_INTEGRATOR, kind);                 \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_GAUSS_INTERGRATOR_INIT_HELPER)
 
 #undef AKANTU_GAUSS_INTERGRATOR_INIT_HELPER
 #undef INIT_INTEGRATOR
   } // namespace details
 } // namespace integrator
 
 template <ElementKind kind, class IntegrationOrderFunctor>
 inline void IntegratorGauss<kind, IntegrationOrderFunctor>::initIntegrator(
     const Array<Real> & nodes, const ElementType & type,
     const GhostType & ghost_type) {
   integrator::details::GaussIntegratorInitHelper<kind>::call(*this, nodes, type,
                                                              ghost_type);
 }
 
 namespace integrator {
   namespace details {
     template <ElementKind kind> struct GaussIntegratorComputeJacobiansHelper {};
 
 #define AKANTU_COMPUTE_JACOBIANS(type)                                         \
   _int.template computeJacobiansOnIntegrationPoints<type>(                     \
       nodes, quad_points, jacobians, ghost_type, filter_elements);
 
 #define AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS(kind)                       \
   template <> struct GaussIntegratorComputeJacobiansHelper<kind> {             \
     template <ElementKind k, class IOF>                                        \
     static void                                                                \
     call(const IntegratorGauss<k, IOF> & _int, const Array<Real> & nodes,      \
          const Matrix<Real> & quad_points, Array<Real> & jacobians,            \
          const ElementType & type, const GhostType & ghost_type,               \
          const Array<UInt> & filter_elements) {                                \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(AKANTU_COMPUTE_JACOBIANS, kind);        \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS)
 
 #undef AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS
 #undef AKANTU_COMPUTE_JACOBIANS
   } // namespace details
 } // namespace integrator
 
 template <ElementKind kind, class IntegrationOrderFunctor>
 void IntegratorGauss<kind, IntegrationOrderFunctor>::
     computeJacobiansOnIntegrationPoints(
         const Array<Real> & nodes, const Matrix<Real> & quad_points,
         Array<Real> & jacobians, const ElementType & type,
         const GhostType & ghost_type,
         const Array<UInt> & filter_elements) const {
   integrator::details::GaussIntegratorComputeJacobiansHelper<kind>::call(
       *this, nodes, quad_points, jacobians, type, ghost_type, filter_elements);
 }
 
 } // namespace akantu
diff --git a/src/fe_engine/shape_cohesive.hh b/src/fe_engine/shape_cohesive.hh
index 2d82583b8..953611653 100644
--- a/src/fe_engine/shape_cohesive.hh
+++ b/src/fe_engine/shape_cohesive.hh
@@ -1,184 +1,185 @@
 /**
  * @file   shape_cohesive.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Tue Feb 15 2011
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  shape functions for cohesive elements description
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 #include "shape_lagrange.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_SHAPE_COHESIVE_HH__
 #define __AKANTU_SHAPE_COHESIVE_HH__
 
 namespace akantu {
 
 struct CohesiveReduceFunctionMean {
   inline Real operator()(Real u_plus, Real u_minus) {
     return .5 * (u_plus + u_minus);
   }
 };
 
 struct CohesiveReduceFunctionOpening {
   inline Real operator()(Real u_plus, Real u_minus) {
     return (u_plus - u_minus);
   }
 };
 
 template <> class ShapeLagrange<_ek_cohesive> : public ShapeLagrangeBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
-  ShapeLagrange(const Mesh & mesh, const ID & id = "shape_cohesive",
+  ShapeLagrange(const Mesh & mesh, UInt spatial_dimension,
+                const ID & id = "shape_cohesive",
                 const MemoryID & memory_id = 0);
 
   ~ShapeLagrange() override = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   inline void initShapeFunctions(const Array<Real> & nodes,
                                  const Matrix<Real> & integration_points,
                                  const ElementType & type,
                                  const GhostType & ghost_type);
 
   /// extract the nodal values and store them per element
   template <ElementType type, class ReduceFunction>
   void extractNodalToElementField(
       const Array<Real> & nodal_f, Array<Real> & elemental_f,
       const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// computes the shape functions derivatives for given interpolation points
   template <ElementType type>
   void computeShapeDerivativesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
       Array<Real> & shape_derivatives, const GhostType & ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   void computeShapeDerivativesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
       Array<Real> & shape_derivatives, const ElementType & type,
       const GhostType & ghost_type,
       const Array<UInt> & filter_elements) const override;
 
   /// pre compute all shapes on the element integration points from natural
   /// coordinates
   template <ElementType type>
   void precomputeShapesOnIntegrationPoints(const Array<Real> & nodes,
                                            GhostType ghost_type);
 
   /// pre compute all shape derivatives on the element integration points from
   /// natural coordinates
   template <ElementType type>
   void precomputeShapeDerivativesOnIntegrationPoints(const Array<Real> & nodes,
                                                      GhostType ghost_type);
 
   /// interpolate nodal values on the integration points
   template <ElementType type, class ReduceFunction>
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       const GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   template <ElementType type>
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       const GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const {
     interpolateOnIntegrationPoints<type, CohesiveReduceFunctionMean>(
         u, uq, nb_degree_of_freedom, ghost_type, filter_elements);
   }
 
   /// compute the gradient of u on the integration points in the natural
   /// coordinates
   template <ElementType type>
   void gradientOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const {
     variationOnIntegrationPoints<type, CohesiveReduceFunctionMean>(
         u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements);
   }
 
   /* ------------------------------------------------------------------------ */
   template <ElementType type>
   void computeBtD(const Array<Real> & /*Ds*/, Array<Real> & /*BtDs*/,
                   GhostType /*ghost_type*/,
                   const Array<UInt> & /*filter_elements*/) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   template <ElementType type>
   void computeBtDB(const Array<Real> & /*Ds*/, Array<Real> & /*BtDBs*/,
                    UInt /*order_d*/, GhostType /*ghost_type*/,
                    const Array<UInt> & /*filter_elements*/) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// multiply a field by shape functions
   template <ElementType type>
   void
   computeNtb(const Array<Real> & /*bs*/, Array<Real> & /*Ntbs*/,
              const GhostType & /*ghost_type*/,
              const Array<UInt> & /*filter_elements*/ = empty_filter) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   /* ------------------------------------------------------------------------ */
   /// compute the gradient of u on the integration points
   template <ElementType type, class ReduceFunction>
   void variationOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// compute the normals to the field u on integration points
   template <ElementType type, class ReduceFunction>
   void computeNormalsOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & normals_u,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 };
 
 /// standard output stream operator
 template <class ShapeFunction>
 inline std::ostream & operator<<(std::ostream & stream,
                                  const ShapeCohesive<ShapeFunction> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "shape_cohesive_inline_impl.cc"
 
 #endif /* __AKANTU_SHAPE_COHESIVE_HH__ */
diff --git a/src/fe_engine/shape_cohesive_inline_impl.cc b/src/fe_engine/shape_cohesive_inline_impl.cc
index 9b570d6ba..35aa13ca2 100644
--- a/src/fe_engine/shape_cohesive_inline_impl.cc
+++ b/src/fe_engine/shape_cohesive_inline_impl.cc
@@ -1,329 +1,331 @@
 /**
  * @file   shape_cohesive_inline_impl.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Fri Feb 03 2012
  * @date last modification: Mon Feb 19 2018
  *
  * @brief  ShapeCohesive inline implementation
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "mesh_iterators.hh"
 #include "shape_cohesive.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_SHAPE_COHESIVE_INLINE_IMPL_CC__
 #define __AKANTU_SHAPE_COHESIVE_INLINE_IMPL_CC__
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 inline ShapeLagrange<_ek_cohesive>::ShapeLagrange(const Mesh & mesh,
+                                                  UInt spatial_dimension,
                                                   const ID & id,
                                                   const MemoryID & memory_id)
-    : ShapeLagrangeBase(mesh, _ek_cohesive, id, memory_id) {}
+    : ShapeLagrangeBase(mesh, spatial_dimension, _ek_cohesive, id, memory_id) {}
 
 #define INIT_SHAPE_FUNCTIONS(type)                                             \
   setIntegrationPointsByType<type>(integration_points, ghost_type);            \
   precomputeShapesOnIntegrationPoints<type>(nodes, ghost_type);                \
   precomputeShapeDerivativesOnIntegrationPoints<type>(nodes, ghost_type);
 
 /* -------------------------------------------------------------------------- */
 inline void ShapeLagrange<_ek_cohesive>::initShapeFunctions(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     const ElementType & type, const GhostType & ghost_type) {
   AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
 }
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void ShapeLagrange<_ek_cohesive>::computeShapeDerivativesOnIntegrationPoints(
     const Array<Real> &, const Matrix<Real> & integration_points,
     Array<Real> & shape_derivatives, const GhostType & ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt size_of_shapesd = ElementClass<type>::getShapeDerivativesSize();
   UInt spatial_dimension = ElementClass<type>::getNaturalSpaceDimension();
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   UInt nb_points = integration_points.cols();
   UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
 
   AKANTU_DEBUG_ASSERT(shape_derivatives.getNbComponent() == size_of_shapesd,
                       "The shapes_derivatives array does not have the correct "
                           << "number of component");
 
   shape_derivatives.resize(nb_element * nb_points);
   Real * shapesd_val = shape_derivatives.storage();
 
   auto compute = [&](const auto & el) {
     auto ptr = shapesd_val + el * nb_points * size_of_shapesd;
     Tensor3<Real> B(ptr, spatial_dimension, nb_nodes_per_element, nb_points);
     ElementClass<type>::computeDNDS(integration_points, B);
   };
 
   for_each_element(nb_element, filter_elements, compute);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 inline void
 ShapeLagrange<_ek_cohesive>::computeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     Array<Real> & shape_derivatives, const ElementType & type,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
 #define AKANTU_COMPUTE_SHAPES(type)                                            \
   computeShapeDerivativesOnIntegrationPoints<type>(                            \
       nodes, integration_points, shape_derivatives, ghost_type,                \
       filter_elements);
 
   AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES);
 
 #undef AKANTU_COMPUTE_SHAPES
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void ShapeLagrange<_ek_cohesive>::precomputeShapesOnIntegrationPoints(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   Matrix<Real> & natural_coords = integration_points(type, ghost_type);
   UInt size_of_shapes = ElementClass<type>::getShapeSize();
 
   Array<Real> & shapes_tmp =
       shapes.alloc(0, size_of_shapes, itp_type, ghost_type);
 
   this->computeShapesOnIntegrationPoints<type>(nodes, natural_coords,
                                                shapes_tmp, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void ShapeLagrange<_ek_cohesive>::precomputeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   Matrix<Real> & natural_coords = integration_points(type, ghost_type);
   UInt size_of_shapesd = ElementClass<type>::getShapeDerivativesSize();
 
   Array<Real> & shapes_derivatives_tmp =
       shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type);
 
   this->computeShapeDerivativesOnIntegrationPoints<type>(
       nodes, natural_coords, shapes_derivatives_tmp, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type, class ReduceFunction>
 void ShapeLagrange<_ek_cohesive>::extractNodalToElementField(
     const Array<Real> & nodal_f, Array<Real> & elemental_f,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_nodes_per_itp_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
   UInt nb_degree_of_freedom = nodal_f.getNbComponent();
   UInt nb_element = this->mesh.getNbElement(type, ghost_type);
 
   const auto & conn_array = this->mesh.getConnectivity(type, ghost_type);
   auto conn = conn_array.begin(conn_array.getNbComponent() / 2, 2);
 
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 
   elemental_f.resize(nb_element);
 
   Array<Real>::matrix_iterator u_it =
       elemental_f.begin(nb_degree_of_freedom, nb_nodes_per_itp_element);
 
   ReduceFunction reduce_function;
 
   auto compute = [&](const auto & el) {
     Matrix<Real> & u = *u_it;
     Matrix<UInt> el_conn(conn[el]);
 
     // compute the average/difference of the nodal field loaded from cohesive
     // element
     for (UInt n = 0; n < el_conn.rows(); ++n) {
       UInt node_plus = el_conn(n, 0);
       UInt node_minus = el_conn(n, 1);
       for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
         Real u_plus = nodal_f(node_plus, d);
         Real u_minus = nodal_f(node_minus, d);
         u(d, n) = reduce_function(u_plus, u_minus);
       }
     }
 
     ++u_it;
   };
 
   for_each_element(nb_element, filter_elements, compute);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type, class ReduceFunction>
 void ShapeLagrange<_ek_cohesive>::interpolateOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
 
   AKANTU_DEBUG_ASSERT(this->shapes.exists(itp_type, ghost_type),
                       "No shapes for the type "
                           << this->shapes.printType(itp_type, ghost_type));
 
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
   Array<Real> u_el(0, nb_degree_of_freedom * nb_nodes_per_element);
   this->extractNodalToElementField<type, ReduceFunction>(in_u, u_el, ghost_type,
                                                          filter_elements);
 
   this->template interpolateElementalFieldOnIntegrationPoints<type>(
       u_el, out_uq, ghost_type, shapes(itp_type, ghost_type), filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type, class ReduceFunction>
 void ShapeLagrange<_ek_cohesive>::variationOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
 
   AKANTU_DEBUG_ASSERT(
       this->shapes_derivatives.exists(itp_type, ghost_type),
       "No shapes for the type "
           << this->shapes_derivatives.printType(itp_type, ghost_type));
 
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
   Array<Real> u_el(0, nb_degree_of_freedom * nb_nodes_per_element);
   this->extractNodalToElementField<type, ReduceFunction>(in_u, u_el, ghost_type,
                                                          filter_elements);
 
   this->template gradientElementalFieldOnIntegrationPoints<type>(
       u_el, nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type),
       filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type, class ReduceFunction>
 void ShapeLagrange<_ek_cohesive>::computeNormalsOnIntegrationPoints(
     const Array<Real> & u, Array<Real> & normals_u, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_element = this->mesh.getNbElement(type, ghost_type);
   UInt nb_points = this->integration_points(type, ghost_type).cols();
   UInt spatial_dimension = this->mesh.getSpatialDimension();
 
   if (filter_elements != empty_filter)
     nb_element = filter_elements.size();
 
   normals_u.resize(nb_points * nb_element);
 
   Array<Real> tangents_u(0, (spatial_dimension * (spatial_dimension - 1)));
 
   if (spatial_dimension > 1) {
     tangents_u.resize(nb_element * nb_points);
     this->template variationOnIntegrationPoints<type, ReduceFunction>(
         u, tangents_u, spatial_dimension, ghost_type, filter_elements);
   }
 
   Real * tangent = tangents_u.storage();
 
   if (spatial_dimension == 3) {
     for (auto & normal : make_view(normals_u, spatial_dimension)) {
       Math::vectorProduct3(tangent, tangent + spatial_dimension,
                            normal.storage());
 
       normal /= normal.norm();
       tangent += spatial_dimension * 2;
     }
   } else if (spatial_dimension == 2) {
     for (auto & normal : make_view(normals_u, spatial_dimension)) {
       Vector<Real> a1(tangent, spatial_dimension);
 
       normal(0) = -a1(1);
       normal(1) = a1(0);
       normal.normalize();
 
       tangent += spatial_dimension;
     }
   } else if (spatial_dimension == 1) {
     const auto facet_type = Mesh::getFacetType(type);
     const auto & mesh_facets = mesh.getMeshFacets();
     const auto & facets = mesh_facets.getSubelementToElement(type, ghost_type);
     const auto & segments =
         mesh_facets.getElementToSubelement(facet_type, ghost_type);
 
     Real values[2];
 
     for (auto el : arange(nb_element)) {
       if (filter_elements != empty_filter)
         el = filter_elements(el);
 
       for (UInt p = 0; p < 2; ++p) {
         Element facet = facets(el, p);
         Element segment = segments(facet.element)[0];
         Vector<Real> barycenter(values + p, 1);
         mesh.getBarycenter(segment, barycenter);
       }
 
       Real difference = values[0] - values[1];
 
       AKANTU_DEBUG_ASSERT(difference != 0.,
                           "Error in normal computation for cohesive elements");
 
       normals_u(el) = difference / std::abs(difference);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
+
 #endif /* __AKANTU_SHAPE_COHESIVE_INLINE_IMPL_CC__ */
diff --git a/src/fe_engine/shape_functions.cc b/src/fe_engine/shape_functions.cc
index 118a8c635..34ebb6586 100644
--- a/src/fe_engine/shape_functions.cc
+++ b/src/fe_engine/shape_functions.cc
@@ -1,235 +1,235 @@
 /**
  * @file   shape_functions.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Aug 09 2017
  * @date last modification: Wed Oct 11 2017
  *
  * @brief  implementation of th shape functions interface
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "shape_functions.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
-ShapeFunctions::ShapeFunctions(const Mesh & mesh, const ID & id,
-                               const MemoryID & memory_id)
+ShapeFunctions::ShapeFunctions(const Mesh & mesh, UInt spatial_dimension,
+                               const ID & id, const MemoryID & memory_id)
     : Memory(id, memory_id), shapes("shapes_generic", id, memory_id),
       shapes_derivatives("shapes_derivatives_generic", id, memory_id),
-      mesh(mesh) {}
+      mesh(mesh), _spatial_dimension(spatial_dimension) {}
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 inline void
 ShapeFunctions::initElementalFieldInterpolationFromIntegrationPoints(
     const Array<Real> & interpolation_points_coordinates,
     ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
     ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
     const Array<Real> & quadrature_points_coordinates,
     const GhostType & ghost_type, const Array<UInt> & element_filter) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = this->mesh.getSpatialDimension();
   UInt nb_element = this->mesh.getNbElement(type, ghost_type);
   UInt nb_element_filter;
 
   if (element_filter == empty_filter)
     nb_element_filter = nb_element;
   else
     nb_element_filter = element_filter.size();
 
   UInt nb_quad_per_element =
       GaussIntegrationElement<type>::getNbQuadraturePoints();
   UInt nb_interpolation_points_per_elem =
       interpolation_points_coordinates.size() / nb_element;
 
   AKANTU_DEBUG_ASSERT(interpolation_points_coordinates.size() % nb_element == 0,
                       "Number of interpolation points should be a multiple of "
                       "total number of elements");
 
   if (!quad_points_coordinates_inv_matrices.exists(type, ghost_type))
     quad_points_coordinates_inv_matrices.alloc(
         nb_element_filter, nb_quad_per_element * nb_quad_per_element, type,
         ghost_type);
   else
     quad_points_coordinates_inv_matrices(type, ghost_type)
         .resize(nb_element_filter);
 
   if (!interpolation_points_coordinates_matrices.exists(type, ghost_type))
     interpolation_points_coordinates_matrices.alloc(
         nb_element_filter,
         nb_interpolation_points_per_elem * nb_quad_per_element, type,
         ghost_type);
   else
     interpolation_points_coordinates_matrices(type, ghost_type)
         .resize(nb_element_filter);
 
   Array<Real> & quad_inv_mat =
       quad_points_coordinates_inv_matrices(type, ghost_type);
   Array<Real> & interp_points_mat =
       interpolation_points_coordinates_matrices(type, ghost_type);
 
   Matrix<Real> quad_coord_matrix(nb_quad_per_element, nb_quad_per_element);
 
   Array<Real>::const_matrix_iterator quad_coords_it =
       quadrature_points_coordinates.begin_reinterpret(
           spatial_dimension, nb_quad_per_element, nb_element_filter);
 
   Array<Real>::const_matrix_iterator points_coords_begin =
       interpolation_points_coordinates.begin_reinterpret(
           spatial_dimension, nb_interpolation_points_per_elem, nb_element);
 
   Array<Real>::matrix_iterator inv_quad_coord_it =
       quad_inv_mat.begin(nb_quad_per_element, nb_quad_per_element);
 
   Array<Real>::matrix_iterator int_points_mat_it = interp_points_mat.begin(
       nb_interpolation_points_per_elem, nb_quad_per_element);
 
   /// loop over the elements of the current material and element type
   for (UInt el = 0; el < nb_element_filter;
        ++el, ++inv_quad_coord_it, ++int_points_mat_it, ++quad_coords_it) {
     /// matrix containing the quadrature points coordinates
     const Matrix<Real> & quad_coords = *quad_coords_it;
     /// matrix to store the matrix inversion result
     Matrix<Real> & inv_quad_coord_matrix = *inv_quad_coord_it;
 
     /// insert the quad coordinates in a matrix compatible with the
     /// interpolation
     buildElementalFieldInterpolationMatrix<type>(quad_coords,
                                                  quad_coord_matrix);
 
     /// invert the interpolation matrix
     inv_quad_coord_matrix.inverse(quad_coord_matrix);
 
     /// matrix containing the interpolation points coordinates
     const Matrix<Real> & points_coords =
         points_coords_begin[element_filter(el)];
     /// matrix to store the interpolation points coordinates
     /// compatible with these functions
     Matrix<Real> & inv_points_coord_matrix = *int_points_mat_it;
 
     /// insert the quad coordinates in a matrix compatible with the
     /// interpolation
     buildElementalFieldInterpolationMatrix<type>(points_coords,
                                                  inv_points_coord_matrix);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void ShapeFunctions::initElementalFieldInterpolationFromIntegrationPoints(
     const ElementTypeMapArray<Real> & interpolation_points_coordinates,
     ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
     ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
     const ElementTypeMapArray<Real> & quadrature_points_coordinates,
     const ElementTypeMapArray<UInt> * element_filter) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = this->mesh.getSpatialDimension();
 
   for (auto ghost_type : ghost_types) {
     auto types_iterable = mesh.elementTypes(spatial_dimension, ghost_type);
     if (element_filter) {
       types_iterable =
           element_filter->elementTypes(spatial_dimension, ghost_type);
     }
 
     for (auto type : types_iterable) {
       UInt nb_element = mesh.getNbElement(type, ghost_type);
       if (nb_element == 0)
         continue;
 
       const Array<UInt> * elem_filter;
       if (element_filter)
         elem_filter = &((*element_filter)(type, ghost_type));
       else
         elem_filter = &(empty_filter);
 
 #define AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS(type)          \
   this->initElementalFieldInterpolationFromIntegrationPoints<type>(            \
       interpolation_points_coordinates(type, ghost_type),                      \
       interpolation_points_coordinates_matrices,                               \
       quad_points_coordinates_inv_matrices,                                    \
       quadrature_points_coordinates(type, ghost_type), ghost_type,             \
       *elem_filter)
 
       AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(
           AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS);
 #undef AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void ShapeFunctions::interpolateElementalFieldFromIntegrationPoints(
     const ElementTypeMapArray<Real> & field,
     const ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
     const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
     ElementTypeMapArray<Real> & result, const GhostType & ghost_type,
     const ElementTypeMapArray<UInt> * element_filter) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = this->mesh.getSpatialDimension();
 
   auto types_iterable = mesh.elementTypes(spatial_dimension, ghost_type);
   if (element_filter) {
     types_iterable =
         element_filter->elementTypes(spatial_dimension, ghost_type);
   }
 
   for (auto type : types_iterable) {
     UInt nb_element = mesh.getNbElement(type, ghost_type);
     if (nb_element == 0)
       continue;
 
     const Array<UInt> * elem_filter;
     if (element_filter)
       elem_filter = &((*element_filter)(type, ghost_type));
     else
       elem_filter = &(empty_filter);
 
 #define AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS(type)                 \
   interpolateElementalFieldFromIntegrationPoints<type>(                        \
       field(type, ghost_type),                                                 \
       interpolation_points_coordinates_matrices(type, ghost_type),             \
       quad_points_coordinates_inv_matrices(type, ghost_type), result,          \
       ghost_type, *elem_filter)
 
     AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(
         AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS);
 #undef AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
diff --git a/src/fe_engine/shape_functions.hh b/src/fe_engine/shape_functions.hh
index ba7d2ef7b..a73dffcf4 100644
--- a/src/fe_engine/shape_functions.hh
+++ b/src/fe_engine/shape_functions.hh
@@ -1,212 +1,215 @@
 /**
  * @file   shape_functions.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  shape function class
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_memory.hh"
 #include "mesh.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_SHAPE_FUNCTIONS_HH__
 #define __AKANTU_SHAPE_FUNCTIONS_HH__
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 class ShapeFunctions : protected Memory {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
-  ShapeFunctions(const Mesh & mesh, const ID & id = "shape",
+  ShapeFunctions(const Mesh & mesh, UInt spatial_dimension, const ID & id = "shape",
                  const MemoryID & memory_id = 0);
   ~ShapeFunctions() override = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// function to print the contain of the class
   virtual void printself(std::ostream & stream, int indent = 0) const {
     std::string space;
     for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
       ;
     stream << space << "Shapes [" << std::endl;
     integration_points.printself(stream, indent + 1);
     // shapes.printself(stream, indent + 1);
     // shapes_derivatives.printself(stream, indent + 1);
     stream << space << "]" << std::endl;
   }
 
   /// set the integration points for a given element
   template <ElementType type>
   void setIntegrationPointsByType(const Matrix<Real> & integration_points,
                                   const GhostType & ghost_type);
 
   /// Build pre-computed matrices for interpolation of field form integration
   /// points at other given positions (interpolation_points)
   void initElementalFieldInterpolationFromIntegrationPoints(
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
       ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       const ElementTypeMapArray<Real> & quadrature_points_coordinates,
       const ElementTypeMapArray<UInt> * element_filter) const;
 
   /// Interpolate field at given position from given values of this field at
   /// integration points (field)
   /// using matrices precomputed with
   /// initElementalFieldInterplationFromIntegrationPoints
   void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> &
           interpolation_points_coordinates_matrices,
       const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       ElementTypeMapArray<Real> & result, const GhostType & ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const;
 
 protected:
   /// interpolate nodal values stored by element on the integration points
   template <ElementType type>
   void interpolateElementalFieldOnIntegrationPoints(
       const Array<Real> & u_el, Array<Real> & uq, const GhostType & ghost_type,
       const Array<Real> & shapes,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// gradient of nodal values stored by element on the control points
   template <ElementType type>
   void gradientElementalFieldOnIntegrationPoints(
       const Array<Real> & u_el, Array<Real> & out_nablauq,
       const GhostType & ghost_type, const Array<Real> & shapes_derivatives,
       const Array<UInt> & filter_elements) const;
 
 protected:
   /// By element versions of non-templated eponym methods
   template <ElementType type>
   inline void interpolateElementalFieldFromIntegrationPoints(
       const Array<Real> & field,
       const Array<Real> & interpolation_points_coordinates_matrices,
       const Array<Real> & quad_points_coordinates_inv_matrices,
       ElementTypeMapArray<Real> & result, const GhostType & ghost_type,
       const Array<UInt> & element_filter) const;
 
   /// Interpolate field at given position from given values of this field at
   /// integration points (field)
   /// using matrices precomputed with
   /// initElementalFieldInterplationFromIntegrationPoints
   template <ElementType type>
   inline void initElementalFieldInterpolationFromIntegrationPoints(
       const Array<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
       ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       const Array<Real> & quadrature_points_coordinates,
       const GhostType & ghost_type, const Array<UInt> & element_filter) const;
 
   /// build matrix for the interpolation of field form integration points
   template <ElementType type>
   inline void buildElementalFieldInterpolationMatrix(
       const Matrix<Real> & coordinates, Matrix<Real> & coordMatrix,
       UInt integration_order =
           ElementClassProperty<type>::polynomial_degree) const;
 
   /// build the so called interpolation matrix (first collumn is 1, then the
   /// other collumns are the traansposed coordinates)
   inline void buildInterpolationMatrix(const Matrix<Real> & coordinates,
                                        Matrix<Real> & coordMatrix,
                                        UInt integration_order) const;
 
 public:
   virtual void onElementsAdded(const Array<Element> &) {
     AKANTU_TO_IMPLEMENT();
   }
   virtual void onElementsRemoved(const Array<Element> &,
                                  const ElementTypeMapArray<UInt> &) {
     AKANTU_TO_IMPLEMENT();
   }
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the size of the shapes returned by the element class
   static inline UInt getShapeSize(const ElementType & type);
 
   /// get the size of the shapes derivatives returned by the element class
   static inline UInt getShapeDerivativesSize(const ElementType & type);
 
   inline const Matrix<Real> &
   getIntegrationPoints(const ElementType & type,
                        const GhostType & ghost_type) const {
     return integration_points(type, ghost_type);
   }
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// get a the shapes vector
   inline const Array<Real> &
   getShapes(const ElementType & el_type,
             const GhostType & ghost_type = _not_ghost) const;
 
   /// get a the shapes derivatives vector
   inline const Array<Real> &
   getShapesDerivatives(const ElementType & el_type,
                        const GhostType & ghost_type = _not_ghost) const;
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// shape functions for all elements
   ElementTypeMapArray<Real, InterpolationType> shapes;
 
   /// shape functions derivatives for all elements
   ElementTypeMapArray<Real, InterpolationType> shapes_derivatives;
 
   /// associated mesh
   const Mesh & mesh;
+  
+  // spatial dimension of the elements to consider
+  UInt _spatial_dimension;
 
   /// shape functions for all elements
   ElementTypeMap<Matrix<Real>> integration_points;
 };
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 /// standard output stream operator
 inline std::ostream & operator<<(std::ostream & stream,
                                  const ShapeFunctions & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 #include "shape_functions_inline_impl.cc"
 
 #endif /* __AKANTU_SHAPE_FUNCTIONS_HH__ */
diff --git a/src/fe_engine/shape_lagrange.hh b/src/fe_engine/shape_lagrange.hh
index 18a08b964..84aef3ff2 100644
--- a/src/fe_engine/shape_lagrange.hh
+++ b/src/fe_engine/shape_lagrange.hh
@@ -1,172 +1,173 @@
 /**
  * @file   shape_lagrange.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Feb 15 2011
  * @date last modification: Mon Jan 29 2018
  *
  * @brief  lagrangian shape functions class
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "shape_lagrange_base.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_SHAPE_LAGRANGE_HH__
 #define __AKANTU_SHAPE_LAGRANGE_HH__
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <class Shape> class ShapeCohesive;
 class ShapeIGFEM;
 
 template <ElementKind kind> class ShapeLagrange : public ShapeLagrangeBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
-  ShapeLagrange(const Mesh & mesh, const ID & id = "shape_lagrange",
+  ShapeLagrange(const Mesh & mesh, UInt spatial_dimension,
+                const ID & id = "shape_lagrange",
                 const MemoryID & memory_id = 0);
   ~ShapeLagrange() override = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// initialization function for structural elements not yet implemented
   inline void initShapeFunctions(const Array<Real> & nodes,
                                  const Matrix<Real> & integration_points,
                                  const ElementType & type,
                                  const GhostType & ghost_type);
 
   /// computes the shape functions derivatives for given interpolation points
   template <ElementType type>
   void computeShapeDerivativesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
       Array<Real> & shape_derivatives, const GhostType & ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   void computeShapeDerivativesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
       Array<Real> & shape_derivatives, const ElementType & type,
       const GhostType & ghost_type,
       const Array<UInt> & filter_elements) const override;
 
   /// pre compute all shapes on the element integration points from natural
   /// coordinates
   template <ElementType type>
   void precomputeShapesOnIntegrationPoints(const Array<Real> & nodes,
                                            const GhostType & ghost_type);
 
   /// pre compute all shape derivatives on the element integration points from
   /// natural coordinates
   template <ElementType type>
   void
   precomputeShapeDerivativesOnIntegrationPoints(const Array<Real> & nodes,
                                                 const GhostType & ghost_type);
 
   /// interpolate nodal values on the integration points
   template <ElementType type>
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   template <ElementType type>
   void interpolateOnIntegrationPoints(
       const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
       const Array<Real> & shapes, GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// interpolate on physical point
   template <ElementType type>
   void interpolate(const Vector<Real> & real_coords, UInt elem,
                    const Matrix<Real> & nodal_values,
                    Vector<Real> & interpolated,
                    const GhostType & ghost_type) const;
 
   /// compute the gradient of u on the integration points
   template <ElementType type>
   void gradientOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   template <ElementType type>
   void computeBtD(const Array<Real> & Ds, Array<Real> & BtDs,
                   GhostType ghost_type,
                   const Array<UInt> & filter_elements) const;
 
   template <ElementType type>
   void computeBtDB(const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
                    GhostType ghost_type,
                    const Array<UInt> & filter_elements) const;
 
   /// multiply a field by shape functions  @f$ fts_{ij} = f_i * \varphi_j @f$
   template <ElementType type>
   void computeNtb(const Array<Real> & bs, Array<Real> & Ntbs,
                   GhostType ghost_type,
                   const Array<UInt> & filter_elements = empty_filter) const;
 
   /// find natural coords from real coords provided an element
   template <ElementType type>
   void inverseMap(const Vector<Real> & real_coords, UInt element,
                   Vector<Real> & natural_coords,
                   const GhostType & ghost_type = _not_ghost) const;
 
   /// return true if the coordinates provided are inside the element, false
   /// otherwise
   template <ElementType type>
   bool contains(const Vector<Real> & real_coords, UInt elem,
                 const GhostType & ghost_type) const;
 
   /// compute the shape on a provided point
   template <ElementType type>
   void computeShapes(const Vector<Real> & real_coords, UInt elem,
                      Vector<Real> & shapes, const GhostType & ghost_type) const;
 
   /// compute the shape derivatives on a provided point
   template <ElementType type>
   void computeShapeDerivatives(const Matrix<Real> & real_coords, UInt elem,
                                Tensor3<Real> & shapes,
                                const GhostType & ghost_type) const;
 
 protected:
   /// compute the shape derivatives on integration points for a given element
   template <ElementType type>
   inline void
   computeShapeDerivativesOnCPointsByElement(const Matrix<Real> & node_coords,
                                             const Matrix<Real> & natural_coords,
                                             Tensor3<Real> & shapesd) const;
 };
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "shape_lagrange_inline_impl.cc"
 
 #endif /* __AKANTU_SHAPE_LAGRANGE_HH__ */
diff --git a/src/fe_engine/shape_lagrange_base.cc b/src/fe_engine/shape_lagrange_base.cc
index 71d0f3942..07a32a9b3 100644
--- a/src/fe_engine/shape_lagrange_base.cc
+++ b/src/fe_engine/shape_lagrange_base.cc
@@ -1,160 +1,164 @@
 /**
  * @file   shape_lagrange_base.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Aug 09 2017
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  common par for the shape lagrange
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "shape_lagrange_base.hh"
 #include "mesh_iterators.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 ShapeLagrangeBase::ShapeLagrangeBase(const Mesh & mesh,
+                                     UInt spatial_dimension,
                                      const ElementKind & kind, const ID & id,
                                      const MemoryID & memory_id)
-    : ShapeFunctions(mesh, id, memory_id), _kind(kind) {}
+    : ShapeFunctions(mesh, spatial_dimension, id, memory_id), _kind(kind) {}
 
 /* -------------------------------------------------------------------------- */
 ShapeLagrangeBase::~ShapeLagrangeBase() = default;
 
 /* -------------------------------------------------------------------------- */
 #define AKANTU_COMPUTE_SHAPES(type)                                            \
   _this.template computeShapesOnIntegrationPoints<type>(                       \
       nodes, integration_points, shapes, ghost_type, filter_elements)
 
 namespace shape_lagrange {
   namespace details {
     template <ElementKind kind> struct Helper {
       template <class S>
       static void call(const S &, const Array<Real> &, const Matrix<Real> &,
                        Array<Real> &, const ElementType &, const GhostType &,
                        const Array<UInt> &) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define AKANTU_COMPUTE_SHAPES_KIND(kind)                                       \
   template <> struct Helper<kind> {                                            \
     template <class S>                                                         \
     static void call(const S & _this, const Array<Real> & nodes,               \
                      const Matrix<Real> & integration_points,                  \
                      Array<Real> & shapes, const ElementType & type,           \
                      const GhostType & ghost_type,                             \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES, kind);           \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_COMPUTE_SHAPES_KIND,
                                AKANTU_FE_ENGINE_LIST_LAGRANGE_BASE)
 
   } // namespace details
 } // namespace shape_lagrange
 
 /* -------------------------------------------------------------------------- */
 void ShapeLagrangeBase::computeShapesOnIntegrationPoints(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     Array<Real> & shapes, const ElementType & type,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
 
   auto kind = Mesh::getKind(type);
 
 #define AKANTU_COMPUTE_SHAPES_KIND_SWITCH(kind)                                \
   shape_lagrange::details::Helper<kind>::call(                                 \
       *this, nodes, integration_points, shapes, type, ghost_type,              \
       filter_elements);
 
   AKANTU_BOOST_LIST_SWITCH(
       AKANTU_COMPUTE_SHAPES_KIND_SWITCH,
       BOOST_PP_LIST_TO_SEQ(AKANTU_FE_ENGINE_LIST_LAGRANGE_BASE), kind);
 
 #undef AKANTU_COMPUTE_SHAPES
 #undef AKANTU_COMPUTE_SHAPES_KIND
 #undef AKANTU_COMPUTE_SHAPES_KIND_SWITCH
 }
 
 /* -------------------------------------------------------------------------- */
 void ShapeLagrangeBase::onElementsAdded(const Array<Element> & new_elements) {
   AKANTU_DEBUG_IN();
   const auto & nodes = mesh.getNodes();
 
   for (auto elements_range : MeshElementsByTypes(new_elements)) {
     auto type = elements_range.getType();
     auto ghost_type = elements_range.getGhostType();
 
+    if(mesh.getSpatialDimension(type) != _spatial_dimension)
+      continue;
+    
     if (mesh.getKind(type) != _kind)
       continue;
 
     auto & elements = elements_range.getElements();
 
     auto itp_type = FEEngine::getInterpolationType(type);
 
     if (not this->shapes_derivatives.exists(itp_type, ghost_type)) {
       auto size_of_shapesd = this->getShapeDerivativesSize(type);
       this->shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type);
     }
 
     if (not shapes.exists(itp_type, ghost_type)) {
       auto size_of_shapes = this->getShapeSize(type);
       this->shapes.alloc(0, size_of_shapes, itp_type, ghost_type);
     }
 
     const auto & natural_coords = integration_points(type, ghost_type);
     computeShapesOnIntegrationPoints(nodes, natural_coords,
                                      shapes(itp_type, ghost_type), type,
                                      ghost_type, elements);
 
     computeShapeDerivativesOnIntegrationPoints(
         nodes, natural_coords, shapes_derivatives(itp_type, ghost_type), type,
         ghost_type, elements);
   }
 #undef INIT_SHAPE_FUNCTIONS
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void ShapeLagrangeBase::onElementsRemoved(
     const Array<Element> &, const ElementTypeMapArray<UInt> & new_numbering) {
   this->shapes.onElementsRemoved(new_numbering);
   this->shapes_derivatives.onElementsRemoved(new_numbering);
 }
 
 /* -------------------------------------------------------------------------- */
 void ShapeLagrangeBase::printself(std::ostream & stream, int indent) const {
   std::string space(indent, AKANTU_INDENT);
 
   stream << space << "Shapes Lagrange [" << std::endl;
   ShapeFunctions::printself(stream, indent + 1);
   shapes.printself(stream, indent + 1);
   shapes_derivatives.printself(stream, indent + 1);
   stream << space << "]" << std::endl;
 }
 
 } // namespace akantu
diff --git a/src/fe_engine/shape_lagrange_base.hh b/src/fe_engine/shape_lagrange_base.hh
index c50b27689..39dae5598 100644
--- a/src/fe_engine/shape_lagrange_base.hh
+++ b/src/fe_engine/shape_lagrange_base.hh
@@ -1,92 +1,92 @@
 /**
  * @file   shape_lagrange_base.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Aug 09 2017
  * @date last modification: Wed Nov 08 2017
  *
  * @brief  Base class for the shape lagrange
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "shape_functions.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_SHAPE_LAGRANGE_BASE_HH__
 #define __AKANTU_SHAPE_LAGRANGE_BASE_HH__
 
 namespace akantu {
 
 class ShapeLagrangeBase : public ShapeFunctions {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
-  ShapeLagrangeBase(const Mesh & mesh, const ElementKind & kind,
-                    const ID & id = "shape_lagrange",
+  ShapeLagrangeBase(const Mesh & mesh, UInt spatial_dimension,
+                    const ElementKind & kind, const ID & id = "shape_lagrange",
                     const MemoryID & memory_id = 0);
   ~ShapeLagrangeBase() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// computes the shape functions for given interpolation points
   virtual void computeShapesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
       Array<Real> & shapes, const ElementType & type,
       const GhostType & ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// computes the shape functions derivatives for given interpolation points
   virtual void computeShapeDerivativesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
       Array<Real> & shape_derivatives, const ElementType & type,
       const GhostType & ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// function to print the containt of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   template <ElementType type>
   void computeShapesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
       Array<Real> & shapes, const GhostType & ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const;
 
 public:
   void onElementsAdded(const Array<Element> & elements) override;
   void
   onElementsRemoved(const Array<Element> & elements,
                     const ElementTypeMapArray<UInt> & new_numbering) override;
 
 protected:
   /// The kind to consider
   ElementKind _kind;
 };
 
 } // namespace akantu
 
 #include "shape_lagrange_base_inline_impl.cc"
 
 #endif /* __AKANTU_SHAPE_LAGRANGE_BASE_HH__ */
diff --git a/src/fe_engine/shape_lagrange_inline_impl.cc b/src/fe_engine/shape_lagrange_inline_impl.cc
index d26190221..817fde29c 100644
--- a/src/fe_engine/shape_lagrange_inline_impl.cc
+++ b/src/fe_engine/shape_lagrange_inline_impl.cc
@@ -1,538 +1,538 @@
 /**
  * @file   shape_lagrange_inline_impl.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Oct 27 2010
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  ShapeLagrange inline implementation
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu 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.
  *
  * Akantu 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  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_iterators.hh"
 #include "aka_voigthelper.hh"
 #include "fe_engine.hh"
 #include "shape_lagrange.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__
 #define __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 #define INIT_SHAPE_FUNCTIONS(type)                                             \
   setIntegrationPointsByType<type>(integration_points, ghost_type);            \
   precomputeShapesOnIntegrationPoints<type>(nodes, ghost_type);                \
   if (ElementClass<type>::getNaturalSpaceDimension() ==                        \
           mesh.getSpatialDimension() ||                                        \
       kind != _ek_regular)                                                     \
     precomputeShapeDerivativesOnIntegrationPoints<type>(nodes, ghost_type);
 
 template <ElementKind kind>
 inline void ShapeLagrange<kind>::initShapeFunctions(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     const ElementType & type, const GhostType & ghost_type) {
   AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
 }
 
 #undef INIT_SHAPE_FUNCTIONS
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 inline void ShapeLagrange<kind>::computeShapeDerivativesOnCPointsByElement(
     const Matrix<Real> & node_coords, const Matrix<Real> & natural_coords,
     Tensor3<Real> & shapesd) const {
   AKANTU_DEBUG_IN();
 
   // compute dnds
   Tensor3<Real> dnds(node_coords.rows(), node_coords.cols(),
                      natural_coords.cols());
   ElementClass<type>::computeDNDS(natural_coords, dnds);
   // compute jacobian
   Tensor3<Real> J(node_coords.rows(), natural_coords.rows(),
                   natural_coords.cols());
   ElementClass<type>::computeJMat(dnds, node_coords, J);
 
   // compute dndx
   ElementClass<type>::computeShapeDerivatives(J, dnds, shapesd);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::inverseMap(const Vector<Real> & real_coords,
                                      UInt elem, Vector<Real> & natural_coords,
                                      const GhostType & ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage();
   Matrix<Real> nodes_coord(spatial_dimension, nb_nodes_per_element);
 
   mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(),
                                      elem_val + elem * nb_nodes_per_element,
                                      nb_nodes_per_element, spatial_dimension);
 
   ElementClass<type>::inverseMap(real_coords, nodes_coord, natural_coords);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 bool ShapeLagrange<kind>::contains(const Vector<Real> & real_coords, UInt elem,
                                    const GhostType & ghost_type) const {
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   Vector<Real> natural_coords(spatial_dimension);
 
   inverseMap<type>(real_coords, elem, natural_coords, ghost_type);
   return ElementClass<type>::contains(natural_coords);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::interpolate(const Vector<Real> & real_coords,
                                       UInt elem,
                                       const Matrix<Real> & nodal_values,
                                       Vector<Real> & interpolated,
                                       const GhostType & ghost_type) const {
   UInt nb_shapes = ElementClass<type>::getShapeSize();
   Vector<Real> shapes(nb_shapes);
   computeShapes<type>(real_coords, elem, shapes, ghost_type);
   ElementClass<type>::interpolate(nodal_values, shapes, interpolated);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeShapes(const Vector<Real> & real_coords,
                                         UInt elem, Vector<Real> & shapes,
                                         const GhostType & ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   Vector<Real> natural_coords(spatial_dimension);
 
   inverseMap<type>(real_coords, elem, natural_coords, ghost_type);
   ElementClass<type>::computeShapes(natural_coords, shapes);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeShapeDerivatives(
     const Matrix<Real> & real_coords, UInt elem, Tensor3<Real> & shapesd,
     const GhostType & ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_points = real_coords.cols();
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == shapesd.size(0) &&
                           nb_nodes_per_element == shapesd.size(1),
                       "Shape size doesn't match");
   AKANTU_DEBUG_ASSERT(nb_points == shapesd.size(2),
                       "Number of points doesn't match shapes size");
 
   Matrix<Real> natural_coords(spatial_dimension, nb_points);
 
   // Creates the matrix of natural coordinates
   for (UInt i = 0; i < nb_points; i++) {
     Vector<Real> real_point = real_coords(i);
     Vector<Real> natural_point = natural_coords(i);
 
     inverseMap<type>(real_point, elem, natural_point, ghost_type);
   }
 
   UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage();
   Matrix<Real> nodes_coord(spatial_dimension, nb_nodes_per_element);
 
   mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(),
                                      elem_val + elem * nb_nodes_per_element,
                                      nb_nodes_per_element, spatial_dimension);
 
   computeShapeDerivativesOnCPointsByElement<type>(nodes_coord, natural_coords,
                                                   shapesd);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
-ShapeLagrange<kind>::ShapeLagrange(const Mesh & mesh, const ID & id,
-                                   const MemoryID & memory_id)
-    : ShapeLagrangeBase(mesh, kind, id, memory_id) {}
+ShapeLagrange<kind>::ShapeLagrange(const Mesh & mesh, UInt spatial_dimension,
+                                   const ID & id, const MemoryID & memory_id)
+    : ShapeLagrangeBase(mesh, spatial_dimension, kind, id, memory_id) {}
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     Array<Real> & shape_derivatives, const GhostType & ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
 
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   UInt nb_points = integration_points.cols();
   UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
 
   UInt size_of_shapesd = ElementClass<type>::getShapeDerivativesSize();
   AKANTU_DEBUG_ASSERT(shape_derivatives.getNbComponent() == size_of_shapesd,
                       "The shapes_derivatives array does not have the correct "
                           << "number of component");
   shape_derivatives.resize(nb_element * nb_points);
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type,
                                        filter_elements);
 
   Real * shapesd_val = shape_derivatives.storage();
   Array<Real>::matrix_iterator x_it =
       x_el.begin(spatial_dimension, nb_nodes_per_element);
 
   if (filter_elements != empty_filter)
     nb_element = filter_elements.size();
 
   for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) {
     if (filter_elements != empty_filter)
       shapesd_val = shape_derivatives.storage() +
                     filter_elements(elem) * size_of_shapesd * nb_points;
 
     Matrix<Real> & X = *x_it;
     Tensor3<Real> B(shapesd_val, spatial_dimension, nb_nodes_per_element,
                     nb_points);
     computeShapeDerivativesOnCPointsByElement<type>(X, integration_points, B);
 
     if (filter_elements == empty_filter)
       shapesd_val += size_of_shapesd * nb_points;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 void ShapeLagrange<kind>::computeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     Array<Real> & shape_derivatives, const ElementType & type,
     const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
 #define AKANTU_COMPUTE_SHAPES(type)                                            \
   computeShapeDerivativesOnIntegrationPoints<type>(                            \
       nodes, integration_points, shape_derivatives, ghost_type,                \
       filter_elements);
 
   AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES);
 
 #undef AKANTU_COMPUTE_SHAPES
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::precomputeShapesOnIntegrationPoints(
     const Array<Real> & nodes, const GhostType & ghost_type) {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   Matrix<Real> & natural_coords = integration_points(type, ghost_type);
   UInt size_of_shapes = ElementClass<type>::getShapeSize();
 
   Array<Real> & shapes_tmp =
       shapes.alloc(0, size_of_shapes, itp_type, ghost_type);
 
   this->computeShapesOnIntegrationPoints<type>(nodes, natural_coords,
                                                shapes_tmp, ghost_type);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::precomputeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, const GhostType & ghost_type) {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   Matrix<Real> & natural_coords = integration_points(type, ghost_type);
   UInt size_of_shapesd = ElementClass<type>::getShapeDerivativesSize();
 
   Array<Real> & shapes_derivatives_tmp =
       shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type);
 
   this->computeShapeDerivativesOnIntegrationPoints<type>(
       nodes, natural_coords, shapes_derivatives_tmp, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::interpolateOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
     const Array<Real> & shapes, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   Array<Real> u_el(0, nb_degree_of_freedom * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
                                        filter_elements);
 
   this->interpolateElementalFieldOnIntegrationPoints<type>(
       u_el, out_uq, ghost_type, shapes, filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::interpolateOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   AKANTU_DEBUG_ASSERT(shapes.exists(itp_type, ghost_type),
                       "No shapes for the type "
                           << shapes.printType(itp_type, ghost_type));
 
   this->interpolateOnIntegrationPoints<type>(in_u, out_uq, nb_degree_of_freedom,
                                              shapes(itp_type, ghost_type),
                                              ghost_type, filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::gradientOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_nablauq,
     UInt nb_degree_of_freedom, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   AKANTU_DEBUG_ASSERT(
       shapes_derivatives.exists(itp_type, ghost_type),
       "No shapes derivatives for the type "
           << shapes_derivatives.printType(itp_type, ghost_type));
 
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   Array<Real> u_el(0, nb_degree_of_freedom * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
                                        filter_elements);
 
   this->gradientElementalFieldOnIntegrationPoints<type>(
       u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type),
       filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeBtD(
     const Array<Real> & Ds, Array<Real> & BtDs, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   auto itp_type = ElementClassProperty<type>::interpolation_type;
   const auto & shapes_derivatives =
       this->shapes_derivatives(itp_type, ghost_type);
 
   auto spatial_dimension = mesh.getSpatialDimension();
   auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
 
   Array<Real> shapes_derivatives_filtered(0,
                                           shapes_derivatives.getNbComponent());
   auto && view =
       make_view(shapes_derivatives, spatial_dimension, nb_nodes_per_element);
   auto B_it = view.begin();
   auto B_end = view.end();
 
   if (filter_elements != empty_filter) {
     FEEngine::filterElementalData(this->mesh, shapes_derivatives,
                                   shapes_derivatives_filtered, type, ghost_type,
                                   filter_elements);
     auto && view = make_view(shapes_derivatives_filtered, spatial_dimension,
                              nb_nodes_per_element);
     B_it = view.begin();
     B_end = view.end();
   }
 
   for (auto && values :
        zip(range(B_it, B_end),
            make_view(Ds, Ds.getNbComponent() / spatial_dimension,
                      spatial_dimension),
            make_view(BtDs, BtDs.getNbComponent() / nb_nodes_per_element,
                      nb_nodes_per_element))) {
     const auto & B = std::get<0>(values);
     const auto & D = std::get<1>(values);
     auto & Bt_D = std::get<2>(values);
     // transposed due to the storage layout of B
     Bt_D.template mul<false, false>(D, B);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeBtDB(
     const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   auto itp_type = ElementClassProperty<type>::interpolation_type;
   const auto & shapes_derivatives =
       this->shapes_derivatives(itp_type, ghost_type);
 
   constexpr auto dim = ElementClass<type>::getSpatialDimension();
   auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
 
   Array<Real> shapes_derivatives_filtered(0,
                                           shapes_derivatives.getNbComponent());
   auto && view = make_view(shapes_derivatives, dim, nb_nodes_per_element);
   auto B_it = view.begin();
   auto B_end = view.end();
 
   if (filter_elements != empty_filter) {
     FEEngine::filterElementalData(this->mesh, shapes_derivatives,
                                   shapes_derivatives_filtered, type, ghost_type,
                                   filter_elements);
     auto && view =
         make_view(shapes_derivatives_filtered, dim, nb_nodes_per_element);
     B_it = view.begin();
     B_end = view.end();
   }
 
   if (order_d == 4) {
     UInt tangent_size = VoigtHelper<dim>::size;
     Matrix<Real> B(tangent_size, dim * nb_nodes_per_element);
     Matrix<Real> Bt_D(dim * nb_nodes_per_element, tangent_size);
 
     for (auto && values :
          zip(range(B_it, B_end), make_view(Ds, tangent_size, tangent_size),
              make_view(BtDBs, dim * nb_nodes_per_element,
                        dim * nb_nodes_per_element))) {
       const auto & Bfull = std::get<0>(values);
       const auto & D = std::get<1>(values);
       auto & Bt_D_B = std::get<2>(values);
 
       VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(Bfull, B,
                                                          nb_nodes_per_element);
       Bt_D.template mul<true, false>(B, D);
       Bt_D_B.template mul<false, false>(Bt_D, B);
     }
   } else if (order_d == 2) {
     Matrix<Real> Bt_D(nb_nodes_per_element, dim);
     for (auto && values :
          zip(range(B_it, B_end), make_view(Ds, dim, dim),
              make_view(BtDBs, nb_nodes_per_element, nb_nodes_per_element))) {
       const auto & B = std::get<0>(values);
       const auto & D = std::get<1>(values);
       auto & Bt_D_B = std::get<2>(values);
       Bt_D.template mul<true, false>(B, D);
       Bt_D_B.template mul<false, false>(Bt_D, B);
     }
   }
 }
 
 template <>
 template <>
 inline void ShapeLagrange<_ek_regular>::computeBtDB<_point_1>(
     const Array<Real> & /*Ds*/, Array<Real> & /*BtDBs*/, UInt /*order_d*/,
     GhostType /*ghost_type*/, const Array<UInt> & /*filter_elements*/) const {
   AKANTU_TO_IMPLEMENT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeNtb(
     const Array<Real> & bs, Array<Real> & Ntbs, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   Ntbs.resize(bs.size());
 
   UInt size_of_shapes = ElementClass<type>::getShapeSize();
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   UInt nb_degree_of_freedom = bs.getNbComponent();
 
   Array<Real> shapes_filtered(0, size_of_shapes);
   auto && view = make_view(shapes(itp_type, ghost_type), 1, size_of_shapes);
   auto N_it = view.begin();
   auto N_end = view.end();
 
   if (filter_elements != empty_filter) {
     FEEngine::filterElementalData(this->mesh, shapes(itp_type, ghost_type),
                                   shapes_filtered, type, ghost_type,
                                   filter_elements);
     auto && view = make_view(shapes_filtered, 1, size_of_shapes);
     N_it = view.begin();
     N_end = view.end();
   }
 
   for (auto && values :
        zip(make_view(bs, nb_degree_of_freedom, 1), range(N_it, N_end),
            make_view(Ntbs, nb_degree_of_freedom, size_of_shapes))) {
     const auto & b = std::get<0>(values);
     const auto & N = std::get<1>(values);
     auto & Ntb = std::get<2>(values);
 
     Ntb.template mul<false, false>(b, N);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
 
 #endif /* __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__ */