diff --git a/src/fe_engine/element_class_structural.hh b/src/fe_engine/element_class_structural.hh
index fca7db751..76732da6a 100644
--- a/src/fe_engine/element_class_structural.hh
+++ b/src/fe_engine/element_class_structural.hh
@@ -1,279 +1,274 @@
 /**
  * @file   element_class_structural.hh
  *
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Damien Spielmann <damien.spielmann@epfl.ch>
  *
  * @date creation: Thu Feb 21 2013
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Specialization of the element classes for structural elements
  *
  *
  * Copyright (©) 2014-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 "element_class.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_ELEMENT_CLASS_STRUCTURAL_HH_
 #define AKANTU_ELEMENT_CLASS_STRUCTURAL_HH_
 
 namespace akantu {
 
 /// Macro to generate the InterpolationProperty structures for different
 /// interpolation types
 #define AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(                  \
     itp_type, itp_geom_type, ndof, nb_stress, nb_dnds_cols)                    \
   template <> struct InterpolationProperty<itp_type> {                         \
     static const InterpolationKind kind{_itk_structural};                      \
     static const UInt nb_nodes_per_element{                                    \
         InterpolationProperty<itp_geom_type>::nb_nodes_per_element};           \
     static const InterpolationType itp_geometry_type{itp_geom_type};           \
     static const UInt natural_space_dimension{                                 \
         InterpolationProperty<itp_geom_type>::natural_space_dimension};        \
     static const UInt nb_degree_of_freedom{ndof};                              \
     static const UInt nb_stress_components{nb_stress};                         \
     static const UInt dnds_columns{nb_dnds_cols};                              \
   }
 
 /* -------------------------------------------------------------------------- */
 template <InterpolationType interpolation_type>
 class InterpolationElement<interpolation_type, _itk_structural> {
 public:
   using interpolation_property = InterpolationProperty<interpolation_type>;
 
   /// compute the shape values for a given set of points in natural coordinates
   static inline void computeShapes(const Matrix<Real> & natural_coord,
                                    const Matrix<Real> & real_coord,
-                                   Tensor3<Real> & N) {
+                                   const Matrix<Real> & T, Tensor3<Real> & Ns) {
     for (UInt i = 0; i < natural_coord.cols(); ++i) {
-      Matrix<Real> n_t = N(i);
-      computeShapes(natural_coord(i), real_coord, n_t);
+      Matrix<Real> N_T = Ns(i);
+      Matrix<Real> N(N_T.rows(), N_T.cols());
+      computeShapes(natural_coord(i), real_coord, N);
+      N_T.mul<false, false>(N, T);
     }
   }
 
   /// compute the shape values for a given point in natural coordinates
   static inline void computeShapes(const Vector<Real> & natural_coord,
                                    const Matrix<Real> & real_coord,
                                    Matrix<Real> & N);
 
-  static inline void computeShapesMass(const Matrix<Real> & natural_coord,
-                                       const Matrix<Real> & real_coord,
-                                       Tensor3<Real> & N) {
-    for (UInt i = 0; i < natural_coord.cols(); ++i) {
-      Matrix<Real> n_t = N(i);
-      computeShapesMass(natural_coord(i), real_coord, n_t);
-    }
-  }
-
-  static inline void computeShapesMass(const Vector<Real> & natural_coords,
-                                       const Matrix<Real> & real_coords,
-                                       Matrix<Real> & N) {
-    Matrix<Real> Ntotal(interpolation_property::nb_degree_of_freedom,
-                        interpolation_property::nb_degree_of_freedom *
-                            interpolation_property::nb_nodes_per_element);
-    computeShapes(natural_coords, real_coords, Ntotal);
+  static inline void computeShapesMass(const Matrix<Real> & natural_coords,
+                                       const Matrix<Real> & xs,
+                                       const Matrix<Real> & T,
+                                       Tensor3<Real> & Ns) {
+    for (UInt i = 0; i < natural_coords.cols(); ++i) {
+      Matrix<Real> N_T = Ns(i);
+      Vector<Real> X = natural_coords(i);
+      Matrix<Real> N(interpolation_property::nb_degree_of_freedom, N_T.cols());
 
-    N = Ntotal.block(0, 0, N.rows(), N.cols());
+      computeShapes(X, xs, N);
+      N_T.mul<false, false>(N.block(0, 0, N_T.rows(), N_T.cols()), T);
+    }
   }
 
   /// compute shape derivatives (input is dxds) for a set of points
   static inline void computeShapeDerivatives(const Tensor3<Real> & Js,
                                              const Tensor3<Real> & DNDSs,
                                              const Matrix<Real> & R,
                                              Tensor3<Real> & Bs) {
     for (UInt i = 0; i < Js.size(2); ++i) {
       Matrix<Real> J = Js(i);
       Matrix<Real> DNDS = DNDSs(i);
       Matrix<Real> DNDX(DNDS.rows(), DNDS.cols());
       auto inv_J = J.inverse();
       DNDX.mul<false, false>(inv_J, DNDS);
       Matrix<Real> B_R = Bs(i);
       Matrix<Real> B(B_R.rows(), B_R.cols());
       arrangeInVoigt(DNDX, B);
       B_R.mul<false, false>(B, R);
     }
   }
 
   /**
    * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of
    * shape functions along with variation of natural coordinates on a given set
    * of points in natural coordinates
    */
   static inline void computeDNDS(const Matrix<Real> & natural_coord,
                                  const Matrix<Real> & real_coord,
                                  Tensor3<Real> & dnds) {
     for (UInt i = 0; i < natural_coord.cols(); ++i) {
       Matrix<Real> dnds_t = dnds(i);
       computeDNDS(natural_coord(i), real_coord, dnds_t);
     }
   }
 
   /**
    * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of
    * shape functions along with
    * variation of natural coordinates on a given point in natural
    * coordinates
    */
   static inline void computeDNDS(const Vector<Real> & natural_coord,
                                  const Matrix<Real> & real_coord,
                                  Matrix<Real> & dnds);
 
   /**
    * arrange B in Voigt notation from DNDS
    */
   static inline void arrangeInVoigt(const Matrix<Real> & dnds,
                                     Matrix<Real> & B) {
     // Default implementation assumes dnds is already in Voigt notation
     B.deepCopy(dnds);
   }
 
 public:
   static inline constexpr auto getNbNodesPerInterpolationElement() {
     return interpolation_property::nb_nodes_per_element;
   }
 
   static inline constexpr auto getShapeSize() {
     return interpolation_property::nb_nodes_per_element *
            interpolation_property::nb_degree_of_freedom *
            interpolation_property::nb_degree_of_freedom;
   }
+  static inline constexpr auto getShapeIndependantSize() {
+    return interpolation_property::nb_nodes_per_element *
+           interpolation_property::nb_degree_of_freedom *
+           interpolation_property::nb_stress_components;
+  }
   static inline constexpr auto getShapeDerivativesSize() {
     return interpolation_property::nb_nodes_per_element *
            interpolation_property::nb_degree_of_freedom *
            interpolation_property::nb_stress_components;
   }
   static inline constexpr auto getNaturalSpaceDimension() {
     return interpolation_property::natural_space_dimension;
   }
   static inline constexpr auto getNbDegreeOfFreedom() {
     return interpolation_property::nb_degree_of_freedom;
   }
   static inline constexpr auto getNbStressComponents() {
     return interpolation_property::nb_stress_components;
   }
-
-  static inline constexpr auto getShapeMassSize() {
-    return interpolation_property::natural_space_dimension *
-           interpolation_property::nb_degree_of_freedom *
-           interpolation_property::nb_nodes_per_element;
-  }
 };
 
 /// Macro to generate the element class structures for different structural
 /// element types
 /* -------------------------------------------------------------------------- */
 #define AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY(                       \
     elem_type, geom_type, interp_type, parent_el_type, elem_kind, sp,          \
     gauss_int_type, min_int_order)                                             \
   template <> struct ElementClassProperty<elem_type> {                         \
     static const GeometricalType geometrical_type{geom_type};                  \
     static const InterpolationType interpolation_type{interp_type};            \
     static const ElementType parent_element_type{parent_el_type};              \
     static const ElementKind element_kind{elem_kind};                          \
     static const UInt spatial_dimension{sp};                                   \
     static const GaussIntegrationType gauss_integration_type{gauss_int_type};  \
     static const UInt polynomial_degree{min_int_order};                        \
   }
 
 /* -------------------------------------------------------------------------- */
 /* ElementClass for structural elements                                       */
 /* -------------------------------------------------------------------------- */
 template <ElementType element_type>
 class ElementClass<element_type, _ek_structural>
     : public GeometricalElement<
           ElementClassProperty<element_type>::geometrical_type>,
       public InterpolationElement<
           ElementClassProperty<element_type>::interpolation_type> {
 protected:
   using geometrical_element =
       GeometricalElement<ElementClassProperty<element_type>::geometrical_type>;
   using interpolation_element = InterpolationElement<
       ElementClassProperty<element_type>::interpolation_type>;
   using parent_element =
       ElementClass<ElementClassProperty<element_type>::parent_element_type>;
 
 public:
   static inline void
   computeRotationMatrix(Matrix<Real> & /*R*/, const Matrix<Real> & /*X*/,
                         const Vector<Real> & /*extra_normal*/) {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// compute jacobian (or integration variable change factor) for a given point
   static inline void computeJMat(const Vector<Real> & natural_coords,
                                  const Matrix<Real> & Xs, Matrix<Real> & J) {
     Matrix<Real> dnds(Xs.rows(), Xs.cols());
     parent_element::computeDNDS(natural_coords, dnds);
     J.mul<false, true>(dnds, Xs);
   }
 
   static inline void computeJMat(const Matrix<Real> & natural_coords,
                                  const Matrix<Real> & Xs, Tensor3<Real> & Js) {
     for (UInt i = 0; i < natural_coords.cols(); ++i) {
       // because non-const l-value reference does not bind to r-value
       Matrix<Real> J = Js(i);
       computeJMat(Vector<Real>(natural_coords(i)), Xs, J);
     }
   }
 
   static inline void computeJacobian(const Matrix<Real> & natural_coords,
                                      const Matrix<Real> & node_coords,
                                      Vector<Real> & jacobians) {
     using itp = typename interpolation_element::interpolation_property;
     Tensor3<Real> Js(itp::natural_space_dimension, itp::natural_space_dimension,
                      natural_coords.cols());
     computeJMat(natural_coords, node_coords, Js);
     for (UInt i = 0; i < natural_coords.cols(); ++i) {
       Matrix<Real> J = Js(i);
       jacobians(i) = J.det();
     }
   }
 
   static inline void computeRotation(const Matrix<Real> & node_coords,
                                      Matrix<Real> & rotation);
 
 public:
   static AKANTU_GET_MACRO_NOT_CONST(Kind, _ek_structural, ElementKind);
   static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, _not_defined, ElementType);
   static AKANTU_GET_MACRO_NOT_CONST(FacetType, _not_defined, ElementType);
   static constexpr auto getFacetType(__attribute__((unused)) UInt t = 0) {
     return _not_defined;
   }
   static constexpr AKANTU_GET_MACRO_NOT_CONST(
       SpatialDimension, ElementClassProperty<element_type>::spatial_dimension,
       UInt);
   static constexpr auto getFacetTypes() {
     return ElementClass<_not_defined>::getFacetTypes();
   }
 };
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 #include "element_class_hermite_inline_impl.hh"
 /* keep order */
 #include "element_class_bernoulli_beam_inline_impl.hh"
 #include "element_class_kirchhoff_shell_inline_impl.hh"
 /* -------------------------------------------------------------------------- */
 
 #endif /* AKANTU_ELEMENT_CLASS_STRUCTURAL_HH_ */
