diff --git a/src/fe_engine/fe_engine.hh b/src/fe_engine/fe_engine.hh
index 1c0516a8f..859d605c3 100644
--- a/src/fe_engine/fe_engine.hh
+++ b/src/fe_engine/fe_engine.hh
@@ -1,393 +1,393 @@
 /**
  * @file   fe_engine.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  FEM class
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_memory.hh"
 #include "element_type_map.hh"
 #include "mesh_events.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_FE_ENGINE_HH__
 #define __AKANTU_FE_ENGINE_HH__
 
 namespace akantu {
 class Mesh;
 class Integrator;
 class ShapeFunctions;
 class DOFManager;
 class Element;
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 /* -------------------------------------------------------------------------- */
 
 /**
  * The  generic  FEEngine class  derived  in  a  FEEngineTemplate class
  * containing  the
  * shape functions and the integration method
  */
 class FEEngine : protected Memory, public MeshEventHandler {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   FEEngine(Mesh & mesh, UInt spatial_dimension = _all_dimensions,
            const ID & id = "fem", MemoryID memory_id = 0);
 
   ~FEEngine() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// pre-compute all the shape functions, their derivatives and the jacobians
   virtual void
   initShapeFunctions(const GhostType & ghost_type = _not_ghost) = 0;
 
   /// extract the nodal values and store them per element
   template <typename T>
   static void extractNodalToElementField(
       const Mesh & mesh, const Array<T> & nodal_f, Array<T> & elemental_f,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter);
 
   /// filter a field
   template <typename T>
   static void
   filterElementalData(const Mesh & mesh, const Array<T> & quad_f,
                       Array<T> & filtered_f, const ElementType & type,
                       const GhostType & ghost_type = _not_ghost,
                       const Array<UInt> & filter_elements = empty_filter);
 
   /* ------------------------------------------------------------------------ */
   /* Integration method bridges                                               */
   /* ------------------------------------------------------------------------ */
   /// integrate f for all elements of type "type"
   virtual void
   integrate(const Array<Real> & f, Array<Real> & intf,
             UInt nb_degree_of_freedom, const ElementType & type,
             const GhostType & ghost_type = _not_ghost,
             const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// integrate a scalar value f on all elements of type "type"
   virtual Real
   integrate(const Array<Real> & f, const ElementType & type,
             const GhostType & ghost_type = _not_ghost,
             const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// integrate f for all integration points of type "type" but don't sum over
   /// all integration points
   virtual void integrateOnIntegrationPoints(
       const Array<Real> & f, Array<Real> & intf, UInt nb_degree_of_freedom,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// integrate one element scalar value on all elements of type "type"
   virtual Real integrate(const Vector<Real> & f, const ElementType & type,
                          UInt index,
                          const GhostType & ghost_type = _not_ghost) const = 0;
 
 /* ------------------------------------------------------------------------ */
 /* compatibility with old FEEngine fashion */
 /* ------------------------------------------------------------------------ */
 #ifndef SWIG
   /// get the number of integration points
   virtual UInt
   getNbIntegrationPoints(const ElementType & type,
                          const GhostType & ghost_type = _not_ghost) const = 0;
   /// get the precomputed shapes
   const virtual Array<Real> &
   getShapes(const ElementType & type, const GhostType & ghost_type = _not_ghost,
             UInt id = 0) const = 0;
 
   /// get the derivatives of shapes
   const virtual Array<Real> &
   getShapesDerivatives(const ElementType & type,
                        const GhostType & ghost_type = _not_ghost,
                        UInt id = 0) const = 0;
 
   /// get integration points
   const virtual Matrix<Real> &
   getIntegrationPoints(const ElementType & type,
                        const GhostType & ghost_type = _not_ghost) const = 0;
 #endif
   /* ------------------------------------------------------------------------ */
   /* Shape method bridges                                                     */
   /* ------------------------------------------------------------------------ */
   /// Compute the gradient nablauq on the integration points of an element type
   /// from nodal values u
   virtual void gradientOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq,
       const UInt nb_degree_of_freedom, const ElementType & type,
       const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// Interpolate a nodal field u at the integration points of an element type
   /// -> uq
   virtual void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// Interpolate a nodal field u at the integration points of many element
   /// types -> uq
   virtual void interpolateOnIntegrationPoints(
       const Array<Real> & u, ElementTypeMapArray<Real> & uq,
       const ElementTypeMapArray<UInt> * filter_elements = nullptr) const = 0;
 
   /// pre multiplies a tensor by the shapes derivaties
   virtual void
   computeBtD(const Array<Real> & Ds, Array<Real> & BtDs,
              const ElementType & type, const GhostType & ghost_type,
              const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// left and right  multiplies a tensor by the shapes derivaties
   virtual void
   computeBtDB(const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
               const ElementType & type, const GhostType & ghost_type,
               const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// Compute the interpolation point position in the global coordinates for
   /// many element types
   virtual void computeIntegrationPointsCoordinates(
       ElementTypeMapArray<Real> & integration_points_coordinates,
       const ElementTypeMapArray<UInt> * filter_elements = nullptr) const = 0;
 
   /// Compute the interpolation point position in the global coordinates for an
   /// element type
   virtual void computeIntegrationPointsCoordinates(
       Array<Real> & integration_points_coordinates, const ElementType & type,
       const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const = 0;
 
   /// Build pre-computed matrices for interpolation of field form integration
   /// points at other given positions (interpolation_points)
   virtual void initElementalFieldInterpolationFromIntegrationPoints(
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
       ElementTypeMapArray<Real> & integration_points_coordinates_inv_matrices,
       const ElementTypeMapArray<UInt> * element_filter) const = 0;
 
   /// interpolate field at given position (interpolation_points) from given
   /// values of this field at integration points (field)
   virtual void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & result, const GhostType ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const = 0;
 
   /// Interpolate field at given position from given values of this field at
   /// integration points (field)
   /// using matrices precomputed with
   /// initElementalFieldInterplationFromIntegrationPoints
   virtual void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> &
           interpolation_points_coordinates_matrices,
       const ElementTypeMapArray<Real> &
           integration_points_coordinates_inv_matrices,
       ElementTypeMapArray<Real> & result, const GhostType ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const = 0;
 
   /// interpolate on a phyiscal point inside an element
   virtual void interpolate(const Vector<Real> & real_coords,
                            const Matrix<Real> & nodal_values,
                            Vector<Real> & interpolated,
                            const Element & element) const = 0;
 
   /// compute the shape on a provided point
   virtual void
   computeShapes(const Vector<Real> & real_coords, UInt elem,
                 const ElementType & type, Vector<Real> & shapes,
                 const GhostType & ghost_type = _not_ghost) const = 0;
 
   /// compute the shape derivatives on a provided point
   virtual void
   computeShapeDerivatives(const Vector<Real> & real__coords, UInt element,
                           const ElementType & type,
                           Matrix<Real> & shape_derivatives,
                           const GhostType & ghost_type = _not_ghost) const = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Other methods                                                            */
   /* ------------------------------------------------------------------------ */
 
   /// pre-compute normals on integration points
   virtual void computeNormalsOnIntegrationPoints(
       const GhostType & ghost_type = _not_ghost) = 0;
 
   /// pre-compute normals on integration points
   virtual void computeNormalsOnIntegrationPoints(
       __attribute__((unused)) const Array<Real> & field,
       __attribute__((unused)) const GhostType & ghost_type = _not_ghost) {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// pre-compute normals on integration points
   virtual void computeNormalsOnIntegrationPoints(
       __attribute__((unused)) const Array<Real> & field,
       __attribute__((unused)) Array<Real> & normal,
       __attribute__((unused)) const ElementType & type,
       __attribute__((unused)) const GhostType & ghost_type = _not_ghost) const {
     AKANTU_TO_IMPLEMENT();
   }
 
-#ifdef AKANTU_STRUCTURAL_MECHANICS
-
-  /// assemble a field as a matrix for structural elements (ex. rho to mass
-  /// matrix)
-  virtual void assembleFieldMatrix(
-      __attribute__((unused)) const Array<Real> & field_1,
-      __attribute__((unused)) UInt nb_degree_of_freedom,
-      __attribute__((unused)) SparseMatrix & M,
-      __attribute__((unused)) Array<Real> * n,
-      __attribute__((unused)) ElementTypeMapArray<Real> & rotation_mat,
-      __attribute__((unused)) ElementType type,
-      __attribute__((unused)) const GhostType & ghost_type) const {
-
-    AKANTU_TO_IMPLEMENT();
-  }
-  /// compute shapes function in a matrix for structural elements
-  virtual void computeShapesMatrix(
-      __attribute__((unused)) const ElementType & type,
-      __attribute__((unused)) UInt nb_degree_of_freedom,
-      __attribute__((unused)) UInt nb_nodes_per_element,
-      __attribute__((unused)) Array<Real> * n, __attribute__((unused)) UInt id,
-      __attribute__((unused)) UInt degree_to_interpolate,
-      __attribute__((unused)) UInt degree_interpolated,
-      __attribute__((unused)) const bool sign,
-      __attribute__((unused)) const GhostType & ghost_type) const {
-
-    AKANTU_TO_IMPLEMENT();
-  }
-
-#endif
+// #ifdef AKANTU_STRUCTURAL_MECHANICS
+
+//   /// assemble a field as a matrix for structural elements (ex. rho to mass
+//   /// matrix)
+//   virtual void assembleFieldMatrix(
+//       __attribute__((unused)) const Array<Real> & field_1,
+//       __attribute__((unused)) UInt nb_degree_of_freedom,
+//       __attribute__((unused)) SparseMatrix & M,
+//       __attribute__((unused)) Array<Real> * n,
+//       __attribute__((unused)) ElementTypeMapArray<Real> & rotation_mat,
+//       __attribute__((unused)) ElementType type,
+//       __attribute__((unused)) const GhostType & ghost_type) const {
+
+//     AKANTU_TO_IMPLEMENT();
+//   }
+//   /// compute shapes function in a matrix for structural elements
+//   virtual void computeShapesMatrix(
+//       __attribute__((unused)) const ElementType & type,
+//       __attribute__((unused)) UInt nb_degree_of_freedom,
+//       __attribute__((unused)) UInt nb_nodes_per_element,
+//       __attribute__((unused)) Array<Real> * n, __attribute__((unused)) UInt id,
+//       __attribute__((unused)) UInt degree_to_interpolate,
+//       __attribute__((unused)) UInt degree_interpolated,
+//       __attribute__((unused)) const bool sign,
+//       __attribute__((unused)) const GhostType & ghost_type) const {
+
+//     AKANTU_TO_IMPLEMENT();
+//   }
+
+// #endif
 
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const;
 
 private:
   /// initialise the class
   void init();
 
   /* ------------------------------------------------------------------------ */
   /* Mesh Event Handler interface                                             */
   /* ------------------------------------------------------------------------ */
 public:
   void onNodesAdded(const Array<UInt> &, const NewNodesEvent &) final {}
   void onNodesRemoved(const Array<UInt> &, const Array<UInt> &,
                       const RemovedNodesEvent &) final {}
   void onElementsAdded(const Array<Element> &,
                        const NewElementsEvent &) override{};
   void onElementsRemoved(const Array<Element> &,
                          const ElementTypeMapArray<UInt> &,
                          const RemovedElementsEvent &) override {}
   void onElementsChanged(const Array<Element> &, const Array<Element> &,
                          const ElementTypeMapArray<UInt> &,
                          const ChangedElementsEvent &) override {}
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   using ElementTypesIteratorHelper =
       ElementTypeMapArray<Real, ElementType>::ElementTypesIteratorHelper;
 
   ElementTypesIteratorHelper elementTypes(UInt dim = _all_dimensions,
                                           GhostType ghost_type = _not_ghost,
                                           ElementKind kind = _ek_regular) const;
 
   /// get the dimension of the element handeled by this fe_engine object
   AKANTU_GET_MACRO(ElementDimension, element_dimension, UInt);
 
   /// get the mesh contained in the fem object
   AKANTU_GET_MACRO(Mesh, mesh, const Mesh &);
   /// get the mesh contained in the fem object
   AKANTU_GET_MACRO_NOT_CONST(Mesh, mesh, Mesh &);
 
   /// get the in-radius of an element
   static inline Real getElementInradius(const Matrix<Real> & coord,
                                         const ElementType & type);
 
   /// get the normals on integration points
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NormalsOnIntegrationPoints,
                                          normals_on_integration_points, Real);
 
   /// get cohesive element type for a given facet type
   static inline ElementType
   getCohesiveElementType(const ElementType & type_facet);
 
   /// get igfem element type for a given regular type
   static inline Vector<ElementType>
   getIGFEMElementTypes(const ElementType & type);
 
   /// get the interpolation element associated to an element type
   static inline InterpolationType
   getInterpolationType(const ElementType & el_type);
 
   /// get the shape function class (probably useless: see getShapeFunction in
   /// fe_engine_template.hh)
   virtual const ShapeFunctions & getShapeFunctionsInterface() const = 0;
   /// get the integrator class (probably useless: see getIntegrator in
   /// fe_engine_template.hh)
   virtual const Integrator & getIntegratorInterface() const = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// spatial dimension of the problem
   UInt element_dimension;
 
   /// the mesh on which all computation are made
   Mesh & mesh;
 
   /// normals at integration points
   ElementTypeMapArray<Real> normals_on_integration_points;
 };
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 
 /// standard output stream operator
 inline std::ostream & operator<<(std::ostream & stream,
                                  const FEEngine & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "fe_engine_inline_impl.cc"
 #include "fe_engine_template.hh"
 
 #endif /* __AKANTU_FE_ENGINE_HH__ */
diff --git a/src/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh
index 4325dbce6..e320f4127 100644
--- a/src/fe_engine/fe_engine_template.hh
+++ b/src/fe_engine/fe_engine_template.hh
@@ -1,407 +1,407 @@
 /**
  * @file   fe_engine_template.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Mon Jan 29 2018
  *
  * @brief  templated class that calls integration and shape objects
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_FE_ENGINE_TEMPLATE_HH__
 #define __AKANTU_FE_ENGINE_TEMPLATE_HH__
 
 /* -------------------------------------------------------------------------- */
 #include "fe_engine.hh"
 #include "integrator.hh"
 #include "shape_functions.hh"
 /* -------------------------------------------------------------------------- */
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 class DOFManager;
 namespace fe_engine {
   namespace details {
     template <ElementKind> struct AssembleLumpedTemplateHelper;
     template <ElementKind> struct AssembleFieldMatrixHelper;
   } // namespace details
 } // namespace fe_engine
 
 template <ElementKind, typename> struct AssembleFieldMatrixStructHelper;
 
 struct DefaultIntegrationOrderFunctor {
   template <ElementType type> static inline constexpr int getOrder() {
     return ElementClassProperty<type>::polynomial_degree;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind = _ek_regular,
           class IntegrationOrderFunctor = DefaultIntegrationOrderFunctor>
 class FEEngineTemplate : public FEEngine {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   using Integ = I<kind, IntegrationOrderFunctor>;
   using Shape = S<kind>;
 
   FEEngineTemplate(Mesh & mesh, UInt spatial_dimension = _all_dimensions,
                    ID id = "fem", MemoryID memory_id = 0);
 
   ~FEEngineTemplate() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// pre-compute all the shape functions, their derivatives and the jacobians
   void initShapeFunctions(const GhostType & ghost_type = _not_ghost) override;
   void initShapeFunctions(const Array<Real> & nodes,
                           const GhostType & ghost_type = _not_ghost);
 
   /* ------------------------------------------------------------------------ */
   /* Integration method bridges                                               */
   /* ------------------------------------------------------------------------ */
   /// integrate f for all elements of type "type"
   void
   integrate(const Array<Real> & f, Array<Real> & intf,
             UInt nb_degree_of_freedom, const ElementType & type,
             const GhostType & ghost_type = _not_ghost,
             const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// integrate a scalar value on all elements of type "type"
   Real
   integrate(const Array<Real> & f, const ElementType & type,
             const GhostType & ghost_type = _not_ghost,
             const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// integrate one element scalar value on all elements of type "type"
   Real integrate(const Vector<Real> & f, const ElementType & type, UInt index,
                  const GhostType & ghost_type = _not_ghost) const override;
 
   /// integrate partially around an integration point (@f$ intf_q = f_q * J_q *
   /// w_q @f$)
   void integrateOnIntegrationPoints(
       const Array<Real> & f, Array<Real> & intf, UInt nb_degree_of_freedom,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// interpolate on a phyiscal point inside an element
   void interpolate(const Vector<Real> & real_coords,
                    const Matrix<Real> & nodal_values,
                    Vector<Real> & interpolated,
                    const Element & element) const override;
 
   /// get the number of integration points
   UInt getNbIntegrationPoints(
       const ElementType & type,
       const GhostType & ghost_type = _not_ghost) const override;
 
   /// get shapes precomputed
   const Array<Real> & getShapes(const ElementType & type,
                                 const GhostType & ghost_type = _not_ghost,
                                 UInt id = 0) const override;
 
   /// get the derivatives of shapes
   const Array<Real> &
   getShapesDerivatives(const ElementType & type,
                        const GhostType & ghost_type = _not_ghost,
                        UInt id = 0) const override;
 
   /// get integration points
   const inline Matrix<Real> & getIntegrationPoints(
       const ElementType & type,
       const GhostType & ghost_type = _not_ghost) const override;
 
   /* ------------------------------------------------------------------------ */
   /* Shape method bridges                                                     */
   /* ------------------------------------------------------------------------ */
 
   /// compute the gradient of a nodal field on the integration points
   void gradientOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq,
       const UInt nb_degree_of_freedom, const ElementType & type,
       const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// interpolate a nodal field on the integration points
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       const ElementType & type, const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// interpolate a nodal field on the integration points given a
   /// by_element_type
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, ElementTypeMapArray<Real> & uq,
       const ElementTypeMapArray<UInt> * filter_elements =
           nullptr) const override;
 
   /// pre multiplies a tensor by the shapes derivaties
   void
   computeBtD(const Array<Real> & Ds, Array<Real> & BtDs,
              const ElementType & type, const GhostType & ghost_type,
              const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// left and right  multiplies a tensor by the shapes derivaties
   void computeBtDB(
       const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
       const ElementType & type, const GhostType & ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// compute the position of integration points given by an element_type_map
   /// from nodes position
   inline void computeIntegrationPointsCoordinates(
       ElementTypeMapArray<Real> & quadrature_points_coordinates,
       const ElementTypeMapArray<UInt> * filter_elements =
           nullptr) const override;
 
   /// compute the position of integration points from nodes position
   inline void computeIntegrationPointsCoordinates(
       Array<Real> & quadrature_points_coordinates, const ElementType & type,
       const GhostType & ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// interpolate field at given position (interpolation_points) from given
   /// values of this field at integration points (field)
   inline void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & result, const GhostType ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const override;
 
   /// Interpolate field at given position from given values of this field at
   /// integration points (field)
   /// using matrices precomputed with
   /// initElementalFieldInterplationFromIntegrationPoints
   inline void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> &
           interpolation_points_coordinates_matrices,
       const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       ElementTypeMapArray<Real> & result, const GhostType ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const override;
 
   /// Build pre-computed matrices for interpolation of field form integration
   /// points at other given positions (interpolation_points)
   inline void initElementalFieldInterpolationFromIntegrationPoints(
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
       ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       const ElementTypeMapArray<UInt> * element_filter =
           nullptr) const override;
 
   /// find natural coords from real coords provided an element
   void inverseMap(const Vector<Real> & real_coords, UInt element,
                   const ElementType & type, Vector<Real> & natural_coords,
                   const GhostType & ghost_type = _not_ghost) const;
 
   /// return true if the coordinates provided are inside the element, false
   /// otherwise
   inline bool contains(const Vector<Real> & real_coords, UInt element,
                        const ElementType & type,
                        const GhostType & ghost_type = _not_ghost) const;
 
   /// compute the shape on a provided point
   inline void
   computeShapes(const Vector<Real> & real_coords, UInt element,
                 const ElementType & type, Vector<Real> & shapes,
                 const GhostType & ghost_type = _not_ghost) const override;
 
   /// compute the shape derivatives on a provided point
   inline void computeShapeDerivatives(
       const Vector<Real> & real__coords, UInt element, const ElementType & type,
       Matrix<Real> & shape_derivatives,
       const GhostType & ghost_type = _not_ghost) const override;
 
   /* ------------------------------------------------------------------------ */
   /* Other methods                                                            */
   /* ------------------------------------------------------------------------ */
   /// pre-compute normals on integration points
   void computeNormalsOnIntegrationPoints(
       const GhostType & ghost_type = _not_ghost) override;
   void computeNormalsOnIntegrationPoints(
       const Array<Real> & field,
       const GhostType & ghost_type = _not_ghost) override;
   void computeNormalsOnIntegrationPoints(
       const Array<Real> & field, Array<Real> & normal, const ElementType & type,
       const GhostType & ghost_type = _not_ghost) const override;
   template <ElementType type>
   void computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                          Array<Real> & normal,
                                          const GhostType & ghost_type) const;
 
 private:
   // To avoid a weird full specialization of a method in a non specalized class
   void
   computeNormalsOnIntegrationPointsPoint1(const Array<Real> &,
                                           Array<Real> & normal,
                                           const GhostType & ghost_type) const;
 
 public:
   /// function to print the contain of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   /// assemble a field as a lumped matrix (ex. rho in lumped mass)
   template <class Functor>
   void assembleFieldLumped(const Functor & field_funct, const ID & matrix_id,
                            const ID & dof_id, DOFManager & dof_manager,
                            ElementType type,
                            const GhostType & ghost_type) const;
 
   /// assemble a field as a matrix (ex. rho to mass matrix)
   template <class Functor>
   void assembleFieldMatrix(const Functor & field_funct, const ID & matrix_id,
                            const ID & dof_id, DOFManager & dof_manager,
                            ElementType type,
                            const GhostType & ghost_type) const;
 
-#ifdef AKANTU_STRUCTURAL_MECHANICS
-  /// assemble a field as a matrix (ex. rho to mass matrix)
-  void assembleFieldMatrix(const Array<Real> & field_1,
-                           UInt nb_degree_of_freedom, SparseMatrix & M,
-                           Array<Real> * n,
-                           ElementTypeMapArray<Real> & rotation_mat,
-                           const ElementType & type,
-                           const GhostType & ghost_type = _not_ghost) const;
-
-  /// compute shapes function in a matrix for structural elements
-  void
-  computeShapesMatrix(const ElementType & type, UInt nb_degree_of_freedom,
-                      UInt nb_nodes_per_element, Array<Real> * n, UInt id,
-                      UInt degree_to_interpolate, UInt degree_interpolated,
-                      const bool sign,
-                      const GhostType & ghost_type = _not_ghost) const override;
-#endif
+// #ifdef AKANTU_STRUCTURAL_MECHANICS
+//   /// assemble a field as a matrix (ex. rho to mass matrix)
+//   void assembleFieldMatrix(const Array<Real> & field_1,
+//                            UInt nb_degree_of_freedom, SparseMatrix & M,
+//                            Array<Real> * n,
+//                            ElementTypeMapArray<Real> & rotation_mat,
+//                            const ElementType & type,
+//                            const GhostType & ghost_type = _not_ghost) const;
+
+//   /// compute shapes function in a matrix for structural elements
+//   void
+//   computeShapesMatrix(const ElementType & type, UInt nb_degree_of_freedom,
+//                       UInt nb_nodes_per_element, Array<Real> * n, UInt id,
+//                       UInt degree_to_interpolate, UInt degree_interpolated,
+//                       const bool sign,
+//                       const GhostType & ghost_type = _not_ghost) const override;
+// #endif
 
 private:
   friend struct fe_engine::details::AssembleLumpedTemplateHelper<kind>;
   friend struct fe_engine::details::AssembleFieldMatrixHelper<kind>;
   friend struct AssembleFieldMatrixStructHelper<kind, void>;
 
   /// templated function to compute the scaling to assemble a lumped matrix
   template <class Functor, ElementType type>
   void assembleFieldLumped(const Functor & field_funct, const ID & matrix_id,
                            const ID & dof_id, DOFManager & dof_manager,
                            const GhostType & ghost_type) const;
 
   /// @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j
   /// dV = \int \rho \varphi_i dV @f$
   template <ElementType type>
   void assembleLumpedRowSum(const Array<Real> & field, const ID & matrix_id,
                             const ID & dof_id, DOFManager & dof_manager,
                             const GhostType & ghost_type) const;
 
   /// @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
   template <ElementType type>
   void assembleLumpedDiagonalScaling(const Array<Real> & field,
                                      const ID & matrix_id, const ID & dof_id,
                                      DOFManager & dof_manager,
                                      const GhostType & ghost_type) const;
 
   /// assemble a field as a matrix (ex. rho to mass matrix)
   template <class Functor, ElementType type>
   void assembleFieldMatrix(const Functor & field_funct, const ID & matrix_id,
                            const ID & dof_id, DOFManager & dof_manager,
                            const GhostType & ghost_type) const;
 
 #ifdef AKANTU_STRUCTURAL_MECHANICS
 
   /// assemble a field as a matrix for structural elements (ex. rho to mass
   /// matrix)
   template <ElementType type>
   void assembleFieldMatrix(const Array<Real> & field_1,
                            UInt nb_degree_of_freedom, SparseMatrix & M,
                            Array<Real> * n,
                            ElementTypeMapArray<Real> & rotation_mat,
                            __attribute__((unused))
                            const GhostType & ghost_type) const;
 
 #endif
 
   /* ------------------------------------------------------------------------ */
   /* Mesh Event Handler interface                                             */
   /* ------------------------------------------------------------------------ */
 public:
   void onElementsAdded(const Array<Element> &,
                        const NewElementsEvent &) override;
   void onElementsRemoved(const Array<Element> &,
                          const ElementTypeMapArray<UInt> &,
                          const RemovedElementsEvent &) override;
   void onElementsChanged(const Array<Element> &, const Array<Element> &,
                          const ElementTypeMapArray<UInt> &,
                          const ChangedElementsEvent &) override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the shape class (probably useless: see getShapeFunction)
   const ShapeFunctions & getShapeFunctionsInterface() const override {
     return shape_functions;
   };
   /// get the shape class
   const Shape & getShapeFunctions() const { return shape_functions; };
 
   /// get the integrator class (probably useless: see getIntegrator)
   const Integrator & getIntegratorInterface() const override {
     return integrator;
   };
   /// get the integrator class
   const Integ & getIntegrator() const { return integrator; };
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   Integ integrator;
   Shape shape_functions;
 };
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "fe_engine_template_tmpl.hh"
 #include "fe_engine_template_tmpl_field.hh"
 /* -------------------------------------------------------------------------- */
 /* Shape Linked specialization                                                */
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 #include "fe_engine_template_tmpl_struct.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 /* Shape IGFEM specialization                                                 */
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_IGFEM)
 #include "fe_engine_template_tmpl_igfem.hh"
 #endif
 
 #endif /* __AKANTU_FE_ENGINE_TEMPLATE_HH__ */
diff --git a/src/fe_engine/fe_engine_template_tmpl_struct.hh b/src/fe_engine/fe_engine_template_tmpl_struct.hh
index e6318d1e5..16854018f 100644
--- a/src/fe_engine/fe_engine_template_tmpl_struct.hh
+++ b/src/fe_engine/fe_engine_template_tmpl_struct.hh
@@ -1,99 +1,99 @@
 /**
  * @file   fe_engine_template_tmpl_struct.hh
  *
  * @author Fabian Barras <fabian.barras@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: Tue Feb 20 2018
  *
  * @brief  Template implementation of FEEngineTemplate for Structural Element
  * Kinds
  *
  * @section LICENSE
  *
  * 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 "shape_structural.hh"
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind, typename = void>
 struct AssembleFieldMatrixStructHelper {};
 
 template <ElementKind kind>
 struct AssembleFieldMatrixStructHelper<
     kind, typename std::enable_if<kind == _ek_structural>::type> {
   template <template <ElementKind, class> class I,
             template <ElementKind> class S, ElementKind k, class IOF>
   static void call(const FEEngineTemplate<I, S, k, IOF> & fem,
                    const Array<Real> & field_1, UInt nb_degree_of_freedom,
                    SparseMatrix & M, Array<Real> * n,
                    ElementTypeMapArray<Real> & rotation_mat,
                    const ElementType & type, const GhostType & ghost_type) {
 #define ASSEMBLE_MATRIX(type)                                                  \
   fem.template assembleFieldMatrix<type>(field_1, nb_degree_of_freedom, M, n,  \
                                          rotation_mat, ghost_type)
 
     AKANTU_BOOST_KIND_ELEMENT_SWITCH(ASSEMBLE_MATRIX, _ek_structural);
 #undef ASSEMBLE_MATRIX
   }
 };
 
-template <template <ElementKind, class> class I, template <ElementKind> class S,
-          ElementKind kind, class IntegrationOrderFunctor>
-inline void
-FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::assembleFieldMatrix(
-    const Array<Real> & field_1, UInt nb_degree_of_freedom, SparseMatrix & M,
-    Array<Real> * n, ElementTypeMapArray<Real> & rotation_mat,
-    const ElementType & type, const GhostType & ghost_type) const {
-  AKANTU_DEBUG_IN();
+// template <template <ElementKind, class> class I, template <ElementKind> class S,
+//           ElementKind kind, class IntegrationOrderFunctor>
+// inline void
+// FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::assembleFieldMatrix(
+//     const Array<Real> & field_1, UInt nb_degree_of_freedom, SparseMatrix & M,
+//     Array<Real> * n, ElementTypeMapArray<Real> & rotation_mat,
+//     const ElementType & type, const GhostType & ghost_type) const {
+//   AKANTU_DEBUG_IN();
 
-  AssembleFieldMatrixStructHelper<kind>::template call(
-      *this, field_1, nb_degree_of_freedom, M, n, rotation_mat, type,
-      ghost_type);
+//   AssembleFieldMatrixStructHelper<kind>::template call(
+//       *this, field_1, nb_degree_of_freedom, M, n, rotation_mat, type,
+//       ghost_type);
 
-  AKANTU_DEBUG_OUT();
-}
-/* -------------------------------------------------------------------------- */
-template <template <ElementKind, class> class I, template <ElementKind> class S,
-          ElementKind kind, class IntegrationOrderFunctor>
-inline void
-FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapesMatrix(
-    const ElementType &, UInt, UInt, Array<Real> *, UInt, UInt, UInt,
-    const bool, const GhostType &) const {
+//   AKANTU_DEBUG_OUT();
+// }
+// /* -------------------------------------------------------------------------- */
+// template <template <ElementKind, class> class I, template <ElementKind> class S,
+//           ElementKind kind, class IntegrationOrderFunctor>
+// inline void
+// FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapesMatrix(
+//     const ElementType &, UInt, UInt, Array<Real> *, UInt, UInt, UInt,
+//     const bool, const GhostType &) const {
 
-  AKANTU_TO_IMPLEMENT();
-}
+//   AKANTU_TO_IMPLEMENT();
+// }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 inline void
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::assembleFieldMatrix(
     const Array<Real> &, UInt, SparseMatrix &, Array<Real> *,
     ElementTypeMapArray<Real> &, const GhostType &) const {
   AKANTU_TO_IMPLEMENT();
 }
 
 } // akantu
diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc
index 29392c28c..9ba97cefe 100644
--- a/src/model/structural_mechanics/structural_mechanics_model.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model.cc
@@ -1,483 +1,484 @@
 /**
  * @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 StucturalMechanics elements
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "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.cc"
 #endif
 /* -------------------------------------------------------------------------- */
 #include "structural_mechanics_model_inline_impl.cc"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 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>("paraview_all", 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, _all_ghost_types = true,
                     _nb_component = [&stress_components](
                         const ElementType & type, const GhostType &) -> 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) {
   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 (!dof_manager.hasDOFs("displacement")) {
     dof_manager.registerDOFs("displacement", *displacement_rotation,
                              _dst_nodal);
     dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
   }
 
   if (time_step_solver_type == _tsst_dynamic ||
       time_step_solver_type == _tsst_dynamic_lumped) {
     this->allocNodalField(velocity, spatial_dimension, "velocity");
     this->allocNodalField(acceleration, spatial_dimension, "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();
 
   getDOFManager().getMatrix("K").clear();
 
   for (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 (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(const 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).clear();
 
   Array<Real> rotations(nb_element,
                         nb_degree_of_freedom * nb_degree_of_freedom);
   rotations.clear();
 
 #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);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 dumper::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;
 
   dumper::Field * field = NULL;
   field = mesh.createNodalField(uint_nodal_fields[field_name], group_name);
   return field;
 }
 
 /* -------------------------------------------------------------------------- */
 dumper::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;
   }
 
   if (field_name == "displacement") {
     return mesh.createStridedNodalField(displacement_rotation, group_name, n, 0,
                                         padding_flag);
   }
 
   if (field_name == "rotation") {
     return mesh.createStridedNodalField(displacement_rotation, group_name,
                                         nb_degree_of_freedom - n, n,
                                         padding_flag);
   }
 
   if (field_name == "force") {
     return mesh.createStridedNodalField(external_force, group_name, n, 0,
                                         padding_flag);
   }
 
   if (field_name == "momentum") {
     return mesh.createStridedNodalField(
         external_force, group_name, nb_degree_of_freedom - n, n, padding_flag);
   }
 
   if (field_name == "internal_force") {
     return mesh.createStridedNodalField(internal_force, group_name, n, 0,
                                         padding_flag);
   }
 
   if (field_name == "internal_momentum") {
     return mesh.createStridedNodalField(
         internal_force, group_name, nb_degree_of_freedom - n, n, padding_flag);
   }
 
   return nullptr;
 }
 
 /* -------------------------------------------------------------------------- */
 dumper::Field * StructuralMechanicsModel::createElementalField(
     const std::string & field_name, const std::string & group_name, bool,
-    const ElementKind & kind, const std::string &) {
+    const UInt & spatial_dimension,
+    const ElementKind & kind) {
 
   dumper::Field * field = NULL;
 
   if (field_name == "element_index_by_material")
     field = mesh.createElementalField<UInt, Vector, dumper::ElementalField>(
-        field_name, group_name, this->spatial_dimension, kind);
+        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();
 }
 
 /// 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->clear();
   computeStresses();
   assembleInternalForce();
   dof_manager.assembleToResidual("displacement", *internal_force, -1);
   dof_manager.assembleToResidual("displacement", *external_force, 1);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /* Virtual methods from MeshEventHandler */
 /* -------------------------------------------------------------------------- */
 
 /// function to implement to react on  akantu::NewNodesEvent
 void StructuralMechanicsModel::onNodesAdded(const Array<UInt> & /*nodes_list*/,
                                             const NewNodesEvent & /*event*/) {}
 /// function to implement to react on  akantu::RemovedNodesEvent
 void StructuralMechanicsModel::onNodesRemoved(
     const Array<UInt> & /*nodes_list*/, const Array<UInt> & /*new_numbering*/,
     const RemovedNodesEvent & /*event*/) {}
 /// function to implement to react on  akantu::NewElementsEvent
 void StructuralMechanicsModel::onElementsAdded(
     const Array<Element> & /*elements_list*/,
     const NewElementsEvent & /*event*/) {}
 /// function to implement to react on  akantu::RemovedElementsEvent
 void StructuralMechanicsModel::onElementsRemoved(
     const Array<Element> & /*elements_list*/,
     const ElementTypeMapArray<UInt> & /*new_numbering*/,
     const RemovedElementsEvent & /*event*/) {}
 /// function to implement to react on  akantu::ChangedElementsEvent
 void StructuralMechanicsModel::onElementsChanged(
     const Array<Element> & /*old_elements_list*/,
     const Array<Element> & /*new_elements_list*/,
     const ElementTypeMapArray<UInt> & /*new_numbering*/,
     const ChangedElementsEvent & /*event*/) {}
 
 /* -------------------------------------------------------------------------- */
 /* 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", _tsst_static);
   }
   case _implicit_dynamic: {
     return std::make_tuple("implicit", _tsst_dynamic);
   }
   default:
     return std::make_tuple("unknown", _tsst_not_defined);
   }
 }
 
 /* ------------------------------------------------------------------------ */
 ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions(
     const TimeStepSolverType & type) const {
   ModelSolverOptions options;
 
   switch (type) {
   case _tsst_static: {
     options.non_linear_solver_type = _nls_linear;
     options.integration_scheme_type["displacement"] = _ist_pseudo_time;
     options.solution_type["displacement"] = IntegrationScheme::_not_defined;
     break;
   }
   default:
     AKANTU_EXCEPTION(type << " is not a valid time step solver type");
   }
 
   return options;
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleInternalForce() {
   for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions,
                    _element_kind = _ek_structural)) {
     assembleInternalForce(type, _not_ghost);
     // assembleInternalForce(type, _ghost);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleInternalForce(const 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);
 }
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu
diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh
index e0ef8f284..2cc61cce7 100644
--- a/src/model/structural_mechanics/structural_mechanics_model.hh
+++ b/src/model/structural_mechanics/structural_mechanics_model.hh
@@ -1,333 +1,333 @@
 /**
  * @file   structural_mechanics_model.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  Particular implementation of the structural elements in the
  * StructuralMechanicsModel
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_named_argument.hh"
 #include "boundary_condition.hh"
 #include "model.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__
 #define __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 class Material;
 class MaterialSelector;
 class DumperIOHelper;
 class NonLocalManager;
 template <ElementKind kind, class IntegrationOrderFunctor>
 class IntegratorGauss;
 template <ElementKind kind> class ShapeStructural;
 } // namespace akantu
 
 namespace akantu {
 
 struct StructuralMaterial {
   Real E{0};
   Real A{1};
   Real I{0};
   Real Iz{0};
   Real Iy{0};
   Real GJ{0};
   Real rho{0};
   Real t{0};
   Real nu{0};
 };
 
 class StructuralMechanicsModel : public Model {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   using MyFEEngineType =
       FEEngineTemplate<IntegratorGauss, ShapeStructural, _ek_structural>;
 
   StructuralMechanicsModel(Mesh & mesh,
                            UInt spatial_dimension = _all_dimensions,
                            const ID & id = "structural_mechanics_model",
                            const MemoryID & memory_id = 0);
 
   virtual ~StructuralMechanicsModel();
 
   /// Init full model
   void initFullImpl(const ModelOptions & options) override;
 
   /// Init boundary FEEngine
   void initFEEngineBoundary() override;
 
   /* ------------------------------------------------------------------------ */
   /* Virtual methods from SolverCallback */
   /* ------------------------------------------------------------------------ */
   /// get the type of matrix needed
   MatrixType getMatrixType(const ID &) override;
 
   /// callback to assemble a Matrix
   void assembleMatrix(const ID &) override;
 
   /// callback to assemble a lumped Matrix
   void assembleLumpedMatrix(const ID &) override;
 
   /// callback to assemble the residual (rhs)
   void assembleResidual() override;
 
   /* ------------------------------------------------------------------------ */
   /* Virtual methods from MeshEventHandler */
   /* ------------------------------------------------------------------------ */
 
   /// function to implement to react on  akantu::NewNodesEvent
   void onNodesAdded(const Array<UInt> & nodes_list,
                     const NewNodesEvent & event) override;
   /// function to implement to react on  akantu::RemovedNodesEvent
   void onNodesRemoved(const Array<UInt> & nodes_list,
                       const Array<UInt> & new_numbering,
                       const RemovedNodesEvent & event) override;
   /// function to implement to react on  akantu::NewElementsEvent
   void onElementsAdded(const Array<Element> & elements_list,
                        const NewElementsEvent & event) override;
   /// function to implement to react on  akantu::RemovedElementsEvent
   void onElementsRemoved(const Array<Element> & elements_list,
                          const ElementTypeMapArray<UInt> & new_numbering,
                          const RemovedElementsEvent & event) override;
   /// function to implement to react on  akantu::ChangedElementsEvent
   void onElementsChanged(const Array<Element> & old_elements_list,
                          const Array<Element> & new_elements_list,
                          const ElementTypeMapArray<UInt> & new_numbering,
                          const ChangedElementsEvent & event) override;
 
   /* ------------------------------------------------------------------------ */
   /* Virtual methods from Model */
   /* ------------------------------------------------------------------------ */
 protected:
   /// get some default values for derived classes
   std::tuple<ID, TimeStepSolverType>
   getDefaultSolverID(const AnalysisMethod & method) override;
 
   ModelSolverOptions
   getDefaultSolverOptions(const TimeStepSolverType & type) const override;
 
   UInt getNbDegreeOfFreedom(const ElementType & type) const;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
   void initSolver(TimeStepSolverType, NonLinearSolverType) override;
 
   /// initialize the model
   void initModel() override;
 
   /// compute the stresses per elements
   void computeStresses();
 
   /// compute the nodal forces
   void assembleInternalForce();
 
   /// compute the nodal forces for an element type
   void assembleInternalForce(const ElementType & type, GhostType gt);
 
   /// assemble the stiffness matrix
   void assembleStiffnessMatrix();
 
   /// assemble the mass matrix for consistent mass resolutions
   void assembleMass();
 
   /// TODO remove
   void computeRotationMatrix(const ElementType & type);
 
 protected:
   /// compute Rotation Matrices
   template <const ElementType type>
   void computeRotationMatrix(__attribute__((unused)) Array<Real> & rotations) {}
 
   /* ------------------------------------------------------------------------ */
   /* Mass (structural_mechanics_model_mass.cc) */
   /* ------------------------------------------------------------------------ */
 
   /// assemble the mass matrix for either _ghost or _not_ghost elements
   void assembleMass(GhostType ghost_type);
 
   /// computes rho
   void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type);
 
   /// finish the computation of residual to solve in increment
   void updateResidualInternal();
 
   /* ------------------------------------------------------------------------ */
 private:
   template <ElementType type> void assembleStiffnessMatrix();
   template <ElementType type> void assembleMass();
   template <ElementType type> void computeStressOnQuad();
   template <ElementType type>
   void computeTangentModuli(Array<Real> & tangent_moduli);
 
   /* ------------------------------------------------------------------------ */
   /* Dumpable interface                                                       */
   /* ------------------------------------------------------------------------ */
 public:
-  virtual dumper::Field * createNodalFieldReal(const std::string & field_name,
-                                               const std::string & group_name,
-                                               bool padding_flag);
-
-  virtual dumper::Field * createNodalFieldBool(const std::string & field_name,
-                                               const std::string & group_name,
-                                               bool padding_flag);
-
-  virtual dumper::Field *
-  createElementalField(const std::string & field_name,
-                       const std::string & group_name, bool padding_flag,
-                       const ElementKind & kind,
-                       const std::string & fe_engine_id = "");
+  dumper::Field * createNodalFieldReal(const std::string & field_name,
+                                       const std::string & group_name,
+                                       bool padding_flag) override;
+
+  dumper::Field * createNodalFieldBool(const std::string & field_name,
+                                       const std::string & group_name,
+                                       bool padding_flag) override;
+
+  dumper::Field * createElementalField(const std::string & field_name,
+                                       const std::string & group_name,
+                                       bool padding_flag,
+                                       const UInt & spatial_dimension,
+                                       const ElementKind & kind) override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// set the value of the time step
   // void setTimeStep(Real time_step, const ID & solver_id = "") override;
 
   /// return the dimension of the system space
   AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
 
   /// get the StructuralMechanicsModel::displacement vector
   AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array<Real> &);
 
   /// get the StructuralMechanicsModel::velocity vector
   AKANTU_GET_MACRO(Velocity, *velocity, Array<Real> &);
 
   /// get    the    StructuralMechanicsModel::acceleration    vector,   updated
   /// by
   /// StructuralMechanicsModel::updateAcceleration
   AKANTU_GET_MACRO(Acceleration, *acceleration, Array<Real> &);
 
   /// get the StructuralMechanicsModel::external_force vector
   AKANTU_GET_MACRO(ExternalForce, *external_force, Array<Real> &);
 
   /// get the StructuralMechanicsModel::internal_force vector (boundary forces)
   AKANTU_GET_MACRO(InternalForce, *internal_force, Array<Real> &);
 
   /// get the StructuralMechanicsModel::boundary vector
   AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<bool> &);
 
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(RotationMatrix, rotation_matrix, Real);
 
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
 
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt);
 
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Set_ID, set_ID, UInt);
 
   void addMaterial(StructuralMaterial & material) {
     materials.push_back(material);
   }
 
   const StructuralMaterial & getMaterial(const Element & element) const {
     return materials[element_material(element)];
   }
 
   /* ------------------------------------------------------------------------ */
   /* Boundaries (structural_mechanics_model_boundary.cc)                      */
   /* ------------------------------------------------------------------------ */
 public:
   /// Compute Linear load function set in global axis
   template <ElementType type>
   void computeForcesByGlobalTractionArray(const Array<Real> & tractions);
 
   /// Compute Linear load function set in local axis
   template <ElementType type>
   void computeForcesByLocalTractionArray(const Array<Real> & tractions);
 
   /// compute force vector from a function(x,y,momentum) that describe stresses
   // template <ElementType type>
   // void computeForcesFromFunction(BoundaryFunction in_function,
   //                                BoundaryFunctionType function_type);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   /// time step
   Real time_step;
 
   /// conversion coefficient form force/mass to acceleration
   Real f_m2a;
 
   /// displacements array
   Array<Real> * displacement_rotation{nullptr};
 
   /// velocities array
   Array<Real> * velocity{nullptr};
 
   /// accelerations array
   Array<Real> * acceleration{nullptr};
 
   /// forces array
   Array<Real> * internal_force{nullptr};
 
   /// forces array
   Array<Real> * external_force{nullptr};
 
   /// lumped mass array
   Array<Real> * mass{nullptr};
 
   /// boundaries array
   Array<bool> * blocked_dofs{nullptr};
 
   /// stress array
   ElementTypeMapArray<Real> stress;
 
   ElementTypeMapArray<UInt> element_material;
 
   // Define sets of beams
   ElementTypeMapArray<UInt> set_ID;
 
   /// number of degre of freedom
   UInt nb_degree_of_freedom;
 
   // Rotation matrix
   ElementTypeMapArray<Real> rotation_matrix;
 
   // /// analysis method check the list in akantu::AnalysisMethod
   // AnalysisMethod method;
 
   /// flag defining if the increment must be computed or not
   bool increment_flag;
 
   /* ------------------------------------------------------------------------ */
   std::vector<StructuralMaterial> materials;
 };
 
 } // namespace akantu
 
 #endif /* __AKANTU_STRUCTURAL_MECHANICS_MODEL_HH__ */
diff --git a/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc b/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc
index ebf76f302..8e7a048cd 100644
--- a/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc
@@ -1,371 +1,371 @@
 /**
  * @file   structural_mechanics_model_inline_impl.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: Tue Feb 20 2018
  *
  * @brief  Implementation of inline functions of StructuralMechanicsModel
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "structural_mechanics_model.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_STRUCTURAL_MECHANICS_MODEL_INLINE_IMPL_CC__
 #define __AKANTU_STRUCTURAL_MECHANICS_MODEL_INLINE_IMPL_CC__
 
 namespace akantu {
 /* -------------------------------------------------------------------------- */
 inline UInt
 StructuralMechanicsModel::getNbDegreeOfFreedom(const ElementType & type) const {
   UInt ndof = 0;
 #define GET_(type) ndof = ElementClass<type>::getNbDegreeOfFreedom()
   AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_, _ek_structural);
 #undef GET_
 
   return ndof;
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void StructuralMechanicsModel::computeTangentModuli(
     Array<Real> & /*tangent_moduli*/) {
   AKANTU_TO_IMPLEMENT();
 }
 } // namespace akantu
 
 #include "structural_element_bernoulli_beam_2.hh"
 #include "structural_element_bernoulli_beam_3.hh"
 #include "structural_element_kirchhoff_shell.hh"
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 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();
 
   Array<Real> & 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) {
+    const Array<Real> & /*tractions*/) {
   AKANTU_DEBUG_IN();
 
 #if 0
   UInt nb_element = getFEEngine().getMesh().getNbElement(type);
   UInt nb_nodes_per_element =
       getFEEngine().getMesh().getNbNodesPerElement(type);
   UInt 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> Nvoigt(nb_element * nb_quad, nb_degree_of_freedom *
                                                nb_degree_of_freedom *
                                                nb_nodes_per_element);
   transferNMatrixToSymVoigtNMatrix<type>(Nvoigt);
 
   Array<Real>::const_matrix_iterator N_it = Nvoigt.begin(
       nb_degree_of_freedom, nb_degree_of_freedom * nb_nodes_per_element);
   Array<Real>::const_matrix_iterator 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 =
       tractions.begin(nb_degree_of_freedom);
 
   Array<Real> funct(nb_element * nb_quad,
                     nb_degree_of_freedom * nb_nodes_per_element, 0.);
   Array<Real>::iterator<Vector<Real>> Fe_it =
       funct.begin(nb_degree_of_freedom * nb_nodes_per_element);
 
   Vector<Real> fe(nb_degree_of_freedom * nb_nodes_per_element);
   for (UInt e = 0; e < nb_element; ++e, ++T_it) {
     const Matrix<Real> & T = *T_it;
     for (UInt q = 0; q < nb_quad; ++q, ++N_it, ++te_it, ++Fe_it) {
       const Matrix<Real> & N = *N_it;
       const Vector<Real> & te = *te_it;
       Vector<Real> & Fe = *Fe_it;
 
       // compute N^t tl
       fe.mul<true>(N, te);
       // turn N^t tl back in the global referential
       Fe.mul<true>(T, fe);
     }
   }
 
   // allocate the vector that will contain the integrated values
   std::stringstream name;
   name << id << type << ":integral_boundary";
   Array<Real> int_funct(nb_element, nb_degree_of_freedom * nb_nodes_per_element,
                         name.str());
 
   // do the integration
   getFEEngine().integrate(funct, int_funct,
                           nb_degree_of_freedom * nb_nodes_per_element, type);
 
   // assemble the result into force vector
   getFEEngine().assembleArray(int_funct, *force_momentum,
                               dof_synchronizer->getLocalDOFEquationNumbers(),
                               nb_degree_of_freedom, type);
 #endif
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementType type>
 void StructuralMechanicsModel::computeForcesByGlobalTractionArray(
-    const Array<Real> & traction_global) {
+    const Array<Real> & /*traction_global*/) {
   AKANTU_DEBUG_IN();
 #if 0
   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());
 
   Array<Real>::const_matrix_iterator T_it =
       rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element,
                                   nb_degree_of_freedom * nb_nodes_per_element);
 
   Array<Real>::const_iterator<Vector<Real>> Te_it =
       traction_global.begin(nb_degree_of_freedom);
   Array<Real>::iterator<Vector<Real>> 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 Matrix<Real> & 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 Vector<Real> & Te = *Te_it;
       Vector<Real> & te = *te_it;
       // turn the traction in the local referential
       te.mul<false>(R, Te);
     }
   }
 
   computeForcesByLocalTractionArray<type>(traction_local);
 #endif
   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.clear();
 
   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.clear();
   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_CC__ */
diff --git a/src/model/time_step_solver.hh b/src/model/time_step_solver.hh
index 4437296ea..914486d93 100644
--- a/src/model/time_step_solver.hh
+++ b/src/model/time_step_solver.hh
@@ -1,138 +1,138 @@
 /**
  * @file   time_step_solver.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  This corresponding to the time step evolution solver
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 #include "aka_memory.hh"
 #include "integration_scheme.hh"
 #include "parameter_registry.hh"
 #include "solver_callback.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_TIME_STEP_SOLVER_HH__
 #define __AKANTU_TIME_STEP_SOLVER_HH__
 
 namespace akantu {
 class DOFManager;
 class NonLinearSolver;
 }
 
 namespace akantu {
 
 class TimeStepSolver : public Memory,
                        public ParameterRegistry,
                        public SolverCallback {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   TimeStepSolver(DOFManager & dof_manager, const TimeStepSolverType & type,
                  NonLinearSolver & non_linear_solver, const ID & id,
                  UInt memory_id);
   ~TimeStepSolver() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// solves on step
   virtual void solveStep(SolverCallback & solver_callback) = 0;
 
   /// register an integration scheme for a given dof
   virtual void
   setIntegrationScheme(const ID & dof_id, const IntegrationSchemeType & type,
                        IntegrationScheme::SolutionType solution_type =
                            IntegrationScheme::_not_defined) = 0;
 
   /// replies if a integration scheme has been set
   virtual bool hasIntegrationScheme(const ID & dof_id) const = 0;
   /* ------------------------------------------------------------------------ */
   /* Solver Callback interface                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// implementation of the SolverCallback::getMatrixType()
   MatrixType getMatrixType(const ID &) final { return _mt_not_defined; }
   /// implementation of the SolverCallback::predictor()
   void predictor() override;
   /// implementation of the SolverCallback::corrector()
   void corrector() override;
   /// implementation of the SolverCallback::assembleJacobian()
   void assembleMatrix(const ID & matrix_id) override;
   /// implementation of the SolverCallback::assembleJacobian()
   void assembleLumpedMatrix(const ID & matrix_id) final;
   /// implementation of the SolverCallback::assembleResidual()
   void assembleResidual() override;
   /// implementation of the SolverCallback::assembleResidual()
   void assembleResidual(const ID & residual_part) override;
 
   void beforeSolveStep() override;
   void afterSolveStep() override;
 
-  bool canSplitResidual() { return solver_callback->canSplitResidual(); }
+  bool canSplitResidual() override { return solver_callback->canSplitResidual(); }
   /* ------------------------------------------------------------------------ */
   /* Accessor                                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   AKANTU_GET_MACRO(TimeStep, time_step, Real);
   AKANTU_SET_MACRO(TimeStep, time_step, Real);
 
   AKANTU_GET_MACRO(NonLinearSolver, non_linear_solver, const NonLinearSolver &);
   AKANTU_GET_MACRO_NOT_CONST(NonLinearSolver, non_linear_solver,
                              NonLinearSolver &);
 
 protected:
   MatrixType getCommonMatrixType();
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// Underlying dof manager containing the dof to treat
   DOFManager & _dof_manager;
 
   /// Type of solver
   TimeStepSolverType type;
 
   /// The time step for this solver
   Real time_step;
 
   /// Temporary storage for solver callback
   SolverCallback * solver_callback;
 
   /// NonLinearSolver used by this tome step solver
   NonLinearSolver & non_linear_solver;
 
   /// List of required matrices
   std::map<std::string, MatrixType> needed_matrices;
 };
 
 } // akantu
 
 #endif /* __AKANTU_TIME_STEP_SOLVER_HH__ */