diff --git a/src/fe_engine/fe_engine_template_tmpl_field.hh b/src/fe_engine/fe_engine_template_tmpl_field.hh
index 4500396b2..bd1b4d5fb 100644
--- a/src/fe_engine/fe_engine_template_tmpl_field.hh
+++ b/src/fe_engine/fe_engine_template_tmpl_field.hh
@@ -1,500 +1,503 @@
 /**
  * @file   fe_engine_template_tmpl_field.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Aug 09 2017
  * @date last modification: Thu Dec 07 2017
  *
  * @brief  implementation of the assemble field s functions
  *
  *
  * 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 "fe_engine_template.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH_
 #define AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH_
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /* Matrix lumping functions                                                   */
 /* -------------------------------------------------------------------------- */
 namespace fe_engine {
   namespace details {
     namespace {
       template <class Functor>
       void fillField(const Functor & field_funct, Array<Real> & field,
                      UInt nb_element, UInt nb_integration_points,
                      ElementType type, GhostType ghost_type) {
         UInt nb_degree_of_freedom = field.getNbComponent();
         field.resize(nb_integration_points * nb_element);
 
         auto field_it = field.begin_reinterpret(
             nb_degree_of_freedom, nb_integration_points, nb_element);
 
         Element el{type, 0, ghost_type};
         for (; el.element < nb_element; ++el.element, ++field_it) {
           field_funct(*field_it, el);
         }
       }
     } // namespace
   }   // namespace details
 } // namespace fe_engine
 
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct AssembleLumpedTemplateHelper {
       template <template <ElementKind, class> class I,
                 template <ElementKind> class S, ElementKind k, class IOF>
       static void call(const FEEngineTemplate<I, S, k, IOF> & /*unused*/,
                        const std::function<void(Matrix<Real> &,
                                                 const Element &)> & /*unused*/,
                        const ID & /*unused*/, const ID & /*unused*/,
                        DOFManager & /*unused*/, ElementType /*unused*/,
                        GhostType /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define ASSEMBLE_LUMPED(type)                                                  \
   fem.template assembleFieldLumped<type>(field_funct, lumped, dof_id,          \
                                          dof_manager, ghost_type)
 
 #define AKANTU_SPECIALIZE_ASSEMBLE_HELPER(kind)                                \
   template <> struct AssembleLumpedTemplateHelper<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 std::function<void(Matrix<Real> &, const Element &)> &          \
              field_funct,                                                      \
          const ID & lumped, const ID & dof_id, DOFManager & dof_manager,       \
          ElementType type, GhostType ghost_type) {                             \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(ASSEMBLE_LUMPED, kind);                 \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_ASSEMBLE_HELPER,
                                AKANTU_FE_ENGINE_LIST_ASSEMBLE_FIELDS)
 
 #undef AKANTU_SPECIALIZE_ASSEMBLE_HELPER
 #undef AKANTU_SPECIALIZE_ASSEMBLE_HELPER_LIST_KIND
 #undef ASSEMBLE_LUMPED
   } // namespace details
 } // namespace fe_engine
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IOF>
 void FEEngineTemplate<I, S, kind, IOF>::assembleFieldLumped(
     const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
     const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
     ElementType type, GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   fe_engine::details::AssembleLumpedTemplateHelper<kind>::call(
       *this, field_funct, matrix_id, dof_id, dof_manager, type, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::assembleFieldLumped(
     const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
     const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
     GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_degree_of_freedom = dof_manager.getDOFs(dof_id).getNbComponent();
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   UInt nb_integration_points = this->getNbIntegrationPoints(type);
 
   Array<Real> field(0, nb_degree_of_freedom);
   fe_engine::details::fillField(field_funct, field, nb_element,
                                 nb_integration_points, type, ghost_type);
 
   switch (type) {
   case _triangle_6:
   case _quadrangle_8:
   case _tetrahedron_10:
   case _hexahedron_20:
   case _pentahedron_15:
     this->template assembleLumpedDiagonalScaling<type>(field, matrix_id, dof_id,
                                                        dof_manager, ghost_type);
     break;
   default:
     this->template assembleLumpedRowSum<type>(field, matrix_id, dof_id,
                                               dof_manager, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV =
  * \int \rho \varphi_i dV @f$
  */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     assembleLumpedRowSum(const Array<Real> & field, const ID & matrix_id,
                          const ID & dof_id, DOFManager & dof_manager,
                          GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   UInt shapes_size = ElementClass<type>::getShapeSize();
   UInt nb_degree_of_freedom = field.getNbComponent();
 
   auto * field_times_shapes =
       new Array<Real>(0, shapes_size * nb_degree_of_freedom);
 
   shape_functions.template computeNtb<type>(field, *field_times_shapes,
                                             ghost_type);
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   auto * int_field_times_shapes = new Array<Real>(
       nb_element, shapes_size * nb_degree_of_freedom, "inte_rho_x_shapes");
 
   integrator.template integrate<type>(
       *field_times_shapes, *int_field_times_shapes,
       nb_degree_of_freedom * shapes_size, ghost_type, empty_filter);
 
   delete field_times_shapes;
 
   dof_manager.assembleElementalArrayToLumpedMatrix(
       dof_id, *int_field_times_shapes, matrix_id, type, ghost_type);
 
   delete int_field_times_shapes;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
  */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     assembleLumpedDiagonalScaling(const Array<Real> & field,
                                   const ID & matrix_id, const ID & dof_id,
                                   DOFManager & dof_manager,
                                   GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   ElementType type_p1 = ElementClass<type>::getP1ElementType();
   UInt nb_nodes_per_element_p1 = Mesh::getNbNodesPerElement(type_p1);
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_degree_of_freedom = field.getNbComponent();
   UInt nb_element = mesh.getNbElement(type, ghost_type);
 
   Vector<Real> nodal_factor(nb_nodes_per_element);
 
 #define ASSIGN_WEIGHT_TO_NODES(corner, mid)                                    \
   {                                                                            \
     for (UInt n = 0; n < nb_nodes_per_element_p1; n++)                         \
       nodal_factor(n) = corner;                                                \
     for (UInt n = nb_nodes_per_element_p1; n < nb_nodes_per_element; n++)      \
       nodal_factor(n) = mid;                                                   \
   }
 
   if (type == _triangle_6)
     ASSIGN_WEIGHT_TO_NODES(1. / 12., 1. / 4.);
   if (type == _tetrahedron_10)
     ASSIGN_WEIGHT_TO_NODES(1. / 32., 7. / 48.);
   if (type == _quadrangle_8)
     ASSIGN_WEIGHT_TO_NODES(
         3. / 76.,
         16. / 76.); /** coeff. derived by scaling
                      * the diagonal terms of the corresponding
                      * consistent mass computed with 3x3 gauss points;
                      * coeff. are (1./36., 8./36.) for 2x2 gauss points */
   if (type == _hexahedron_20)
     ASSIGN_WEIGHT_TO_NODES(
         7. / 248., 16. / 248.); /** coeff. derived by scaling
                                  * the diagonal terms of the corresponding
                                  * consistent mass computed with 3x3x3 gauss
                                  * points; coeff. are (1./40.,
                                  * 1./15.) for 2x2x2 gauss points */
   if (type == _pentahedron_15) {
     // coefficients derived by scaling the diagonal terms of the corresponding
     // consistent mass computed with 8 gauss points;
     for (UInt n = 0; n < nb_nodes_per_element_p1; n++) {
       nodal_factor(n) = 51. / 2358.;
     }
 
     Real mid_triangle = 192. / 2358.;
     Real mid_quadrangle = 300. / 2358.;
 
     nodal_factor(6) = mid_triangle;
     nodal_factor(7) = mid_triangle;
     nodal_factor(8) = mid_triangle;
     nodal_factor(9) = mid_quadrangle;
     nodal_factor(10) = mid_quadrangle;
     nodal_factor(11) = mid_quadrangle;
     nodal_factor(12) = mid_triangle;
     nodal_factor(13) = mid_triangle;
     nodal_factor(14) = mid_triangle;
   }
 
   if (nb_element == 0) {
     AKANTU_DEBUG_OUT();
     return;
   }
 
 #undef ASSIGN_WEIGHT_TO_NODES
   /// compute @f$ \int \rho dV = \rho V @f$ for each element
   auto int_field = std::make_unique<Array<Real>>(
       field.size(), nb_degree_of_freedom, "inte_rho_x");
   integrator.template integrate<type>(field, *int_field, nb_degree_of_freedom,
                                       ghost_type, empty_filter);
 
   /// distribute the mass of the element to the nodes
   auto lumped_per_node = std::make_unique<Array<Real>>(
       nb_element, nb_degree_of_freedom * nb_nodes_per_element, "mass_per_node");
   auto int_field_it = int_field->begin(nb_degree_of_freedom);
   auto lumped_per_node_it =
       lumped_per_node->begin(nb_degree_of_freedom, nb_nodes_per_element);
 
   for (UInt e = 0; e < nb_element; ++e) {
     for (UInt n = 0; n < nb_nodes_per_element; ++n) {
       Vector<Real> l = (*lumped_per_node_it)(n);
       l = *int_field_it;
       l *= nodal_factor(n);
     }
     ++int_field_it;
     ++lumped_per_node_it;
   }
 
   dof_manager.assembleElementalArrayToLumpedMatrix(dof_id, *lumped_per_node,
                                                    matrix_id, 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 AssembleFieldMatrixHelper {
       template <template <ElementKind, class> class I,
                 template <ElementKind> class S, ElementKind k, class IOF>
       static void call(const FEEngineTemplate<I, S, k, IOF> & /*unused*/,
                        const std::function<void(Matrix<Real> &,
                                                 const Element &)> & /*unused*/,
                        const ID & /*unused*/, const ID & /*unused*/,
                        DOFManager & /*unused*/, ElementType /*unused*/,
                        GhostType /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define ASSEMBLE_MATRIX(type)                                                  \
   fem.template assembleFieldMatrix<type>(field_funct, matrix_id, dof_id,       \
                                          dof_manager, ghost_type)
 
 #define AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER(kind)                   \
   template <> struct AssembleFieldMatrixHelper<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 std::function<void(Matrix<Real> &, const Element &)> &          \
              field_funct,                                                      \
          const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,    \
          ElementType type, GhostType ghost_type) {                             \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(ASSEMBLE_MATRIX, kind);                 \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER,
                                AKANTU_FE_ENGINE_LIST_ASSEMBLE_FIELDS)
 
 #undef AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER
 #undef ASSEMBLE_MATRIX
   } // namespace details
 } // namespace fe_engine
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IOF>
 void FEEngineTemplate<I, S, kind, IOF>::assembleFieldMatrix(
     const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
     const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
     ElementType type, GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   fe_engine::details::AssembleFieldMatrixHelper<kind>::template call(
       *this, field_funct, matrix_id, dof_id, dof_manager, type, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ShapesForMassHelper {
       template <ElementType type, class ShapeFunctions>
       static auto getShapes(ShapeFunctions & shape_functions,
                             const Matrix<Real> & integration_points,
                             const Array<Real> & nodes,
                             UInt & nb_degree_of_freedom, UInt nb_element,
                             GhostType ghost_type) {
 
         UInt shapes_size = ElementClass<type>::getShapeSize();
         Array<Real> shapes(0, shapes_size);
 
         shape_functions.template computeShapesOnIntegrationPoints<type>(
             nodes, integration_points, shapes, ghost_type);
 
         UInt nb_integration_points = integration_points.cols();
         UInt vect_size = nb_integration_points * nb_element;
         UInt lmat_size = nb_degree_of_freedom * shapes_size;
 
         // Extending the shape functions
         /// \todo move this in the shape functions as Voigt format shapes to
         /// have the code in common with the structural elements
         auto shapes_voigt = std::make_unique<Array<Real>>(
             vect_size, lmat_size * nb_degree_of_freedom, 0.);
         auto mshapes_it = shapes_voigt->begin(nb_degree_of_freedom, lmat_size);
         auto shapes_it = shapes.begin(shapes_size);
 
         for (UInt q = 0; q < vect_size; ++q, ++mshapes_it, ++shapes_it) {
           for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
             for (UInt s = 0; s < shapes_size; ++s) {
               (*mshapes_it)(d, s * nb_degree_of_freedom + d) = (*shapes_it)(s);
             }
           }
         }
 
         return shapes_voigt;
       }
     };
 
     template <> struct ShapesForMassHelper<_ek_structural> {
       template <ElementType type, class ShapeFunctions>
       static auto getShapes(ShapeFunctions & shape_functions,
                             const Matrix<Real> & integration_points,
                             const Array<Real> & nodes,
                             UInt & nb_degree_of_freedom, UInt /*nb_element*/,
                             GhostType ghost_type) {
 
-        UInt shapes_size = ElementClass<type>::getShapeMassSize();
-        auto shapes = std::make_unique<Array<Real>>(0, shapes_size);
-        nb_degree_of_freedom = ElementClass<type>::getNaturalSpaceDimension();
+        auto nb_unknown = ElementClass<type>::getNbStressComponents();
+        auto nb_degree_of_freedom_ = ElementClass<type>::getNbDegreeOfFreedom();
+        auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
+        auto shapes = std::make_unique<Array<Real>>(
+            0, nb_unknown * nb_nodes_per_element * nb_degree_of_freedom_);
+        nb_degree_of_freedom = nb_unknown;
         shape_functions.template computeShapesMassOnIntegrationPoints<type>(
             nodes, integration_points, *shapes, ghost_type);
 
         return shapes;
       }
     };
   } // namespace details
 } // namespace fe_engine
   //
 /* -------------------------------------------------------------------------- */
 /**
  * @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV =
  * \int \rho \varphi_i dV @f$
  */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::assembleFieldMatrix(
     const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
     const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
     GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   // \int N * N  so degree 2 * degree of N
   const UInt polynomial_degree =
       2 * ElementClassProperty<type>::polynomial_degree;
 
   // getting the integration points
   Matrix<Real> integration_points =
       integrator.template getIntegrationPoints<type, polynomial_degree>();
 
   UInt nb_degree_of_freedom = dof_manager.getDOFs(dof_id).getNbComponent();
   UInt nb_element = mesh.getNbElement(type, ghost_type);
 
   // getting the shapes on the integration points
   auto shapes_voigt =
       fe_engine::details::ShapesForMassHelper<kind>::template getShapes<type>(
           shape_functions, integration_points, mesh.getNodes(),
           nb_degree_of_freedom, nb_element, ghost_type);
 
   auto vect_size = shapes_voigt->size();
 
   // getting the value to assemble on the integration points
   Array<Real> field(vect_size, nb_degree_of_freedom);
   fe_engine::details::fillField(field_funct, field, nb_element,
                                 integration_points.cols(), type, ghost_type);
 
   auto lmat_size = shapes_voigt->getNbComponent() / nb_degree_of_freedom;
 
   // computing \rho * N
   Array<Real> local_mat(vect_size, lmat_size * lmat_size);
   auto N_it = shapes_voigt->begin(nb_degree_of_freedom, lmat_size);
   auto lmat_it = local_mat.begin(lmat_size, lmat_size);
   auto field_it = field.begin_reinterpret(nb_degree_of_freedom, field.size());
 
   for (UInt q = 0; q < vect_size; ++q, ++lmat_it, ++N_it, ++field_it) {
     const auto & rho = *field_it;
     const auto & N = *N_it;
     auto & mat = *lmat_it;
 
     Matrix<Real> Nt = N.transpose();
     for (UInt d = 0; d < Nt.cols(); ++d) {
       Nt(d) *= rho(d);
     }
 
     mat.template mul<false, false>(Nt, N);
   }
 
   // integrate the elemental values
   Array<Real> int_field_times_shapes(nb_element, lmat_size * lmat_size,
                                      "inte_rho_x_shapes");
   this->integrator.template integrate<type, polynomial_degree>(
       local_mat, int_field_times_shapes, lmat_size * lmat_size, ghost_type);
 
   // assemble the elemental values to the matrix
   dof_manager.assembleElementalMatricesToMatrix(
       matrix_id, dof_id, int_field_times_shapes, type, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
 
 #endif /* AKANTU_FE_ENGINE_TEMPLATE_TMPL_FIELD_HH_ */
diff --git a/src/fe_engine/shape_structural_inline_impl.hh b/src/fe_engine/shape_structural_inline_impl.hh
index fae5b65c7..820898908 100644
--- a/src/fe_engine/shape_structural_inline_impl.hh
+++ b/src/fe_engine/shape_structural_inline_impl.hh
@@ -1,478 +1,509 @@
 /**
  * @file   shape_structural_inline_impl.hh
  *
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Oct 11 2017
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  ShapeStructural inline implementation
  *
  *
  * 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 "mesh_iterators.hh"
 #include "shape_structural.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_
 #define AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_
 
 namespace akantu {
 
 namespace {
   /// Extract nodal coordinates per elements
   template <ElementType type>
   std::unique_ptr<Array<Real>> getNodesPerElement(const Mesh & mesh,
                                                   const Array<Real> & nodes,
                                                   GhostType ghost_type) {
     const auto dim = ElementClass<type>::getSpatialDimension();
     const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
 
     auto nodes_per_element =
         std::make_unique<Array<Real>>(0, dim * nb_nodes_per_element);
     FEEngine::extractNodalToElementField(mesh, nodes, *nodes_per_element, type,
                                          ghost_type);
     return nodes_per_element;
   }
 } // namespace
 
 template <ElementKind kind>
 inline void ShapeStructural<kind>::initShapeFunctions(
     const Array<Real> & /* unused */, const Matrix<Real> & /* unused */,
     ElementType /* unused */, GhostType /* unused */) {
   AKANTU_TO_IMPLEMENT();
 }
 
 /* -------------------------------------------------------------------------- */
 #define INIT_SHAPE_FUNCTIONS(type)                                             \
   setIntegrationPointsByType<type>(integration_points, ghost_type);            \
   precomputeRotationMatrices<type>(nodes, ghost_type);                         \
   precomputeShapesOnIntegrationPoints<type>(nodes, ghost_type);                \
   precomputeShapeDerivativesOnIntegrationPoints<type>(nodes, ghost_type);
 
 template <>
 inline void ShapeStructural<_ek_structural>::initShapeFunctions(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     ElementType type, GhostType ghost_type) {
   AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
 }
 
 #undef INIT_SHAPE_FUNCTIONS
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeStructural<kind>::computeShapesOnIntegrationPointsInternal(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     Array<Real> & shapes, GhostType ghost_type,
     const Array<UInt> & filter_elements, bool mass) const {
 
-  UInt nb_points = integration_points.cols();
-  UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
+  auto nb_points = integration_points.cols();
+  auto nb_element = mesh.getConnectivity(type, ghost_type).size();
+  auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
 
   shapes.resize(nb_element * nb_points);
 
-  UInt ndof = ElementClass<type>::getNbDegreeOfFreedom();
+  auto nb_dofs = ElementClass<type>::getNbDegreeOfFreedom();
+  auto nb_rows = nb_dofs;
+  if (mass) {
+    nb_rows = ElementClass<type>::getNbStressComponents();
+  }
 
 #if !defined(AKANTU_NDEBUG)
-  UInt size_of_shapes = ElementClass<type>::getShapeSize();
+  UInt size_of_shapes = nb_rows * nb_dofs * nb_nodes_per_element;
   AKANTU_DEBUG_ASSERT(shapes.getNbComponent() == size_of_shapes,
                       "The shapes array does not have the correct "
                           << "number of component");
 #endif
 
-  auto nb_rows = ndof;
-  if (mass) {
-    nb_rows = ElementClass<type>::getNaturalSpaceDimension();
-  }
 
   auto shapes_it = shapes.begin_reinterpret(
-      nb_rows, ElementClass<type>::getNbNodesPerInterpolationElement() * ndof,
+      nb_rows, ElementClass<type>::getNbNodesPerInterpolationElement() * nb_dofs,
       nb_points, nb_element);
-
   auto shapes_begin = shapes_it;
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 
   auto nodes_per_element = getNodesPerElement<type>(mesh, nodes, ghost_type);
   auto nodes_it = nodes_per_element->begin(mesh.getSpatialDimension(),
                                            Mesh::getNbNodesPerElement(type));
   auto nodes_begin = nodes_it;
+  auto rot_matrix_it =
+      make_view(rotation_matrices(type, ghost_type), nb_dofs, nb_dofs).begin();
+  auto rot_matrix_begin = rot_matrix_it;
 
   for (UInt elem = 0; elem < nb_element; ++elem) {
     if (filter_elements != empty_filter) {
       shapes_it = shapes_begin + filter_elements(elem);
       nodes_it = nodes_begin + filter_elements(elem);
+      rot_matrix_it = rot_matrix_begin + filter_elements(elem);
     }
 
     Tensor3<Real> & N = *shapes_it;
     auto & real_coord = *nodes_it;
 
-    if(not mass) {
-      ElementClass<type>::computeShapes(integration_points, real_coord, N);
+    auto & RDOFs = *rot_matrix_it;
+
+    Matrix<Real> T(N.size(1), N.size(1), 0);
+
+    for (UInt i = 0; i < nb_nodes_per_element; ++i) {
+      T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
+    }
+
+    if (not mass) {
+      ElementClass<type>::computeShapes(integration_points, real_coord, T, N);
     } else {
-      ElementClass<type>::computeShapesMass(integration_points, real_coord, N);
+      ElementClass<type>::computeShapesMass(integration_points, real_coord, T,
+                                            N);
     }
     if (filter_elements == empty_filter) {
       ++shapes_it;
       ++nodes_it;
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeStructural<kind>::precomputeRotationMatrices(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   const auto spatial_dimension = mesh.getSpatialDimension();
   const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   const auto nb_element = mesh.getNbElement(type, ghost_type);
   const auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
 
   if (not this->rotation_matrices.exists(type, ghost_type)) {
     this->rotation_matrices.alloc(0, nb_dof * nb_dof, type, ghost_type);
   }
 
   auto & rot_matrices = this->rotation_matrices(type, ghost_type);
   rot_matrices.resize(nb_element);
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type);
 
   bool has_extra_normal = mesh.hasData<Real>("extra_normal", type, ghost_type);
   Array<Real>::const_vector_iterator extra_normal;
   if (has_extra_normal) {
     extra_normal = mesh.getData<Real>("extra_normal", type, ghost_type)
                        .begin(spatial_dimension);
   }
 
   for (auto && tuple :
        zip(make_view(x_el, spatial_dimension, nb_nodes_per_element),
            make_view(rot_matrices, nb_dof, nb_dof))) {
     // compute shape derivatives
     auto & X = std::get<0>(tuple);
     auto & R = std::get<1>(tuple);
 
     if (has_extra_normal) {
       ElementClass<type>::computeRotationMatrix(R, X, *extra_normal);
       ++extra_normal;
     } else {
       ElementClass<type>::computeRotationMatrix(
           R, X, Vector<Real>(spatial_dimension));
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <ElementKind kind>
 template <ElementType type>
 void ShapeStructural<kind>::precomputeShapesOnIntegrationPoints(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   const auto & natural_coords = integration_points(type, ghost_type);
   auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   auto nb_points = integration_points(type, ghost_type).cols();
   auto nb_element = mesh.getNbElement(type, ghost_type);
   auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
   const auto dim = ElementClass<type>::getSpatialDimension();
+  const auto spatial_dimension = mesh.getSpatialDimension();
+  const auto natural_spatial_dimension =
+      ElementClass<type>::getNaturalSpaceDimension();
 
   auto itp_type = FEEngine::getInterpolationType(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);
   }
 
+  auto & rot_matrices = this->rotation_matrices(type, ghost_type);
   auto & shapes_ = this->shapes(itp_type, ghost_type);
   shapes_.resize(nb_element * nb_points);
 
   auto nodes_per_element = getNodesPerElement<type>(mesh, nodes, ghost_type);
 
   for (auto && tuple :
        zip(make_view(shapes_, nb_dof, nb_dof * nb_nodes_per_element, nb_points),
-           make_view(*nodes_per_element, dim, nb_nodes_per_element))) {
+           make_view(*nodes_per_element, dim, nb_nodes_per_element),
+           make_view(rot_matrices, nb_dof, nb_dof))) {
     auto & N = std::get<0>(tuple);
-    auto & real_coord = std::get<1>(tuple);
-    ElementClass<type>::computeShapes(natural_coords, real_coord, N);
+    auto & X = std::get<1>(tuple);
+    auto & RDOFs = std::get<2>(tuple);
+
+    Matrix<Real> T(N.size(1), N.size(1), 0);
+
+    for (UInt i = 0; i < nb_nodes_per_element; ++i) {
+      T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
+    }
+
+    auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension);
+    // Rotate to local basis
+    auto x =
+        (R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element);
+
+    ElementClass<type>::computeShapes(natural_coords, x, T, N);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeStructural<kind>::precomputeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   const auto & natural_coords = integration_points(type, ghost_type);
   const auto spatial_dimension = mesh.getSpatialDimension();
   const auto natural_spatial_dimension =
       ElementClass<type>::getNaturalSpaceDimension();
   const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   const auto nb_points = natural_coords.cols();
   const auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
   const auto nb_element = mesh.getNbElement(type, ghost_type);
   const auto nb_stress_components = ElementClass<type>::getNbStressComponents();
 
   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);
   }
 
   auto & rot_matrices = this->rotation_matrices(type, ghost_type);
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type);
 
   auto & shapesd = this->shapes_derivatives(itp_type, ghost_type);
   shapesd.resize(nb_element * nb_points);
 
   for (auto && tuple :
        zip(make_view(x_el, spatial_dimension, nb_nodes_per_element),
            make_view(shapesd, nb_stress_components,
                      nb_nodes_per_element * nb_dof, nb_points),
            make_view(rot_matrices, nb_dof, nb_dof))) {
     // compute shape derivatives
     auto & X = std::get<0>(tuple);
     auto & B = std::get<1>(tuple);
     auto & RDOFs = std::get<2>(tuple);
 
     Tensor3<Real> dnds(natural_spatial_dimension,
                        ElementClass<type>::interpolation_property::dnds_columns,
                        B.size(2));
     ElementClass<type>::computeDNDS(natural_coords, X, dnds);
 
     Tensor3<Real> J(natural_spatial_dimension, natural_spatial_dimension,
                     natural_coords.cols());
 
     // Computing the coordinates of the element in the natural space
     auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension);
     Matrix<Real> T(B.size(1), B.size(1), 0);
 
     for (UInt i = 0; i < nb_nodes_per_element; ++i) {
       T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
     }
 
     // Rotate to local basis
     auto x =
         (R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element);
 
     ElementClass<type>::computeJMat(natural_coords, x, J);
     ElementClass<type>::computeShapeDerivatives(J, dnds, T, B);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeStructural<kind>::interpolateOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_dof,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(out_uq.getNbComponent() == nb_dof,
                       "The output array shape is not correct");
 
   auto itp_type = FEEngine::getInterpolationType(type);
   const auto & shapes_ = shapes(itp_type, ghost_type);
 
   auto nb_element = mesh.getNbElement(type, ghost_type);
   auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
   auto nb_quad_points_per_element = integration_points(type, ghost_type).cols();
 
   Array<Real> u_el(0, nb_nodes_per_element * nb_dof);
   FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
                                        filter_elements);
 
   auto nb_quad_points = nb_quad_points_per_element * u_el.size();
   out_uq.resize(nb_quad_points);
 
   auto out_it = out_uq.begin_reinterpret(nb_dof, 1, nb_quad_points_per_element,
                                          u_el.size());
   auto shapes_it =
       shapes_.begin_reinterpret(nb_dof, nb_dof * nb_nodes_per_element,
                                 nb_quad_points_per_element, nb_element);
   auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1,
                                      nb_quad_points_per_element, u_el.size());
 
   for_each_element(nb_element, filter_elements, [&](auto && el) {
     auto & uq = *out_it;
     const auto & u = *u_it;
     auto N = Tensor3<Real>(shapes_it[el]);
 
     for (auto && q : arange(uq.size(2))) {
       auto uq_q = Matrix<Real>(uq(q));
       auto u_q = Matrix<Real>(u(q));
       auto N_q = Matrix<Real>(N(q));
 
       uq_q.mul<false, false>(N_q, u_q);
     }
 
     ++out_it;
     ++u_it;
   });
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeStructural<kind>::gradientOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_nablauq, UInt nb_dof,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   auto itp_type = FEEngine::getInterpolationType(type);
   const auto & shapesd = shapes_derivatives(itp_type, ghost_type);
 
   auto nb_element = mesh.getNbElement(type, ghost_type);
   auto element_dimension = ElementClass<type>::getSpatialDimension();
   auto nb_quad_points_per_element = integration_points(type, ghost_type).cols();
   auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
 
   Array<Real> u_el(0, nb_nodes_per_element * nb_dof);
   FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
                                        filter_elements);
 
   auto nb_quad_points = nb_quad_points_per_element * u_el.size();
   out_nablauq.resize(nb_quad_points);
 
   auto out_it = out_nablauq.begin_reinterpret(
       element_dimension, 1, nb_quad_points_per_element, u_el.size());
   auto shapesd_it = shapesd.begin_reinterpret(
       element_dimension, nb_dof * nb_nodes_per_element,
       nb_quad_points_per_element, nb_element);
   auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1,
                                      nb_quad_points_per_element, u_el.size());
 
   for_each_element(nb_element, filter_elements, [&](auto && el) {
     auto & nablau = *out_it;
     const auto & u = *u_it;
     auto B = Tensor3<Real>(shapesd_it[el]);
 
     for (auto && q : arange(nablau.size(2))) {
       auto nablau_q = Matrix<Real>(nablau(q));
       auto u_q = Matrix<Real>(u(q));
       auto B_q = Matrix<Real>(B(q));
 
       nablau_q.mul<false, false>(B_q, u_q);
     }
 
     ++out_it;
     ++u_it;
   });
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 template <ElementType type>
 void ShapeStructural<_ek_structural>::computeBtD(
     const Array<Real> & Ds, Array<Real> & BtDs, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   auto itp_type = ElementClassProperty<type>::interpolation_type;
 
   auto nb_stress = ElementClass<type>::getNbStressComponents();
   auto nb_dof_per_element = ElementClass<type>::getNbDegreeOfFreedom() *
                             mesh.getNbNodesPerElement(type);
 
   const auto & shapes_derivatives =
       this->shapes_derivatives(itp_type, ghost_type);
 
   Array<Real> shapes_derivatives_filtered(0,
                                           shapes_derivatives.getNbComponent());
   auto && view = make_view(shapes_derivatives, nb_stress, nb_dof_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, nb_stress, nb_dof_per_element);
     B_it = view.begin();
     B_end = view.end();
   }
 
   for (auto && values : zip(range(B_it, B_end), make_view(Ds, nb_stress),
                             make_view(BtDs, BtDs.getNbComponent()))) {
     const auto & B = std::get<0>(values);
     const auto & D = std::get<1>(values);
     auto & Bt_D = std::get<2>(values);
     Bt_D.template mul<true>(B, D);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 template <ElementType type>
 void ShapeStructural<_ek_structural>::computeNtb(
     const Array<Real> & bs, Array<Real> & Ntbs, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   auto itp_type = ElementClassProperty<type>::interpolation_type;
 
   auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
   auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
 
   const auto & shapes = this->shapes(itp_type, ghost_type);
 
   Array<Real> shapes_filtered(0, shapes.getNbComponent());
   auto && view = make_view(shapes, nb_dof, nb_dof * nb_nodes_per_element);
   auto N_it = view.begin();
   auto N_end = view.end();
 
   if (filter_elements != empty_filter) {
     FEEngine::filterElementalData(this->mesh, shapes, shapes_filtered, type,
                                   ghost_type, filter_elements);
     auto && view =
         make_view(shapes_filtered, nb_dof, nb_dof * nb_nodes_per_element);
     N_it = view.begin();
     N_end = view.end();
   }
 
   for (auto && values : zip(range(N_it, N_end), make_view(bs, nb_dof),
                             make_view(Ntbs, nb_dof * nb_nodes_per_element))) {
     const auto & N = std::get<0>(values);
     const auto & b = std::get<1>(values);
     auto & Nt_b = std::get<2>(values);
     Nt_b.template mul<true>(N, b);
   }
 }
 
 } // namespace akantu
 
 #endif /* AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_ */
diff --git a/src/model/model.cc b/src/model/model.cc
index 8071891b8..ac974fd7e 100644
--- a/src/model/model.cc
+++ b/src/model/model.cc
@@ -1,339 +1,341 @@
 /**
  * @file   model.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Mon Oct 03 2011
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  implementation of model common parts
  *
  *
  * 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 "model.hh"
 #include "communicator.hh"
 #include "data_accessor.hh"
 #include "element_group.hh"
 #include "element_synchronizer.hh"
 #include "synchronizer_registry.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 Model::Model(Mesh & mesh, const ModelType & type, UInt dim, const ID & id,
              const MemoryID & memory_id)
     : Memory(id, memory_id), ModelSolver(mesh, type, id, memory_id), mesh(mesh),
       spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension()
                                                : dim),
       parser(getStaticParser()) {
   AKANTU_DEBUG_IN();
 
   this->mesh.registerEventHandler(*this, _ehp_model);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 Model::~Model() = default;
 
 /* -------------------------------------------------------------------------- */
 void Model::initFullImpl(const ModelOptions & options) {
   AKANTU_DEBUG_IN();
 
   method = options.analysis_method;
   if (!this->hasDefaultSolver()) {
     this->initNewSolver(this->method);
   }
   initModel();
 
   initFEEngineBoundary();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::initNewSolver(const AnalysisMethod & method) {
   ID solver_name;
   TimeStepSolverType tss_type;
   std::tie(solver_name, tss_type) = this->getDefaultSolverID(method);
 
   if (not this->hasSolver(solver_name)) {
     ModelSolverOptions options = this->getDefaultSolverOptions(tss_type);
     this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type);
 
     for (auto && is_type : options.integration_scheme_type) {
       if (!this->hasIntegrationScheme(solver_name, is_type.first)) {
         this->setIntegrationScheme(solver_name, is_type.first, is_type.second,
                                    options.solution_type[is_type.first]);
       }
     }
   }
 
   this->method = method;
   this->setDefaultSolver(solver_name);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::initFEEngineBoundary() {
   FEEngine & fem_boundary = getFEEngineBoundary();
   fem_boundary.initShapeFunctions(_not_ghost);
   fem_boundary.initShapeFunctions(_ghost);
 
   fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost);
   fem_boundary.computeNormalsOnIntegrationPoints(_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::dumpGroup(const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.dump();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::dumpGroup(const std::string & group_name,
                       const std::string & dumper_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.dump(dumper_name);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::dumpGroup() {
   for (auto & group : mesh.iterateElementGroups()) {
     group.dump();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::setGroupDirectory(const std::string & directory) {
   for (auto & group : mesh.iterateElementGroups()) {
     group.setDirectory(directory);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::setGroupDirectory(const std::string & directory,
                               const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.setDirectory(directory);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::setGroupBaseName(const std::string & basename,
                              const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.setBaseName(basename);
 }
 
 /* -------------------------------------------------------------------------- */
 DumperIOHelper & Model::getGroupDumper(const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   return group.getDumper();
 }
 
 /* -------------------------------------------------------------------------- */
 // DUMPER stuff
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldToDumper(const std::string & field_id,
                                       std::shared_ptr<dumpers::Field> field,
                                       DumperIOHelper & dumper) {
 #ifdef AKANTU_USE_IOHELPER
   dumper.registerField(field_id, std::move(field));
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpField(const std::string & field_id) {
 
   this->addDumpFieldToDumper(mesh.getDefaultDumperName(), field_id);
 }
 /* -------------------------------------------------------------------------- */
 
 void Model::addDumpFieldVector(const std::string & field_id) {
 
   this->addDumpFieldVectorToDumper(mesh.getDefaultDumperName(), field_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpFieldTensor(const std::string & field_id) {
 
   this->addDumpFieldTensorToDumper(mesh.getDefaultDumperName(), field_id);
 }
 
 /* -------------------------------------------------------------------------- */
 
 void Model::setBaseName(const std::string & field_id) {
 
   mesh.setBaseName(field_id);
 }
 /* -------------------------------------------------------------------------- */
 
 void Model::setBaseNameToDumper(const std::string & dumper_name,
                                 const std::string & basename) {
   mesh.setBaseNameToDumper(dumper_name, basename);
 }
 /* -------------------------------------------------------------------------- */
 
 void Model::addDumpFieldToDumper(const std::string & dumper_name,
                                  const std::string & field_id) {
 
   this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
                                   false);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupField(const std::string & field_id,
                               const std::string & group_name) {
 
   ElementGroup & group = mesh.getElementGroup(group_name);
   this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id,
                                   group_name, _ek_regular, false);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::removeDumpGroupField(const std::string & field_id,
                                  const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id,
                                        group_name);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::removeDumpGroupFieldFromDumper(const std::string & dumper_name,
                                            const std::string & field_id,
                                            const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.removeDumpFieldFromDumper(dumper_name, field_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpFieldVectorToDumper(const std::string & dumper_name,
                                        const std::string & field_id) {
   this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
                                   true);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldVector(const std::string & field_id,
                                     const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id,
                                         group_name);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldVectorToDumper(const std::string & dumper_name,
                                             const std::string & field_id,
                                             const std::string & group_name) {
   this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name,
                                   _ek_regular, true);
 }
 /* -------------------------------------------------------------------------- */
 
 void Model::addDumpFieldTensorToDumper(const std::string & dumper_name,
                                        const std::string & field_id) {
   this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
                                   true);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldToDumper(const std::string & dumper_name,
                                       const std::string & field_id,
                                       const std::string & group_name,
                                       ElementKind element_kind,
                                       bool padding_flag) {
   this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name,
                                   this->spatial_dimension, element_kind,
                                   padding_flag);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldToDumper(const std::string & dumper_name,
                                       const std::string & field_id,
                                       const std::string & group_name,
                                       UInt spatial_dimension,
                                       ElementKind element_kind,
                                       bool padding_flag) {
 
 #ifdef AKANTU_USE_IOHELPER
   std::shared_ptr<dumpers::Field> field;
 
   if (!field) {
     field = this->createNodalFieldReal(field_id, group_name, padding_flag);
   }
   if (!field) {
     field = this->createNodalFieldUInt(field_id, group_name, padding_flag);
   }
   if (!field) {
     field = this->createNodalFieldBool(field_id, group_name, padding_flag);
   }
   if (!field) {
     field = this->createElementalField(field_id, group_name, padding_flag,
                                        spatial_dimension, element_kind);
   }
   if (!field) {
     field = this->mesh.createFieldFromAttachedData<UInt>(field_id, group_name,
                                                          element_kind);
   }
   if (!field) {
     field = this->mesh.createFieldFromAttachedData<Real>(field_id, group_name,
                                                          element_kind);
   }
 
 #ifndef AKANTU_NDEBUG
   if (!field) {
     AKANTU_DEBUG_WARNING("No field could be found based on name: " << field_id);
   }
 #endif
   if (field) {
     DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name, group_name);
     this->addDumpGroupFieldToDumper(field_id, field, dumper);
   }
 
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
-
 void Model::dump() { mesh.dump(); }
 
 /* -------------------------------------------------------------------------- */
+void Model::dump(UInt step) { mesh.dump(step); }
+
+/* -------------------------------------------------------------------------- */
+void Model::dump(Real time, UInt step) { mesh.dump(time, step); }
 
+/* -------------------------------------------------------------------------- */
 void Model::setDirectory(const std::string & directory) {
   mesh.setDirectory(directory);
 }
 
 /* -------------------------------------------------------------------------- */
-
 void Model::setDirectoryToDumper(const std::string & dumper_name,
                                  const std::string & directory) {
   mesh.setDirectoryToDumper(dumper_name, directory);
 }
 
 /* -------------------------------------------------------------------------- */
-
 void Model::setTextModeToDumper() { mesh.setTextModeToDumper(); }
 
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu
diff --git a/src/model/model.hh b/src/model/model.hh
index 08d099325..f70087c26 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,349 +1,351 @@
 /**
  * @file   model.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Interface of a model
  *
  *
  * 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 "aka_memory.hh"
 #include "aka_named_argument.hh"
 #include "fe_engine.hh"
 #include "mesh.hh"
 #include "model_options.hh"
 #include "model_solver.hh"
 /* -------------------------------------------------------------------------- */
 #include <typeindex>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_MODEL_HH_
 #define AKANTU_MODEL_HH_
 
 namespace akantu {
 class SynchronizerRegistry;
 class Parser;
 class DumperIOHelper;
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 
 class Model : public Memory, public ModelSolver, public MeshEventHandler {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   Model(Mesh & mesh, const ModelType & type, UInt dim = _all_dimensions,
         const ID & id = "model", const MemoryID & memory_id = 0);
 
   ~Model() override;
 
   using FEEngineMap = std::map<std::string, std::unique_ptr<FEEngine>>;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 protected:
   virtual void initFullImpl(const ModelOptions & options);
 
 public:
   template <typename... pack>
   std::enable_if_t<are_named_argument<pack...>::value>
   initFull(pack &&... _pack) {
     switch (this->model_type) {
 #ifdef AKANTU_SOLID_MECHANICS
     case ModelType::_solid_mechanics_model:
       this->initFullImpl(SolidMechanicsModelOptions{
           use_named_args, std::forward<decltype(_pack)>(_pack)...});
       break;
 #endif
 #ifdef AKANTU_COHESIVE_ELEMENT
     case ModelType::_solid_mechanics_model_cohesive:
       this->initFullImpl(SolidMechanicsModelCohesiveOptions{
           use_named_args, std::forward<decltype(_pack)>(_pack)...});
       break;
 #endif
 #ifdef AKANTU_HEAT_TRANSFER
     case ModelType::_heat_transfer_model:
       this->initFullImpl(HeatTransferModelOptions{
           use_named_args, std::forward<decltype(_pack)>(_pack)...});
       break;
 #endif
 #ifdef AKANTU_EMBEDDED
     case ModelType::_embedded_model:
       this->initFullImpl(EmbeddedInterfaceModelOptions{
           use_named_args, std::forward<decltype(_pack)>(_pack)...});
       break;
 #endif
     default:
       this->initFullImpl(ModelOptions{use_named_args,
                                       std::forward<decltype(_pack)>(_pack)...});
     }
   }
 
   template <typename... pack>
   std::enable_if_t<not are_named_argument<pack...>::value>
   initFull(pack &&... _pack) {
     this->initFullImpl(std::forward<decltype(_pack)>(_pack)...);
   }
 
   /// initialize a new solver if needed
   void initNewSolver(const AnalysisMethod & method);
 
 protected:
   /// get some default values for derived classes
   virtual std::tuple<ID, TimeStepSolverType>
   getDefaultSolverID(const AnalysisMethod & method) = 0;
 
   virtual void initModel() = 0;
 
   virtual void initFEEngineBoundary();
 
   /// function to print the containt of the class
   void printself(std::ostream & /*stream*/,
                  int /*indent*/ = 0) const override{};
 
 public:
   /* ------------------------------------------------------------------------ */
   /* Access to the dumpable interface of the boundaries                       */
   /* ------------------------------------------------------------------------ */
   /// Dump the data for a given group
   void dumpGroup(const std::string & group_name);
   void dumpGroup(const std::string & group_name,
                  const std::string & dumper_name);
   /// Dump the data for all boundaries
   void dumpGroup();
   /// Set the directory for a given group
   void setGroupDirectory(const std::string & directory,
                          const std::string & group_name);
   /// Set the directory for all boundaries
   void setGroupDirectory(const std::string & directory);
   /// Set the base name for a given group
   void setGroupBaseName(const std::string & basename,
                         const std::string & group_name);
   /// Get the internal dumper of a given group
   DumperIOHelper & getGroupDumper(const std::string & group_name);
 
   /* ------------------------------------------------------------------------ */
   /* Function for non local capabilities                                      */
   /* ------------------------------------------------------------------------ */
   virtual void updateDataForNonLocalCriterion(__attribute__((unused))
                                               ElementTypeMapReal & criterion) {
     AKANTU_TO_IMPLEMENT();
   }
 
 protected:
   template <typename T>
   void allocNodalField(Array<T> *& array, UInt nb_component,
                        const ID & name);
 
   template <typename T>
   void allocNodalField(std::unique_ptr<Array<T>> & array, UInt nb_component,
                        const ID & name) const;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors */
   /* ------------------------------------------------------------------------ */
 public:
   /// get id of model
   AKANTU_GET_MACRO(ID, id, const ID &)
 
   /// get the number of surfaces
   AKANTU_GET_MACRO(Mesh, mesh, Mesh &)
 
   /// synchronize the boundary in case of parallel run
   virtual void synchronizeBoundaries(){};
 
   /// return the fem object associated with a provided name
   inline FEEngine & getFEEngine(const ID & name = "") const;
 
   /// return the fem boundary object associated with a provided name
   virtual FEEngine & getFEEngineBoundary(const ID & name = "");
 
   /// register a fem object associated with name
   template <typename FEEngineClass>
   inline void registerFEEngineObject(const std::string & name, Mesh & mesh,
                                      UInt spatial_dimension);
   /// unregister a fem object associated with name
   inline void unRegisterFEEngineObject(const std::string & name);
 
   /// return the synchronizer registry
   SynchronizerRegistry & getSynchronizerRegistry();
 
   /// return the fem object associated with a provided name
   template <typename FEEngineClass>
   inline FEEngineClass & getFEEngineClass(std::string name = "") const;
 
   /// return the fem boundary object associated with a provided name
   template <typename FEEngineClass>
   inline FEEngineClass & getFEEngineClassBoundary(std::string name = "");
 
   /// Get the type of analysis method used
   AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod);
 
   /* ------------------------------------------------------------------------ */
   /* Pack and unpack helper functions                                         */
   /* ------------------------------------------------------------------------ */
 public:
   inline UInt getNbIntegrationPoints(const Array<Element> & elements,
                                      const ID & fem_id = ID()) const;
 
   /* ------------------------------------------------------------------------ */
   /* Dumpable interface (kept for convenience) and dumper relative functions  */
   /* ------------------------------------------------------------------------ */
   void setTextModeToDumper();
 
   virtual void addDumpGroupFieldToDumper(const std::string & field_id,
                                          std::shared_ptr<dumpers::Field> field,
                                          DumperIOHelper & dumper);
 
   virtual void addDumpField(const std::string & field_id);
 
   virtual void addDumpFieldVector(const std::string & field_id);
 
   virtual void addDumpFieldToDumper(const std::string & dumper_name,
                                     const std::string & field_id);
 
   virtual void addDumpFieldVectorToDumper(const std::string & dumper_name,
                                           const std::string & field_id);
 
   virtual void addDumpFieldTensorToDumper(const std::string & dumper_name,
                                           const std::string & field_id);
 
   virtual void addDumpFieldTensor(const std::string & field_id);
 
   virtual void setBaseName(const std::string & field_id);
 
   virtual void setBaseNameToDumper(const std::string & dumper_name,
                                    const std::string & basename);
 
   virtual void addDumpGroupField(const std::string & field_id,
                                  const std::string & group_name);
 
   virtual void addDumpGroupFieldToDumper(const std::string & dumper_name,
                                          const std::string & field_id,
                                          const std::string & group_name,
                                          ElementKind element_kind,
                                          bool padding_flag);
 
   virtual void addDumpGroupFieldToDumper(const std::string & dumper_name,
                                          const std::string & field_id,
                                          const std::string & group_name,
                                          UInt spatial_dimension,
                                          ElementKind element_kind,
                                          bool padding_flag);
 
   virtual void removeDumpGroupField(const std::string & field_id,
                                     const std::string & group_name);
   virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name,
                                               const std::string & field_id,
                                               const std::string & group_name);
 
   virtual void addDumpGroupFieldVector(const std::string & field_id,
                                        const std::string & group_name);
 
   virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name,
                                                const std::string & field_id,
                                                const std::string & group_name);
 
   virtual std::shared_ptr<dumpers::Field>
   createNodalFieldReal(const std::string & /*field_name*/,
                        const std::string & /*group_name*/,
                        bool /*padding_flag*/) {
     return nullptr;
   }
 
   virtual std::shared_ptr<dumpers::Field>
   createNodalFieldUInt(const std::string & /*field_name*/,
                        const std::string & /*group_name*/,
                        bool /*padding_flag*/) {
     return nullptr;
   }
 
   virtual std::shared_ptr<dumpers::Field>
   createNodalFieldBool(const std::string & /*field_name*/,
                        const std::string & /*group_name*/,
                        bool /*padding_flag*/) {
     return nullptr;
   }
 
   virtual std::shared_ptr<dumpers::Field>
   createElementalField(const std::string & /*field_name*/,
                        const std::string & /*group_name*/,
                        bool /*padding_flag*/,
                        UInt /*spatial_dimension*/,
                        ElementKind /*kind*/) {
     return nullptr;
   }
 
   void setDirectory(const std::string & directory);
   void setDirectoryToDumper(const std::string & dumper_name,
                             const std::string & directory);
 
   virtual void dump();
+  virtual void dump(UInt step);
+  virtual void dump(Real time, UInt step);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   friend std::ostream & operator<<(std::ostream & /*stream*/,
                                    const Model & /*_this*/);
 
   /// analysis method check the list in akantu::AnalysisMethod
   AnalysisMethod method;
 
   /// Mesh
   Mesh & mesh;
 
   /// Spatial dimension of the problem
   UInt spatial_dimension;
 
   /// the main fem object present in all  models
   FEEngineMap fems;
 
   /// the fem object present in all  models for boundaries
   FEEngineMap fems_boundary;
 
   /// default fem object
   std::string default_fem;
 
   /// parser to the pointer to use
   Parser & parser;
 };
 
 /// standard output stream operator
 inline std::ostream & operator<<(std::ostream & stream, const Model & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "model_inline_impl.hh"
 
 #endif /* AKANTU_MODEL_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh
index a5d896289..b59c314dd 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model.hh
@@ -1,581 +1,577 @@
 /**
  * @file   solid_mechanics_model.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Jul 27 2010
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  Model of Solid Mechanics
  *
  *
  * 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 "boundary_condition.hh"
 #include "data_accessor.hh"
 #include "fe_engine.hh"
 #include "model.hh"
 #include "non_local_manager_callback.hh"
 #include "solid_mechanics_model_event_handler.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_
 #define AKANTU_SOLID_MECHANICS_MODEL_HH_
 
 namespace akantu {
 class Material;
 class MaterialSelector;
 class DumperIOHelper;
 class NonLocalManager;
 template <ElementKind kind, class IntegrationOrderFunctor>
 class IntegratorGauss;
 template <ElementKind kind> class ShapeLagrange;
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 class SolidMechanicsModel
     : public Model,
       public DataAccessor<Element>,
       public DataAccessor<UInt>,
       public BoundaryCondition<SolidMechanicsModel>,
       public NonLocalManagerCallback,
       public EventHandlerManager<SolidMechanicsModelEventHandler> {
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   class NewMaterialElementsEvent : public NewElementsEvent {
   public:
     AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array<UInt> &);
     AKANTU_GET_MACRO(MaterialList, material, const Array<UInt> &);
 
   protected:
     Array<UInt> material;
   };
 
   using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
 
 protected:
   using EventManager = EventHandlerManager<SolidMechanicsModelEventHandler>;
 
 public:
   SolidMechanicsModel(
       Mesh & mesh, UInt dim = _all_dimensions,
       const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0,
       ModelType model_type = ModelType::_solid_mechanics_model);
 
   ~SolidMechanicsModel() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 protected:
   /// initialize completely the model
   void initFullImpl(
       const ModelOptions & options = SolidMechanicsModelOptions()) override;
 
 public:
   /// initialize all internal arrays for materials
   virtual void initMaterials();
 
 protected:
   /// initialize the model
   void initModel() override;
 
   /// function to print the containt of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   /// get some default values for derived classes
   std::tuple<ID, TimeStepSolverType>
   getDefaultSolverID(const AnalysisMethod & method) override;
 
   /* ------------------------------------------------------------------------ */
   /* Solver interface                                                         */
   /* ------------------------------------------------------------------------ */
 public:
   /// assembles the stiffness matrix,
   virtual void assembleStiffnessMatrix();
   /// assembles the internal forces in the array internal_forces
   virtual void assembleInternalForces();
 
 protected:
   /// callback for the solver, this adds f_{ext} - f_{int} to the residual
   void assembleResidual() override;
 
   /// callback for the solver, this adds f_{ext} or  f_{int} to the residual
   void assembleResidual(const ID & residual_part) override;
   bool canSplitResidual() override { return true; }
 
   /// get the type of matrix needed
   MatrixType getMatrixType(const ID & matrix_id) override;
 
   /// callback for the solver, this assembles different matrices
   void assembleMatrix(const ID & matrix_id) override;
 
   /// callback for the solver, this assembles the stiffness matrix
   void assembleLumpedMatrix(const ID & matrix_id) override;
 
   /// callback for the solver, this is called at beginning of solve
   void predictor() override;
   /// callback for the solver, this is called at end of solve
   void corrector() override;
 
   /// callback for the solver, this is called at beginning of solve
   void beforeSolveStep() override;
   /// callback for the solver, this is called at end of solve
   void afterSolveStep(bool converged = true) override;
 
   /// Callback for the model to instantiate the matricees when needed
   void initSolver(TimeStepSolverType time_step_solver_type,
                   NonLinearSolverType non_linear_solver_type) override;
 
 protected:
   /* ------------------------------------------------------------------------ */
   TimeStepSolverType getDefaultSolverType() const override;
   /* ------------------------------------------------------------------------ */
   ModelSolverOptions
   getDefaultSolverOptions(const TimeStepSolverType & type) const override;
 
 public:
   bool isDefaultSolverExplicit() {
     return method == _explicit_lumped_mass ||
            method == _explicit_consistent_mass;
   }
 
 protected:
   /// update the current position vector
   void updateCurrentPosition();
 
   /* ------------------------------------------------------------------------ */
   /* Materials (solid_mechanics_model_material.cc)                            */
   /* ------------------------------------------------------------------------ */
 public:
   /// register an empty material of a given type
   Material & registerNewMaterial(const ID & mat_name, const ID & mat_type,
                                  const ID & opt_param);
 
   /// reassigns materials depending on the material selector
   virtual void reassignMaterial();
 
   /// apply a constant eigen_grad_u on all quadrature points of a given material
   virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
                                const ID & material_name,
                                GhostType ghost_type = _not_ghost);
 
 protected:
   /// register a material in the dynamic database
   Material & registerNewMaterial(const ParserSection & mat_section);
 
   /// read the material files to instantiate all the materials
   void instantiateMaterials();
 
   /// set the element_id_by_material and add the elements to the good materials
   virtual void
   assignMaterialToElements(const ElementTypeMapArray<UInt> * filter = nullptr);
 
   /* ------------------------------------------------------------------------ */
   /* Mass (solid_mechanics_model_mass.cc)                                     */
   /* ------------------------------------------------------------------------ */
 public:
   /// assemble the lumped mass matrix
   void assembleMassLumped();
 
   /// assemble the mass matrix for consistent mass resolutions
   void assembleMass();
 
 protected:
   /// assemble the lumped mass matrix for local and ghost elements
   void assembleMassLumped(GhostType ghost_type);
 
   /// assemble the mass matrix for either _ghost or _not_ghost elements
   void assembleMass(GhostType ghost_type);
 
   /// fill a vector of rho
   void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type);
 
   /// compute the kinetic energy
   Real getKineticEnergy();
   Real getKineticEnergy(ElementType type, UInt index);
 
   /// compute the external work (for impose displacement, the velocity should be
   /// given too)
   Real getExternalWork();
 
   /* ------------------------------------------------------------------------ */
   /* NonLocalManager inherited members                                        */
   /* ------------------------------------------------------------------------ */
 protected:
   void initializeNonLocal() override;
 
   void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override;
 
   void computeNonLocalStresses(GhostType ghost_type) override;
 
   void
   insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override;
 
   /// update the values of the non local internal
   void updateLocalInternal(ElementTypeMapReal & internal_flat,
                            GhostType ghost_type,
                            ElementKind kind) override;
 
   /// copy the results of the averaging in the materials
   void updateNonLocalInternal(ElementTypeMapReal & internal_flat,
                               GhostType ghost_type,
                               ElementKind kind) override;
 
   /* ------------------------------------------------------------------------ */
   /* Data Accessor inherited members                                          */
   /* ------------------------------------------------------------------------ */
 public:
   UInt getNbData(const Array<Element> & elements,
                  const SynchronizationTag & tag) const override;
 
   void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
                 const SynchronizationTag & tag) const override;
 
   void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
                   const SynchronizationTag & tag) override;
 
   UInt getNbData(const Array<UInt> & dofs,
                  const SynchronizationTag & tag) const override;
 
   void packData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
                 const SynchronizationTag & tag) const override;
 
   void unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
                   const SynchronizationTag & tag) override;
 
 protected:
   void
   splitElementByMaterial(const Array<Element> & elements,
                          std::vector<Array<Element>> & elements_per_mat) const;
 
   template <typename Operation>
   void splitByMaterial(const Array<Element> & elements, Operation && op) const;
 
   /* ------------------------------------------------------------------------ */
   /* Mesh Event Handler inherited members                                     */
   /* ------------------------------------------------------------------------ */
 protected:
   void onNodesAdded(const Array<UInt> & nodes_list,
                     const NewNodesEvent & event) override;
   void onNodesRemoved(const Array<UInt> & element_list,
                       const Array<UInt> & new_numbering,
                       const RemovedNodesEvent & event) override;
   void onElementsAdded(const Array<Element> & element_list,
                        const NewElementsEvent & event) override;
   void onElementsRemoved(const Array<Element> & element_list,
                          const ElementTypeMapArray<UInt> & new_numbering,
                          const RemovedElementsEvent & event) override;
   void onElementsChanged(const Array<Element> & /*unused*/,
                          const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const ChangedElementsEvent & /*unused*/) override{};
 
   /* ------------------------------------------------------------------------ */
   /* Dumpable interface (kept for convenience) and dumper relative functions  */
   /* ------------------------------------------------------------------------ */
 public:
   virtual void onDump();
 
   //! decide wether a field is a material internal or not
   bool isInternal(const std::string & field_name,
                   ElementKind element_kind);
   //! give the amount of data per element
   virtual ElementTypeMap<UInt>
   getInternalDataPerElem(const std::string & field_name,
                          ElementKind kind);
 
   //! flatten a given material internal field
   ElementTypeMapArray<Real> &
   flattenInternal(const std::string & field_name, ElementKind kind,
                   GhostType ghost_type = _not_ghost);
   //! flatten all the registered material internals
   void flattenAllRegisteredInternals(ElementKind kind);
 
   std::shared_ptr<dumpers::Field>
   createNodalFieldReal(const std::string & field_name,
                        const std::string & group_name,
                        bool padding_flag) override;
 
   std::shared_ptr<dumpers::Field>
   createNodalFieldBool(const std::string & field_name,
                        const std::string & group_name,
                        bool padding_flag) override;
 
   std::shared_ptr<dumpers::Field>
   createElementalField(const std::string & field_name,
                        const std::string & group_name, bool padding_flag,
                        UInt spatial_dimension,
                        ElementKind kind) override;
 
   virtual void dump(const std::string & dumper_name);
-
   virtual void dump(const std::string & dumper_name, UInt step);
-
   virtual void dump(const std::string & dumper_name, Real time, UInt step);
 
   void dump() override;
-
-  virtual void dump(UInt step);
-
-  virtual void dump(Real time, UInt step);
+  void dump(UInt step) override;
+  void dump(Real time, UInt step) override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// return the dimension of the system space
   AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt);
 
   /// set the value of the time step
   void setTimeStep(Real time_step, const ID & solver_id = "") override;
 
   /// get the value of the conversion from forces/ mass to acceleration
   AKANTU_GET_MACRO(F_M2A, f_m2a, Real);
 
   /// set the value of the conversion from forces/ mass to acceleration
   AKANTU_SET_MACRO(F_M2A, f_m2a, Real);
 
   /// get the SolidMechanicsModel::displacement array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement);
   /// get the SolidMechanicsModel::displacement array
   AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement);
 
   /// get the SolidMechanicsModel::previous_displacement array
   AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement);
 
   /// get the SolidMechanicsModel::current_position array
   const Array<Real> & getCurrentPosition();
 
   /// get  the SolidMechanicsModel::displacement_increment  array
   AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment);
   /// get  the SolidMechanicsModel::displacement_increment  array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment);
 
   /// get the lumped SolidMechanicsModel::mass array
   AKANTU_GET_MACRO_DEREF_PTR(Mass, mass);
 
   /// get the SolidMechanicsModel::velocity array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity);
   /// get the SolidMechanicsModel::velocity array
   AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity);
 
   /// get    the    SolidMechanicsModel::acceleration   array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration);
   /// get    the    SolidMechanicsModel::acceleration   array
   AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration);
 
   /// get the SolidMechanicsModel::external_force array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force);
   /// get the SolidMechanicsModel::external_force array
   AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force);
 
   /// get the SolidMechanicsModel::force array (external forces)
   [[deprecated("Use getExternalForce instead of this function")]] Array<Real> &
   getForce() {
     return getExternalForce();
   }
 
   /// get the SolidMechanicsModel::internal_force array (internal forces)
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force);
   /// get the SolidMechanicsModel::internal_force array (internal forces)
   AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force);
 
   /// get the SolidMechanicsModel::blocked_dofs array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs);
   /// get the SolidMechanicsModel::blocked_dofs array
   AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs);
 
   /// get an iterable on the materials
   inline decltype(auto) getMaterials();
 
   /// get an iterable on the materials
   inline decltype(auto) getMaterials() const;
 
   /// get a particular material (by numerical material index)
   inline Material & getMaterial(UInt mat_index);
 
   /// get a particular material (by numerical material index)
   inline const Material & getMaterial(UInt mat_index) const;
 
   /// get a particular material (by material name)
   inline Material & getMaterial(const std::string & name);
 
   /// get a particular material (by material name)
   inline const Material & getMaterial(const std::string & name) const;
 
   /// get a particular material id from is name
   inline UInt getMaterialIndex(const std::string & name) const;
 
   /// give the number of materials
   inline UInt getNbMaterials() const { return materials.size(); }
 
   /// give the material internal index from its id
   Int getInternalIndexFromID(const ID & id) const;
 
   /// compute the stable time step
   Real getStableTimeStep();
 
   /// get the energies
   Real getEnergy(const std::string & energy_id);
 
   /// compute the energy for energy
   Real getEnergy(const std::string & energy_id, ElementType type,
                  UInt index);
 
   AKANTU_GET_MACRO(MaterialByElement, material_index,
                    const ElementTypeMapArray<UInt> &);
   AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering,
                    const ElementTypeMapArray<UInt> &);
 
   /// vectors containing local material element index for each global element
   /// index
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index,
                                          UInt);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering,
                                          material_local_numbering, UInt);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering,
                                    material_local_numbering, UInt);
 
   AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector,
                              MaterialSelector &);
   void setMaterialSelector(std::shared_ptr<MaterialSelector> material_selector) {
     this->material_selector = std::move(material_selector);
   }
 
   /// Access the non_local_manager interface
   AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &);
 
   /// get the FEEngine object to integrate or interpolate on the boundary
   FEEngine & getFEEngineBoundary(const ID & name = "") override;
 
 protected:
   /// compute the stable time step
   Real getStableTimeStep(GhostType ghost_type);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   /// release version of the displacement array
   UInt displacement_release{0};
 
   /// release version of the current_position array
   UInt current_position_release{0};
 
   /// Check if materials need to recompute the mass array
   bool need_to_reassemble_lumped_mass{true};
   /// Check if materials need to recompute the mass matrix
   bool need_to_reassemble_mass{true};
 
   /// mapping between material name and material internal id
   std::map<std::string, UInt> materials_names_to_id;
 
 protected:
   /// conversion coefficient form force/mass to acceleration
   Real f_m2a{1.0};
 
   /// displacements array
   std::unique_ptr<Array<Real>> displacement;
 
   /// displacements array at the previous time step (used in finite deformation)
   std::unique_ptr<Array<Real>> previous_displacement;
 
   /// increment of displacement
   std::unique_ptr<Array<Real>> displacement_increment;
 
   /// lumped mass array
   std::unique_ptr<Array<Real>> mass;
 
   /// velocities array
   std::unique_ptr<Array<Real>> velocity;
 
   /// accelerations array
   std::unique_ptr<Array<Real>> acceleration;
 
   /// external forces array
   std::unique_ptr<Array<Real>> external_force;
 
   /// internal forces array
   std::unique_ptr<Array<Real>> internal_force;
 
   /// array specifing if a degree of freedom is blocked or not
   std::unique_ptr<Array<bool>> blocked_dofs;
 
   /// array of current position used during update residual
   std::unique_ptr<Array<Real>> current_position;
 
   /// Arrays containing the material index for each element
   ElementTypeMapArray<UInt> material_index;
 
   /// Arrays containing the position in the element filter of the material
   /// (material's local numbering)
   ElementTypeMapArray<UInt> material_local_numbering;
 
   /// list of used materials
   std::vector<std::unique_ptr<Material>> materials;
 
   /// class defining of to choose a material
   std::shared_ptr<MaterialSelector> material_selector;
 
   using flatten_internal_map =
       std::map<std::pair<std::string, ElementKind>,
                std::unique_ptr<ElementTypeMapArray<Real>>>;
 
   /// map a registered internals to be flattened for dump purposes
   flatten_internal_map registered_internals;
 
   /// non local manager
   std::unique_ptr<NonLocalManager> non_local_manager;
 
   /// tells if the material are instantiated
   bool are_materials_instantiated{false};
 };
 
 /* -------------------------------------------------------------------------- */
 namespace BC {
   namespace Neumann {
     using FromStress = FromHigherDim;
     using FromTraction = FromSameDim;
   } // namespace Neumann
 } // namespace BC
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "material.hh"
 #include "parser.hh"
 
 #include "solid_mechanics_model_inline_impl.hh"
 #include "solid_mechanics_model_tmpl.hh"
 /* -------------------------------------------------------------------------- */
 
 #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_io.cc b/src/model/solid_mechanics/solid_mechanics_model_io.cc
index b7603d9cb..a2a87b0fc 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_io.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_io.cc
@@ -1,341 +1,338 @@
 /**
  * @file   solid_mechanics_model_io.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jul 09 2017
  * @date last modification: Sun Dec 03 2017
  *
  * @brief  Dumpable part of the SolidMechnicsModel
  *
  *
  * 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 "solid_mechanics_model.hh"
 
 #include "group_manager_inline_impl.hh"
 
 #include "dumpable_inline_impl.hh"
 #ifdef AKANTU_USE_IOHELPER
 #include "dumper_element_partition.hh"
 #include "dumper_elemental_field.hh"
 #include "dumper_field.hh"
 #include "dumper_homogenizing_field.hh"
 #include "dumper_internal_material_field.hh"
 #include "dumper_iohelper.hh"
 #include "dumper_material_padders.hh"
 #include "dumper_paraview.hh"
 #endif
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 bool SolidMechanicsModel::isInternal(const std::string & field_name,
                                      ElementKind element_kind) {
   /// check if at least one material contains field_id as an internal
   for (auto & material : materials) {
     bool is_internal = material->isInternal<Real>(field_name, element_kind);
     if (is_internal) {
       return true;
     }
   }
 
   return false;
 }
 
 /* -------------------------------------------------------------------------- */
 ElementTypeMap<UInt>
 SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name,
                                             ElementKind element_kind) {
 
   if (!(this->isInternal(field_name, element_kind))) {
     AKANTU_EXCEPTION("unknown internal " << field_name);
   }
 
   for (auto & material : materials) {
     if (material->isInternal<Real>(field_name, element_kind)) {
       return material->getInternalDataPerElem<Real>(field_name, element_kind);
     }
   }
 
   return ElementTypeMap<UInt>();
 }
 
 /* -------------------------------------------------------------------------- */
 ElementTypeMapArray<Real> &
 SolidMechanicsModel::flattenInternal(const std::string & field_name,
                                      ElementKind kind,
                                      const GhostType ghost_type) {
   auto key = std::make_pair(field_name, kind);
 
   ElementTypeMapArray<Real> * internal_flat;
 
   auto it = this->registered_internals.find(key);
   if (it == this->registered_internals.end()) {
     auto internal = std::make_unique<ElementTypeMapArray<Real>>(
         field_name, this->id, this->memory_id);
 
     internal_flat = internal.get();
     this->registered_internals[key] = std::move(internal);
   } else {
     internal_flat = it->second.get();
   }
 
   for (auto type :
        mesh.elementTypes(Model::spatial_dimension, ghost_type, kind)) {
     if (internal_flat->exists(type, ghost_type)) {
       auto & internal = (*internal_flat)(type, ghost_type);
       internal.resize(0);
     }
   }
 
   for (auto & material : materials) {
     if (material->isInternal<Real>(field_name, kind)) {
       material->flattenInternal(field_name, *internal_flat, ghost_type, kind);
     }
   }
 
   return *internal_flat;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::flattenAllRegisteredInternals(
     ElementKind kind) {
   ElementKind _kind;
   ID _id;
 
   for (auto & internal : this->registered_internals) {
     std::tie(_id, _kind) = internal.first;
     if (kind == _kind) {
       this->flattenInternal(_id, kind);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onDump() {
   this->flattenAllRegisteredInternals(_ek_regular);
 }
 
 /* -------------------------------------------------------------------------- */
 #ifdef AKANTU_USE_IOHELPER
 std::shared_ptr<dumpers::Field> SolidMechanicsModel::createElementalField(
     const std::string & field_name, const std::string & group_name,
     bool padding_flag, UInt spatial_dimension,
     ElementKind kind) {
 
   std::shared_ptr<dumpers::Field> field;
 
   if (field_name == "partitions") {
     field = mesh.createElementalField<UInt, dumpers::ElementPartitionField>(
         mesh.getConnectivities(), group_name, spatial_dimension, kind);
   } else if (field_name == "material_index") {
     field = mesh.createElementalField<UInt, Vector, dumpers::ElementalField>(
         material_index, group_name, spatial_dimension, kind);
   } else {
     // this copy of field_name is used to compute derivated data such as
     // strain and von mises stress that are based on grad_u and stress
     std::string field_name_copy(field_name);
 
     if (field_name == "strain" || field_name == "Green strain" ||
         field_name == "principal strain" ||
         field_name == "principal Green strain") {
       field_name_copy = "grad_u";
     } else if (field_name == "Von Mises stress") {
       field_name_copy = "stress";
     }
 
     bool is_internal = this->isInternal(field_name_copy, kind);
 
     if (is_internal) {
       auto nb_data_per_elem =
           this->getInternalDataPerElem(field_name_copy, kind);
       auto & internal_flat = this->flattenInternal(field_name_copy, kind);
 
       field = mesh.createElementalField<Real, dumpers::InternalMaterialField>(
           internal_flat, group_name, spatial_dimension, kind, nb_data_per_elem);
 
       std::unique_ptr<dumpers::ComputeFunctorInterface> func;
       if (field_name == "strain") {
         func = std::make_unique<dumpers::ComputeStrain<false>>(*this);
       } else if (field_name == "Von Mises stress") {
         func = std::make_unique<dumpers::ComputeVonMisesStress>(*this);
       } else if (field_name == "Green strain") {
         func = std::make_unique<dumpers::ComputeStrain<true>>(*this);
       } else if (field_name == "principal strain") {
         func = std::make_unique<dumpers::ComputePrincipalStrain<false>>(*this);
       } else if (field_name == "principal Green strain") {
         func = std::make_unique<dumpers::ComputePrincipalStrain<true>>(*this);
       }
 
       if (func) {
         field = dumpers::FieldComputeProxy::createFieldCompute(field,
                                                               std::move(func));
       }
       // treat the paddings
       if (padding_flag) {
         if (field_name == "stress") {
           if (spatial_dimension == 2) {
             auto foo = std::make_unique<dumpers::StressPadder<2>>(*this);
             field = dumpers::FieldComputeProxy::createFieldCompute(
                 field, std::move(foo));
           }
         } else if (field_name == "strain" || field_name == "Green strain") {
           if (spatial_dimension == 2) {
             auto foo = std::make_unique<dumpers::StrainPadder<2>>(*this);
             field = dumpers::FieldComputeProxy::createFieldCompute(
                 field, std::move(foo));
           }
         }
       }
 
       // homogenize the field
       auto foo = dumpers::HomogenizerProxy::createHomogenizer(*field);
 
       field =
           dumpers::FieldComputeProxy::createFieldCompute(field, std::move(foo));
     }
   }
   return field;
 }
 
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field>
 SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
                                           const std::string & group_name,
                                           bool padding_flag) {
 
   std::map<std::string, Array<Real> *> real_nodal_fields;
   real_nodal_fields["displacement"] = this->displacement.get();
   real_nodal_fields["mass"] = this->mass.get();
   real_nodal_fields["velocity"] = this->velocity.get();
   real_nodal_fields["acceleration"] = this->acceleration.get();
   real_nodal_fields["external_force"] = this->external_force.get();
   real_nodal_fields["internal_force"] = this->internal_force.get();
   real_nodal_fields["increment"] = this->displacement_increment.get();
 
   if (field_name == "force") {
     AKANTU_EXCEPTION("The 'force' field has been renamed in 'external_force'");
   } else if (field_name == "residual") {
     AKANTU_EXCEPTION(
         "The 'residual' field has been replaced by 'internal_force'");
   }
 
   std::shared_ptr<dumpers::Field> field;
   if (padding_flag) {
     field = this->mesh.createNodalField(real_nodal_fields[field_name],
                                         group_name, 3);
   } else {
     field =
         this->mesh.createNodalField(real_nodal_fields[field_name], group_name);
   }
 
   return field;
 }
 
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field> SolidMechanicsModel::createNodalFieldBool(
     const std::string & field_name, const std::string & group_name,
     __attribute__((unused)) bool padding_flag) {
 
   std::map<std::string, Array<bool> *> uint_nodal_fields;
   uint_nodal_fields["blocked_dofs"] = blocked_dofs.get();
 
   std::shared_ptr<dumpers::Field> field;
   field = mesh.createNodalField(uint_nodal_fields[field_name], group_name);
   return field;
 }
 /* -------------------------------------------------------------------------- */
 #else
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field>
 SolidMechanicsModel::createElementalField(const std::string &,
                                           const std::string &, bool,
                                           const UInt &, ElementKind) {
   return nullptr;
 }
 /* --------------------------------------------------------------------------
  */
 std::shaed_ptr<dumpers::Field>
 SolidMechanicsModel::createNodalFieldReal(const std::string &,
                                           const std::string &, bool) {
   return nullptr;
 }
 
 /* --------------------------------------------------------------------------
  */
 std::shared_ptr<dumpers::Field>
 SolidMechanicsModel::createNodalFieldBool(const std::string &,
                                           const std::string &, bool) {
   return nullptr;
 }
 
 #endif
-/* --------------------------------------------------------------------------
- */
+/* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::dump(const std::string & dumper_name) {
   this->onDump();
   EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
   mesh.dump(dumper_name);
 }
 
-/* --------------------------------------------------------------------------
- */
+/* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) {
   this->onDump();
   EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
   mesh.dump(dumper_name, step);
 }
 
-/* -------------------------------------------------------------------------
- */
+/* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::dump(const std::string & dumper_name, Real time,
                                UInt step) {
   this->onDump();
   EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
   mesh.dump(dumper_name, time, step);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::dump() {
   this->onDump();
   EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
   mesh.dump();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::dump(UInt step) {
   this->onDump();
   EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
   mesh.dump(step);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::dump(Real time, UInt step) {
   this->onDump();
   EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
   mesh.dump(time, step);
 }
 
 } // namespace akantu
diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc
index 6c59ef6ba..cbc68ab76 100644
--- a/src/model/structural_mechanics/structural_mechanics_model.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model.cc
@@ -1,501 +1,507 @@
 /**
  * @file   structural_mechanics_model.cc
  *
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Damien Spielmann <damien.spielmann@epfl.ch>
  *
  * @date creation: Fri Jul 15 2011
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  Model implementation for Structural Mechanics elements
  *
  *
  * 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 "structural_mechanics_model.hh"
 #include "dof_manager.hh"
 #include "integrator_gauss.hh"
 #include "mesh.hh"
 #include "shape_structural.hh"
 #include "sparse_matrix.hh"
 #include "time_step_solver.hh"
 /* -------------------------------------------------------------------------- */
 #ifdef AKANTU_USE_IOHELPER
 #include "dumpable_inline_impl.hh"
 #include "dumper_elemental_field.hh"
 #include "dumper_iohelper_paraview.hh"
 #include "group_manager_inline_impl.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 #include "structural_element_bernoulli_beam_2.hh"
 #include "structural_element_bernoulli_beam_3.hh"
 #include "structural_element_kirchhoff_shell.hh"
 /* -------------------------------------------------------------------------- */
 //#include "structural_mechanics_model_inline_impl.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 inline UInt StructuralMechanicsModel::getNbDegreeOfFreedom(ElementType type) {
   UInt ndof = 0;
 #define GET_(type) ndof = ElementClass<type>::getNbDegreeOfFreedom()
   AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_, _ek_structural);
 #undef GET_
 
   return ndof;
 }
 
 /* -------------------------------------------------------------------------- */
 StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim,
                                                    const ID & id,
                                                    const MemoryID & memory_id)
     : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id),
       time_step(NAN), f_m2a(1.0), stress("stress", id, memory_id),
       element_material("element_material", id, memory_id),
       set_ID("beam sets", id, memory_id),
       rotation_matrix("rotation_matices", id, memory_id) {
   AKANTU_DEBUG_IN();
 
   registerFEEngineObject<MyFEEngineType>("StructuralMechanicsFEEngine", mesh,
                                          spatial_dimension);
 
   if (spatial_dimension == 2) {
     nb_degree_of_freedom = 3;
   } else if (spatial_dimension == 3) {
     nb_degree_of_freedom = 6;
   } else {
     AKANTU_TO_IMPLEMENT();
   }
 
 #ifdef AKANTU_USE_IOHELPER
   this->mesh.registerDumper<DumperParaview>("structural_mechanics_model", id,
                                             true);
 #endif
   this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural);
 
   this->initDOFManager();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 StructuralMechanicsModel::~StructuralMechanicsModel() = default;
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) {
   // <<<< This is the SolidMechanicsModel implementation for future ref >>>>
   // material_index.initialize(mesh, _element_kind = _ek_not_defined,
   //                           _default_value = UInt(-1), _with_nb_element =
   //                           true);
   // material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined,
   //                                     _with_nb_element = true);
 
   // Model::initFullImpl(options);
 
   // // initialize pbc
   // if (this->pbc_pair.size() != 0)
   //   this->initPBC();
 
   // // initialize the materials
   // if (this->parser.getLastParsedFile() != "") {
   //   this->instantiateMaterials();
   // }
 
   // this->initMaterials();
   // this->initBC(*this, *displacement, *displacement_increment,
   // *external_force);
 
   // <<<< END >>>>
 
   Model::initFullImpl(options);
 
   // Initializing stresses
   ElementTypeMap<UInt> stress_components;
 
   /// TODO this is ugly af, maybe add a function to FEEngine
   for (auto && type : mesh.elementTypes(_spatial_dimension = _all_dimensions,
                       _element_kind = _ek_structural)) {
     UInt nb_components = 0;
 
 // Getting number of components for each element type
 #define GET_(type) nb_components = ElementClass<type>::getNbStressComponents()
     AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_);
 #undef GET_
 
     stress_components(nb_components, type);
   }
 
   stress.initialize(
       getFEEngine(), _spatial_dimension = _all_dimensions,
       _element_kind = _ek_structural,
       _nb_component = [&stress_components](ElementType type,
                                            GhostType /*unused*/) -> UInt {
         return stress_components(type);
       });
 }
 
 /* -------------------------------------------------------------------------- */
 
 void StructuralMechanicsModel::initFEEngineBoundary() {
   /// TODO: this function should not be reimplemented
   /// we're just avoiding a call to Model::initFEEngineBoundary()
 }
 
 /* -------------------------------------------------------------------------- */
 // void StructuralMechanicsModel::setTimeStep(Real time_step) {
 //   this->time_step = time_step;
 
 // #if defined(AKANTU_USE_IOHELPER)
 //   this->mesh.getDumper().setTimeStep(time_step);
 // #endif
 // }
 
 /* -------------------------------------------------------------------------- */
 /* Initialisation                                                             */
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::initSolver(
     TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) {
   AKANTU_DEBUG_IN();
 
   this->allocNodalField(displacement_rotation, nb_degree_of_freedom,
                         "displacement");
   this->allocNodalField(external_force, nb_degree_of_freedom, "external_force");
   this->allocNodalField(internal_force, nb_degree_of_freedom, "internal_force");
   this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs");
 
   auto & dof_manager = this->getDOFManager();
 
   if (not dof_manager.hasDOFs("displacement")) {
     dof_manager.registerDOFs("displacement", *displacement_rotation,
                              _dst_nodal);
     dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
   }
 
   if (time_step_solver_type == TimeStepSolverType::_dynamic ||
       time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
     this->allocNodalField(velocity, nb_degree_of_freedom, "velocity");
     this->allocNodalField(acceleration, nb_degree_of_freedom, "acceleration");
 
     if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
       dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
       dof_manager.registerDOFsDerivative("displacement", 2,
                                          *this->acceleration);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::initModel() {
   for (auto && type : mesh.elementTypes(_element_kind = _ek_structural)) {
     // computeRotationMatrix(type);
     element_material.alloc(mesh.getNbElement(type), 1, type);
   }
 
   getFEEngine().initShapeFunctions(_not_ghost);
   getFEEngine().initShapeFunctions(_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleStiffnessMatrix() {
   AKANTU_DEBUG_IN();
 
+  if (not getDOFManager().hasMatrix("K")) {
+    getDOFManager().getNewMatrix("K", getMatrixType("K"));
+  }
+
   getDOFManager().getMatrix("K").zero();
 
   for (const auto & type :
        mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) {
 #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix<type>();
 
     AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX);
 #undef ASSEMBLE_STIFFNESS_MATRIX
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::computeStresses() {
   AKANTU_DEBUG_IN();
 
   for (const auto & type :
        mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) {
 #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad<type>();
 
     AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD);
 #undef COMPUTE_STRESS_ON_QUAD
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::computeRotationMatrix(ElementType type) {
   Mesh & mesh = getFEEngine().getMesh();
 
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_element = mesh.getNbElement(type);
 
   if (!rotation_matrix.exists(type)) {
     rotation_matrix.alloc(nb_element,
                           nb_degree_of_freedom * nb_nodes_per_element *
                               nb_degree_of_freedom * nb_nodes_per_element,
                           type);
   } else {
     rotation_matrix(type).resize(nb_element);
   }
   rotation_matrix(type).zero();
 
   Array<Real> rotations(nb_element,
                         nb_degree_of_freedom * nb_degree_of_freedom);
   rotations.zero();
 
 #define COMPUTE_ROTATION_MATRIX(type) computeRotationMatrix<type>(rotations);
 
   AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX);
 #undef COMPUTE_ROTATION_MATRIX
 
   auto R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom);
   auto T_it =
       rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element,
                                   nb_degree_of_freedom * nb_nodes_per_element);
 
   for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) {
     auto & T = *T_it;
     auto & R = *R_it;
     for (UInt k = 0; k < nb_nodes_per_element; ++k) {
       for (UInt i = 0; i < nb_degree_of_freedom; ++i) {
         for (UInt j = 0; j < nb_degree_of_freedom; ++j) {
           T(k * nb_degree_of_freedom + i, k * nb_degree_of_freedom + j) =
               R(i, j);
         }
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field> StructuralMechanicsModel::createNodalFieldBool(
     const std::string & field_name, const std::string & group_name,
     __attribute__((unused)) bool padding_flag) {
 
   std::map<std::string, Array<bool> *> uint_nodal_fields;
   uint_nodal_fields["blocked_dofs"] = blocked_dofs;
 
   return mesh.createNodalField(uint_nodal_fields[field_name], group_name);
 }
 
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field>
 StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name,
                                                const std::string & group_name,
                                                bool padding_flag) {
 
   UInt n;
   if (spatial_dimension == 2) {
     n = 2;
   } else {
     n = 3;
   }
 
   UInt padding_size = 0;
   if (padding_flag) {
     padding_size = 3;
   }
 
   if (field_name == "displacement") {
     return mesh.createStridedNodalField(displacement_rotation, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "velocity") {
     return mesh.createStridedNodalField(velocity, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "acceleration") {
     return mesh.createStridedNodalField(acceleration, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "rotation") {
     return mesh.createStridedNodalField(displacement_rotation, group_name,
                                         nb_degree_of_freedom - n, n,
                                         padding_size);
   }
 
   if (field_name == "force") {
     return mesh.createStridedNodalField(external_force, group_name, n, 0,
                                         padding_size);
   }
   if (field_name == "external_force") {
     return mesh.createStridedNodalField(external_force, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "momentum") {
     return mesh.createStridedNodalField(
         external_force, group_name, nb_degree_of_freedom - n, n, padding_size);
   }
 
   if (field_name == "internal_force") {
     return mesh.createStridedNodalField(internal_force, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "internal_momentum") {
     return mesh.createStridedNodalField(
         internal_force, group_name, nb_degree_of_freedom - n, n, padding_size);
   }
 
   return nullptr;
 }
 
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field> StructuralMechanicsModel::createElementalField(
     const std::string & field_name, const std::string & group_name,
     bool /*unused*/, UInt spatial_dimension, ElementKind kind) {
 
   std::shared_ptr<dumpers::Field> field;
 
   if (field_name == "element_index_by_material") {
     field = mesh.createElementalField<UInt, Vector, dumpers::ElementalField>(
         field_name, group_name, spatial_dimension, kind);
   }
   return field;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Virtual methods from SolverCallback */
 /* -------------------------------------------------------------------------- */
 /// get the type of matrix needed
 MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) {
   return _symmetric;
 }
 
 /// callback to assemble a Matrix
 void StructuralMechanicsModel::assembleMatrix(const ID & id) {
   if (id == "K") {
     assembleStiffnessMatrix();
   } else if (id == "M") {
     assembleMassMatrix();
   }
 }
 
 /// callback to assemble a lumped Matrix
 void StructuralMechanicsModel::assembleLumpedMatrix(const ID & /*id*/) {}
 
 /// callback to assemble the residual StructuralMechanicsModel::(rhs)
 void StructuralMechanicsModel::assembleResidual() {
   AKANTU_DEBUG_IN();
 
   auto & dof_manager = getDOFManager();
 
-  internal_force->zero();
-  computeStresses();
   assembleInternalForce();
-  dof_manager.assembleToResidual("displacement", *internal_force, -1);
+
   dof_manager.assembleToResidual("displacement", *external_force, 1);
+  dof_manager.assembleToResidual("displacement", *internal_force, 1);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /* Virtual methods from Model                                                 */
 /* -------------------------------------------------------------------------- */
 /// get some default values for derived classes
 std::tuple<ID, TimeStepSolverType>
 StructuralMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
   switch (method) {
   case _static: {
     return std::make_tuple("static", TimeStepSolverType::_static);
   }
   case _implicit_dynamic: {
     return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
   }
   default:
     return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
   }
 }
 
 /* ------------------------------------------------------------------------ */
 ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions(
     const TimeStepSolverType & type) const {
   ModelSolverOptions options;
 
   switch (type) {
   case TimeStepSolverType::_static: {
-    options.non_linear_solver_type = NonLinearSolverType::_linear;
+    options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_pseudo_time;
     options.solution_type["displacement"] = IntegrationScheme::_not_defined;
     break;
   }
   case TimeStepSolverType::_dynamic: {
     options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_trapezoidal_rule_2;
     options.solution_type["displacement"] = IntegrationScheme::_displacement;
     break;
   }
   default:
     AKANTU_EXCEPTION(type << " is not a valid time step solver type");
   }
 
   return options;
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleInternalForce() {
+  internal_force->zero();
+  computeStresses();
+
   for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions,
-                   _element_kind = _ek_structural)) {
+                                     _element_kind = _ek_structural)) {
     assembleInternalForce(type, _not_ghost);
     // assembleInternalForce(type, _ghost);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleInternalForce(ElementType type,
                                                      GhostType gt) {
   auto & fem = getFEEngine();
   auto & sigma = stress(type, gt);
   auto ndof = getNbDegreeOfFreedom(type);
   auto nb_nodes = mesh.getNbNodesPerElement(type);
   auto ndof_per_elem = ndof * nb_nodes;
 
   Array<Real> BtSigma(fem.getNbIntegrationPoints(type) *
                           mesh.getNbElement(type),
                       ndof_per_elem, "BtSigma");
   fem.computeBtD(sigma, BtSigma, type, gt);
 
   Array<Real> intBtSigma(0, ndof_per_elem, "intBtSigma");
   fem.integrate(BtSigma, intBtSigma, ndof_per_elem, type, gt);
   BtSigma.resize(0);
 
   getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force,
-                                                   type, gt, 1);
+                                                   type, gt, -1);
 }
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu
diff --git a/src/model/structural_mechanics/structural_mechanics_model_inline_impl.hh b/src/model/structural_mechanics/structural_mechanics_model_inline_impl.hh
index b7f3da972..e6d61f94a 100644
--- a/src/model/structural_mechanics/structural_mechanics_model_inline_impl.hh
+++ b/src/model/structural_mechanics/structural_mechanics_model_inline_impl.hh
@@ -1,342 +1,323 @@
 /**
  * @file   structural_mechanics_model_inline_impl.hh
  *
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Damien Spielmann <damien.spielmann@epfl.ch>
  *
  * @date creation: Fri Jul 15 2011
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Implementation of inline functions of StructuralMechanicsModel
  *
  *
  * 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 "structural_mechanics_model.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_STRUCTURAL_MECHANICS_MODEL_INLINE_IMPL_HH_
 #define AKANTU_STRUCTURAL_MECHANICS_MODEL_INLINE_IMPL_HH_
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void StructuralMechanicsModel::computeTangentModuli(
     Array<Real> & /*tangent_moduli*/) {
   AKANTU_TO_IMPLEMENT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void StructuralMechanicsModel::assembleStiffnessMatrix() {
   AKANTU_DEBUG_IN();
 
   auto nb_element = getFEEngine().getMesh().getNbElement(type);
   auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   auto nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
   auto tangent_size = ElementClass<type>::getNbStressComponents();
 
   auto tangent_moduli = std::make_unique<Array<Real>>(
       nb_element * nb_quadrature_points, tangent_size * tangent_size,
       "tangent_stiffness_matrix");
   computeTangentModuli<type>(*tangent_moduli);
 
   /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
   UInt bt_d_b_size = nb_degree_of_freedom * nb_nodes_per_element;
 
   auto bt_d_b = std::make_unique<Array<Real>>(
       nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B");
 
   const auto & b = getFEEngine().getShapesDerivatives(type);
 
   Matrix<Real> BtD(bt_d_b_size, tangent_size);
 
   for (auto && tuple :
        zip(make_view(b, tangent_size, bt_d_b_size),
            make_view(*tangent_moduli, tangent_size, tangent_size),
            make_view(*bt_d_b, bt_d_b_size, bt_d_b_size))) {
     auto & B = std::get<0>(tuple);
     auto & D = std::get<1>(tuple);
     auto & BtDB = std::get<2>(tuple);
     BtD.mul<true, false>(B, D);
     BtDB.template mul<false, false>(BtD, B);
   }
 
   /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
   auto int_bt_d_b = std::make_unique<Array<Real>>(
       nb_element, bt_d_b_size * bt_d_b_size, "int_B^t*D*B");
 
   getFEEngine().integrate(*bt_d_b, *int_bt_d_b, bt_d_b_size * bt_d_b_size,
                           type);
 
   getDOFManager().assembleElementalMatricesToMatrix(
       "K", "displacement", *int_bt_d_b, type, _not_ghost, _symmetric);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void StructuralMechanicsModel::computeStressOnQuad() {
   AKANTU_DEBUG_IN();
 
   auto & sigma = stress(type, _not_ghost);
 
   auto nb_element = mesh.getNbElement(type);
   auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   auto nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
   auto tangent_size = ElementClass<type>::getNbStressComponents();
 
   auto tangent_moduli = std::make_unique<Array<Real>>(
       nb_element * nb_quadrature_points, tangent_size * tangent_size,
       "tangent_stiffness_matrix");
 
   computeTangentModuli<type>(*tangent_moduli);
 
   /// compute DB
   auto d_b_size = nb_degree_of_freedom * nb_nodes_per_element;
 
   auto d_b = std::make_unique<Array<Real>>(nb_element * nb_quadrature_points,
                                            d_b_size * tangent_size, "D*B");
 
   const auto & b = getFEEngine().getShapesDerivatives(type);
 
   auto B = b.begin(tangent_size, d_b_size);
   auto D = tangent_moduli->begin(tangent_size, tangent_size);
   auto D_B = d_b->begin(tangent_size, d_b_size);
 
   for (UInt e = 0; e < nb_element; ++e) {
     for (UInt q = 0; q < nb_quadrature_points; ++q, ++B, ++D, ++D_B) {
       D_B->template mul<false, false>(*D, *B);
     }
   }
 
   /// compute DBu
   D_B = d_b->begin(tangent_size, d_b_size);
   auto DBu = sigma.begin(tangent_size);
 
   Array<Real> u_el(0, d_b_size);
   FEEngine::extractNodalToElementField(mesh, *displacement_rotation, u_el,
                                        type);
 
   auto ug = u_el.begin(d_b_size);
 
   // No need to rotate because B is post-multiplied
   for (UInt e = 0; e < nb_element; ++e, ++ug) {
     for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_B, ++DBu) {
       DBu->template mul<false>(*D_B, *ug);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void StructuralMechanicsModel::computeForcesByLocalTractionArray(
     const Array<Real> & tractions) {
   AKANTU_DEBUG_IN();
 
   auto nb_element = getFEEngine().getMesh().getNbElement(type);
   auto nb_nodes_per_element =
       getFEEngine().getMesh().getNbNodesPerElement(type);
   auto nb_quad = getFEEngine().getNbIntegrationPoints(type);
 
   // check dimension match
   AKANTU_DEBUG_ASSERT(
       Mesh::getSpatialDimension(type) == getFEEngine().getElementDimension(),
       "element type dimension does not match the dimension of boundaries : "
           << getFEEngine().getElementDimension()
           << " != " << Mesh::getSpatialDimension(type));
 
   // check size of the vector
   AKANTU_DEBUG_ASSERT(
       tractions.size() == nb_quad * nb_element,
       "the size of the vector should be the total number of quadrature points");
 
   // check number of components
   AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == nb_degree_of_freedom,
                       "the number of components should be the spatial "
                       "dimension of the problem");
 
   Array<Real> Ntbs(nb_element * nb_quad,
                    nb_degree_of_freedom * nb_nodes_per_element);
-  Array<Real> TtNtbs(nb_element * nb_quad,
-                     nb_degree_of_freedom * nb_nodes_per_element);
 
   auto & fem = getFEEngine();
   fem.computeNtb(tractions, Ntbs, type);
 
-  auto T_it =
-      rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element,
-                                  nb_degree_of_freedom * nb_nodes_per_element);
-  auto Ntb_it = Ntbs.begin(nb_degree_of_freedom * nb_nodes_per_element);
-  auto TtNtb_it = TtNtbs.begin(nb_degree_of_freedom * nb_nodes_per_element);
-
-  for (UInt e = 0; e < nb_element; ++e, ++T_it) {
-    const auto & T = *T_it;
-    for (UInt q = 0; q < nb_quad; ++q, ++Ntb_it, ++TtNtb_it) {
-      const auto & Ntb = *Ntb_it;
-      auto & TtNtb = *TtNtb_it;
-
-      // turn N^t tl back in the global referential
-      TtNtb.template mul<true>(T, Ntb);
-    }
-  }
-
   // allocate the vector that will contain the integrated values
   auto name = id + std::to_string(type) + ":integral_boundary";
   Array<Real> int_funct(nb_element, nb_degree_of_freedom * nb_nodes_per_element,
                         name);
 
   // do the integration
-  getFEEngine().integrate(TtNtbs, int_funct,
+  getFEEngine().integrate(Ntbs, int_funct,
                           nb_degree_of_freedom * nb_nodes_per_element, type);
 
   // assemble the result into force vector
   getDOFManager().assembleElementalArrayLocalArray(int_funct, *external_force,
                                                    type, _not_ghost, 1);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void StructuralMechanicsModel::computeForcesByGlobalTractionArray(
     const Array<Real> & traction_global) {
   AKANTU_DEBUG_IN();
 
   UInt nb_element = mesh.getNbElement(type);
   UInt nb_quad = getFEEngine().getNbIntegrationPoints(type);
   UInt nb_nodes_per_element =
       getFEEngine().getMesh().getNbNodesPerElement(type);
 
   std::stringstream name;
   name << id << ":structuralmechanics:imposed_linear_load";
   Array<Real> traction_local(nb_element * nb_quad, nb_degree_of_freedom,
                              name.str());
 
   auto T_it =
       rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element,
                                   nb_degree_of_freedom * nb_nodes_per_element);
 
   auto Te_it = traction_global.begin(nb_degree_of_freedom);
   auto te_it = traction_local.begin(nb_degree_of_freedom);
 
   Matrix<Real> R(nb_degree_of_freedom, nb_degree_of_freedom);
   for (UInt e = 0; e < nb_element; ++e, ++T_it) {
     const auto & T = *T_it;
     for (UInt i = 0; i < nb_degree_of_freedom; ++i) {
       for (UInt j = 0; j < nb_degree_of_freedom; ++j) {
         R(i, j) = T(i, j);
       }
     }
 
     for (UInt q = 0; q < nb_quad; ++q, ++Te_it, ++te_it) {
       const auto & Te = *Te_it;
       auto & te = *te_it;
       // turn the traction in the local referential
       te.template mul<false>(R, Te);
     }
   }
 
   computeForcesByLocalTractionArray<type>(traction_local);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * @param myf pointer  to a function that fills a  vector/tensor with respect to
  * passed coordinates
  */
 #if 0
 template <ElementType type>
 inline void StructuralMechanicsModel::computeForcesFromFunction(
     BoundaryFunction myf, BoundaryFunctionType function_type) {
   /** function type is
    ** _bft_forces : linear load is given
    ** _bft_stress : stress function is given -> Not already done for this kind
    *of model
    */
 
   std::stringstream name;
   name << id << ":structuralmechanics:imposed_linear_load";
   Array<Real> lin_load(0, nb_degree_of_freedom, name.str());
   name.zero();
 
   UInt offset = nb_degree_of_freedom;
 
   // prepare the loop over element types
   UInt nb_quad = getFEEngine().getNbIntegrationPoints(type);
   UInt nb_element = getFEEngine().getMesh().getNbElement(type);
 
   name.zero();
   name << id << ":structuralmechanics:quad_coords";
   Array<Real> quad_coords(nb_element * nb_quad, spatial_dimension,
                           "quad_coords");
 
   getFEEngineClass<MyFEEngineType>()
       .getShapeFunctions()
       .interpolateOnIntegrationPoints<type>(getFEEngine().getMesh().getNodes(),
                                             quad_coords, spatial_dimension);
   getFEEngineClass<MyFEEngineType>()
       .getShapeFunctions()
       .interpolateOnIntegrationPoints<type>(
           getFEEngine().getMesh().getNodes(), quad_coords, spatial_dimension,
           _not_ghost, empty_filter, true, 0, 1, 1);
   if (spatial_dimension == 3)
     getFEEngineClass<MyFEEngineType>()
         .getShapeFunctions()
         .interpolateOnIntegrationPoints<type>(
             getFEEngine().getMesh().getNodes(), quad_coords, spatial_dimension,
             _not_ghost, empty_filter, true, 0, 2, 2);
   lin_load.resize(nb_element * nb_quad);
   Real * imposed_val = lin_load.storage();
 
   /// sigma/load on each quadrature points
   Real * qcoord = quad_coords.storage();
   for (UInt el = 0; el < nb_element; ++el) {
     for (UInt q = 0; q < nb_quad; ++q) {
       myf(qcoord, imposed_val, NULL, 0);
       imposed_val += offset;
       qcoord += spatial_dimension;
     }
   }
 
   switch (function_type) {
   case _bft_traction_local:
     computeForcesByLocalTractionArray<type>(lin_load);
     break;
   case _bft_traction:
     computeForcesByGlobalTractionArray<type>(lin_load);
     break;
   default:
     break;
   }
 }
 #endif
 } // namespace akantu
 
 #endif /* AKANTU_STRUCTURAL_MECHANICS_MODEL_INLINE_IMPL_HH_ */
diff --git a/src/model/structural_mechanics/structural_mechanics_model_mass.cc b/src/model/structural_mechanics/structural_mechanics_model_mass.cc
index d7f15095e..9430c9c85 100644
--- a/src/model/structural_mechanics/structural_mechanics_model_mass.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model_mass.cc
@@ -1,78 +1,82 @@
 /**
  * @file   structural_mechanics_model_mass.cc
  *
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Mon Jul 07 2014
  * @date last modification: Fri Dec 15 2017
  *
  * @brief  function handling mass computation
  *
  *
  * Copyright (©) 2014-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_gauss.hh"
 #include "material.hh"
 #include "shape_structural.hh"
 #include "structural_mechanics_model.hh"
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 
 class ComputeRhoFunctorStruct {
 public:
   explicit ComputeRhoFunctorStruct(const StructuralMechanicsModel & model)
       : model(model){};
 
   void operator()(Matrix<Real> & rho, const Element & element) const {
     Real mat_rho = model.getMaterial(element).rho;
     rho.set(mat_rho);
   }
 
 private:
   const StructuralMechanicsModel & model;
 };
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleMassMatrix() {
   AKANTU_DEBUG_IN();
 
+  if (not getDOFManager().hasMatrix("M")) {
+    getDOFManager().getNewMatrix("M", getMatrixType("M"));
+  }
+
   assembleMassMatrix(_not_ghost);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleMassMatrix(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
   auto & fem = getFEEngineClass<MyFEEngineType>();
   ComputeRhoFunctorStruct compute_rho(*this);
 
   for (auto type :
        mesh.elementTypes(spatial_dimension, ghost_type, _ek_structural)) {
     fem.assembleFieldMatrix(compute_rho, "M", "displacement",
                             this->getDOFManager(), type, ghost_type);
   }
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc
index 8c42d1ff5..cd0a99bfa 100644
--- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc
+++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc
@@ -1,213 +1,258 @@
 /**
  * @file   test_structural_mechanics_model_bernoulli_beam_dynamics.cc
  *
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  *
  * @date creation: Mon Jul 07 2014
  * @date last modification: Wed Feb 03 2016
  *
  * @brief  Test for _bernouilli_beam in dynamic
  *
  *
  * Copyright (©) 2014-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 <fstream>
-#include <limits>
-
-/* -------------------------------------------------------------------------- */
-#include "aka_common.hh"
-#include "material.hh"
-#include "mesh.hh"
 #include "mesh_accessor.hh"
-#include "mesh_io.hh"
-#include "mesh_io_msh_struct.hh"
+#include "non_linear_solver_newton_raphson.hh"
 #include "structural_mechanics_model.hh"
-
+/* -------------------------------------------------------------------------- */
+#include <fstream>
 #include <iostream>
-using namespace akantu;
+#include <limits>
 /* -------------------------------------------------------------------------- */
 
-#define TYPE _bernoulli_beam_2
+using namespace akantu;
 
+/* -------------------------------------------------------------------------- */
 static Real analytical_solution(Real time, Real L, Real rho, Real E,
                                 __attribute__((unused)) Real A, Real I,
                                 Real F) {
   Real omega = M_PI * M_PI / L / L * sqrt(E * I / rho);
   Real sum = 0.;
   UInt i = 5;
   for (UInt n = 1; n <= i; n += 2) {
     sum += (1. - cos(n * n * omega * time)) / pow(n, 4);
   }
 
   return 2. * F * pow(L, 3) / pow(M_PI, 4) / E / I * sum;
 }
 
 // load
 const Real F = 0.5e4;
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
-
 int main(int argc, char * argv[]) {
   initialize(argc, argv);
-  Mesh beams(2);
+
+  /* ------------------------------------------------------------------------ */
+  UInt spatial_dimension = 2;
+  // auto method = _static;
+  auto method = _implicit_dynamic;
+  UInt nb_element = 100;
+  /* ------------------------------------------------------------------------ */
+
+  Mesh beams(spatial_dimension);
   debug::setDebugLevel(dblWarning);
-  const ElementType type = _bernoulli_beam_2;
 
-  /* --------------------------------------------------------------------------
-   */
+
+  ElementType type = _bernoulli_beam_2;
+  if (spatial_dimension == 3) {
+    type = _bernoulli_beam_3;
+  }
+  /* ------------------------------------------------------------------------ */
   // Mesh
-  UInt nb_element = 9;
+//  UInt nb_element = 8;
+  //nb_element *= 2;
+
   UInt nb_nodes = nb_element + 1;
-  Real total_length = 10.;
+  Real total_length = 2.;
   Real length = total_length / nb_element;
-  Real heigth = 1.;
 
   MeshAccessor mesh_accessor(beams);
   auto & nodes = mesh_accessor.getNodes();
   nodes.resize(nb_nodes);
 
   beams.addConnectivityType(type);
   auto & connectivity = mesh_accessor.getConnectivity(type);
   connectivity.resize(nb_element);
 
   beams.initNormals();
   auto & normals =
-      mesh_accessor.getData<Real>("extra_normal", type, _not_ghost, 2);
+      mesh_accessor.getData<Real>("extra_normal", type, _not_ghost, spatial_dimension);
   normals.resize(nb_element);
 
   for (UInt i = 0; i < nb_nodes; ++i) {
-    nodes(i, 0) = i * length;
-    nodes(i, 1) = 0;
+    nodes(i, _x) = i * length;
+    nodes(i, _y) = 0;
+    if (spatial_dimension == 3) {
+      nodes(i, _z) = 0;
+    }
   }
 
   for (UInt i = 0; i < nb_element; ++i) {
     connectivity(i, 0) = i;
     connectivity(i, 1) = i + 1;
+
+    if (spatial_dimension == 3) {
+      normals(i, _x) = 0;
+      normals(i, _y) = 0;
+      normals(i, _z) = 1;
+    }
   }
 
   mesh_accessor.makeReady();
 
-  /* --------------------------------------------------------------------------
-   */
+  /* ------------------------------------------------------------------------ */
   // Materials
-
   StructuralMechanicsModel model(beams);
 
   StructuralMaterial mat1;
   mat1.E = 120e6;
   mat1.rho = 1000;
-  mat1.A = heigth;
-  mat1.I = heigth * heigth * heigth / 12;
+
+  mat1.A = 0.003;
+
+  // Real heigth = .3;
+  // Real width  = .1;
+  mat1.I = 0.000225;
+  mat1.Iz = 0.000225;
+  mat1.Iy = 0.000025;
+  mat1.GJ = 0.0000790216;
+
   model.addMaterial(mat1);
 
-  /* --------------------------------------------------------------------------
-   */
+  /* ------------------------------------------------------------------------ */
   // Forces
-  // model.initFull();
-  model.initFull(_analysis_method = _implicit_dynamic);
+  //model.initFull(_analysis_method = _static);
+  model.initFull(_analysis_method = method);
 
   const auto & position = beams.getNodes();
   auto & forces = model.getExternalForce();
   auto & displacement = model.getDisplacement();
   auto & boundary = model.getBlockedDOFs();
 
-  UInt node_to_print = -1;
-
   forces.zero();
   displacement.zero();
   //  boundary.zero();
   // model.getElementMaterial(type)(i,0) = 0;
   // model.getElementMaterial(type)(i,0) = 1;
+  UInt node_to_print = -1;
   for (UInt i = 0; i < nb_element; ++i) {
     model.getElementMaterial(type)(i, 0) = 0;
   }
 
   for (UInt n = 0; n < nb_nodes; ++n) {
     Real x = position(n, _x);
     //    Real y = position(n, 1);
 
-    if (Math::are_float_equal(x, total_length / 2.)) {
+    if (Math::are_float_equal(x, total_length/2)) {
       forces(n, _y) = F;
       node_to_print = n;
     }
   }
 
   std::ofstream pos;
   pos.open("position.csv");
   if (not pos.good()) {
     AKANTU_ERROR("Cannot open file \"position.csv\"");
   }
 
   pos << "id,time,position,solution" << std::endl;
-  // model.computeForcesFromFunction<type>(load, _bft_traction)
   /* --------------------------------------------------------------------------
    */
   // Boundary conditions
   // true ~ displacement is blocked
-  boundary(0, 0) = true;
-  boundary(0, 1) = true;
-  boundary(nb_nodes - 1, 1) = true;
-  /* --------------------------------------------------------------------------
-   */
-  // "Solve"
+  boundary(0, _x) = true;
+  boundary(0, _y) = true;
 
+  boundary(nb_nodes - 1, _y) = true;
+
+  if(spatial_dimension == 3) {
+    boundary(0, _z) = true;
+    boundary(nb_nodes - 1, _z) = true;
+  }
+
+  /* ------------------------------------------------------------------------ */
+  // "Solve"
   Real time = 0;
   // model.assembleStiffnessMatrix();
   // model.assembleMass();
   model.addDumpFieldVector("displacement");
   model.addDumpFieldVector("force");
-  model.addDumpFieldVector("velocity");
-  model.addDumpFieldVector("acceleration");
+  model.addDumpFieldVector("internal_force");
+  if (method == _implicit_dynamic) {
+    model.addDumpFieldVector("velocity");
+    model.addDumpFieldVector("acceleration");
+  }
   model.dump();
+
   Real time_step = 1e-4;
   model.setTimeStep(time_step);
 
-  std::cout << "Time  |   Mid-Span Displacement" << std::endl;
+  model.assembleMatrix("K");
+  model.getDOFManager().getMatrix("K").saveMatrix("K.mtx");
 
-  /// time loop
-  for (UInt s = 1; time < 0.64; ++s) {
+  model.assembleMatrix("M");
+  model.getDOFManager().getMatrix("M").saveMatrix("M.mtx");
+
+  auto & solver = model.getNonLinearSolver();
+  solver.set("max_iterations", 100);
+  solver.set("threshold", 1e-8);
+  solver.set("convergence_type", SolveConvergenceCriteria::_solution);
 
+  if (method == _static) {
     model.solveStep();
+    model.dump();
+    return 0;
+  }
 
+  std::cout << "Time  |   Mid-Span Displacement" << std::endl;
+  /// time loop
+  for (UInt s = 1; time < 0.64; ++s) {
+    try {
+      model.solveStep();
+    } catch (debug::NLSNotConvergedException & e) {
+      std::cerr << "Did not converge after " << e.niter << " error: " << e.error
+                << " > " << e.threshold << std::endl;
+      return 1;
+    }
     pos << s << "," << time << "," << displacement(node_to_print, 1) << ","
         << analytical_solution(s * time_step, total_length, mat1.rho, mat1.E,
                                mat1.A, mat1.I, F)
         << std::endl;
     //    pos << s << "," << time << "," << displacement(node_to_print, 1) <<
     //    "," << analytical_solution(s*time_step) << std::endl;
 
+    Int n_iter = solver.get("nb_iterations");
+
     time += time_step;
-    if (s % 100 == 0)
-      std::cout << time << "  |   " << displacement(node_to_print, 1)
+    if (s % 100 == 0) {
+      std::cout << std::setw(5) << time << "  |   "
+                << std::setw(10) <<displacement(node_to_print, 1) << " |  " << n_iter
                 << std::endl;
-    model.dump();
-  }
-
-  model.getDOFManager().getMatrix("K").saveMatrix("K.mtx");
-  model.getDOFManager().getMatrix("M").saveMatrix("M.mtx");
+      model.dump(time, s);
+    }
 
+  }
 
   pos.close();
 
-  finalize();
-
   return EXIT_SUCCESS;
 }