diff --git a/src/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh
index bbc36a5ff..fa51a0442 100644
--- a/src/fe_engine/fe_engine_template.hh
+++ b/src/fe_engine/fe_engine_template.hh
@@ -1,433 +1,426 @@
 /**
  * @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
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "fe_engine.hh"
 #include "integrator.hh"
 #include "shape_functions.hh"
 /* -------------------------------------------------------------------------- */
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
-
 /* -------------------------------------------------------------------------- */
 #ifndef AKANTU_FE_ENGINE_TEMPLATE_HH_
 #define AKANTU_FE_ENGINE_TEMPLATE_HH_
 
 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,
                    const ID & id = "fem");
 
   ~FEEngineTemplate() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// pre-compute all the shape functions, their derivatives and the jacobians
   void initShapeFunctions(GhostType ghost_type = _not_ghost) override;
   void initShapeFunctions(const Array<Real> & nodes,
                           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, ElementType type,
             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, ElementType type,
             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, ElementType type, UInt index,
                  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,
       ElementType type, 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(
-      ElementType type,
-      GhostType ghost_type = _not_ghost) const override;
+  UInt getNbIntegrationPoints(ElementType type,
+                              GhostType ghost_type = _not_ghost) const override;
 
   /// get shapes precomputed
   const Array<Real> & getShapes(ElementType type,
                                 GhostType ghost_type = _not_ghost,
                                 UInt id = 0) const override;
 
   /// get the derivatives of shapes
-  const Array<Real> &
-  getShapesDerivatives(ElementType type,
-                       GhostType ghost_type = _not_ghost,
-                       UInt id = 0) const override;
+  const Array<Real> & getShapesDerivatives(ElementType type,
+                                           GhostType ghost_type = _not_ghost,
+                                           UInt id = 0) const override;
 
   /// get integration points
-  const inline Matrix<Real> & getIntegrationPoints(
-      ElementType type,
-      GhostType ghost_type = _not_ghost) const override;
+  const inline Matrix<Real> &
+  getIntegrationPoints(ElementType type,
+                       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,
-      UInt nb_degree_of_freedom, ElementType type,
-      GhostType ghost_type = _not_ghost,
+      const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
+      ElementType type, 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,
       ElementType type, 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,
-             ElementType type, GhostType ghost_type,
+  computeBtD(const Array<Real> & Ds, Array<Real> & BtDs, ElementType type,
+             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,
       ElementType type, GhostType ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// left multiples a vector by the shape functions
   void
   computeNtb(const Array<Real> & bs, Array<Real> & Ntbs, ElementType type,
              GhostType ghost_type,
              const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// left and right  multiplies a tensor by the shapes
   void computeNtbN(
       const Array<Real> & bs, Array<Real> & NtbNs, ElementType type,
       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, ElementType type,
       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, 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, 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,
                   ElementType type, Vector<Real> & natural_coords,
                   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,
                        ElementType type,
                        GhostType ghost_type = _not_ghost) const;
 
   /// compute the shape on a provided point
-  inline void
-  computeShapes(const Vector<Real> & real_coords, UInt element,
-                ElementType type, Vector<Real> & shapes,
-                GhostType ghost_type = _not_ghost) const override;
+  inline void computeShapes(const Vector<Real> & real_coords, UInt element,
+                            ElementType type, Vector<Real> & shapes,
+                            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, ElementType type,
-      Matrix<Real> & shape_derivatives,
-      GhostType ghost_type = _not_ghost) const override;
+  inline void
+  computeShapeDerivatives(const Vector<Real> & real_coords, UInt element,
+                          ElementType type, Matrix<Real> & shape_derivatives,
+                          GhostType ghost_type = _not_ghost) const override;
 
   /* ------------------------------------------------------------------------ */
   /* Other methods                                                            */
   /* ------------------------------------------------------------------------ */
   /// pre-compute normals on integration points
-  void computeNormalsOnIntegrationPoints(
-      GhostType ghost_type = _not_ghost) override;
-  void computeNormalsOnIntegrationPoints(
-      const Array<Real> & field,
-      GhostType ghost_type = _not_ghost) override;
+  void
+  computeNormalsOnIntegrationPoints(GhostType ghost_type = _not_ghost) override;
+  void
+  computeNormalsOnIntegrationPoints(const Array<Real> & field,
+                                    GhostType ghost_type = _not_ghost) override;
   void computeNormalsOnIntegrationPoints(
       const Array<Real> & field, Array<Real> & normal, ElementType type,
       GhostType ghost_type = _not_ghost) const override;
   template <ElementType type>
   void computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                          Array<Real> & normal,
                                          GhostType ghost_type) const;
 
 private:
   // To avoid a weird full specialization of a method in a non specalized class
-  void
-  computeNormalsOnIntegrationPointsPoint1(const Array<Real> & /*unused*/,
-                                          Array<Real> & normal,
-                                          GhostType ghost_type) const;
+  void computeNormalsOnIntegrationPointsPoint1(const Array<Real> & /*unused*/,
+                                               Array<Real> & normal,
+                                               GhostType ghost_type) const;
 
 public:
   /// function to print the contain of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   void assembleFieldLumped(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       ElementType type, GhostType ghost_type) const override;
 
   /// assemble a field as a matrix (ex. rho to mass matrix)
   void assembleFieldMatrix(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       ElementType type, GhostType ghost_type) const 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,
   //                          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,
   //                          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,
   //                            ElementType type,
   //                            GhostType ghost_type = _not_ghost)
   //                            const;
 
   //   /// compute shapes function in a matrix for structural elements
   //   void
   //   computeShapesMatrix(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,
   //                       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 <ElementType type>
   void assembleFieldLumped(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       GhostType ghost_type) const;
 
   /// @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,
                             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,
                                      GhostType ghost_type) const;
 
   /// assemble a field as a matrix (ex. rho to mass matrix)
   template <ElementType type>
   void assembleFieldMatrix(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       GhostType ghost_type) const;
 
 #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))
-                           GhostType ghost_type) const;
+                           __attribute__((unused)) GhostType ghost_type) const;
 
 #endif
 
   /* ------------------------------------------------------------------------ */
   /* Mesh Event Handler interface                                             */
   /* ------------------------------------------------------------------------ */
 public:
   void onElementsAdded(const Array<Element> & /*new_elements*/,
                        const NewElementsEvent & /*unused*/) override;
   void onElementsRemoved(const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const RemovedElementsEvent & /*unused*/) override;
   void onElementsChanged(const Array<Element> & /*unused*/,
                          const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const ChangedElementsEvent & /*unused*/) 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.hh b/src/fe_engine/fe_engine_template_tmpl.hh
index 185c93dab..dead0cf5d 100644
--- a/src/fe_engine/fe_engine_template_tmpl.hh
+++ b/src/fe_engine/fe_engine_template_tmpl.hh
@@ -1,1456 +1,1444 @@
 /**
  * @file   fe_engine_template_tmpl.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Mauro Corrado <mauro.corrado@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Tue Feb 15 2011
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Template implementation of FEEngineTemplate
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 #include "dof_manager.hh"
 #include "fe_engine_template.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::FEEngineTemplate(
     Mesh & mesh, UInt spatial_dimension, const ID & id)
     : FEEngine(mesh, spatial_dimension, id),
       integrator(mesh, spatial_dimension, id),
       shape_functions(mesh, spatial_dimension, id) {}
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::~FEEngineTemplate() =
     default;
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GradientOnIntegrationPointsHelper {
       template <class S>
       static void call(const S & /*unused*/, Mesh & /*unused*/,
                        const Array<Real> & /*unused*/, Array<Real> & /*unused*/,
                        const UInt /*unused*/, ElementType /*unused*/,
-                       GhostType /*unused*/,
-                       const Array<UInt> & /*unused*/) {
+                       GhostType /*unused*/, const Array<UInt> & /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_GRADIENT(type)                                                 \
   if (element_dimension == ElementClass<type>::getSpatialDimension())          \
     shape_functions.template gradientOnIntegrationPoints<type>(                \
         u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER(kind)          \
   template <> struct GradientOnIntegrationPointsHelper<kind> {                 \
     template <class S>                                                         \
     static void call(const S & shape_functions, Mesh & mesh,                   \
                      const Array<Real> & u, Array<Real> & nablauq,             \
-                     const UInt nb_degree_of_freedom,                          \
-                     ElementType type, GhostType ghost_type,   \
+                     const UInt nb_degree_of_freedom, ElementType type,        \
+                     GhostType ghost_type,                                     \
                      const Array<UInt> & filter_elements) {                    \
       UInt element_dimension = mesh.getSpatialDimension(type);                 \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_GRADIENT, kind);                \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER,
         AKANTU_FE_ENGINE_LIST_GRADIENT_ON_INTEGRATION_POINTS)
 
 #undef AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER
 #undef COMPUTE_GRADIENT
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     gradientOnIntegrationPoints(const Array<Real> & u, Array<Real> & nablauq,
                                 const UInt nb_degree_of_freedom,
-                                ElementType type,
-                                GhostType ghost_type,
+                                ElementType type, GhostType ghost_type,
                                 const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
   UInt nb_points =
       shape_functions.getIntegrationPoints(type, ghost_type).cols();
 
 #ifndef AKANTU_NDEBUG
 
   UInt element_dimension = mesh.getSpatialDimension(type);
 
   AKANTU_DEBUG_ASSERT(u.size() == mesh.getNbNodes(),
                       "The vector u(" << u.getID()
                                       << ") has not the good size.");
   AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom,
                       "The vector u("
                           << u.getID()
                           << ") has not the good number of component.");
 
   AKANTU_DEBUG_ASSERT(
       nablauq.getNbComponent() == nb_degree_of_freedom * element_dimension,
       "The vector nablauq(" << nablauq.getID()
                             << ") has not the good number of component.");
 
 // AKANTU_DEBUG_ASSERT(nablauq.size() == nb_element * nb_points,
 //                  "The vector nablauq(" << nablauq.getID()
 //                  << ") has not the good size.");
 #endif
 
   nablauq.resize(nb_element * nb_points);
 
   fe_engine::details::GradientOnIntegrationPointsHelper<kind>::call(
       shape_functions, mesh, u, nablauq, nb_degree_of_freedom, type, ghost_type,
       filter_elements);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::initShapeFunctions(
     GhostType ghost_type) {
   initShapeFunctions(mesh.getNodes(), ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::initShapeFunctions(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) {
     integrator.initIntegrator(nodes, type, ghost_type);
     const auto & control_points = getIntegrationPoints(type, ghost_type);
     shape_functions.initShapeFunctions(nodes, control_points, type, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct IntegrateHelper {};
 
 #define INTEGRATE(type)                                                        \
   integrator.template integrate<type>(f, intf, nb_degree_of_freedom,           \
                                       ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_INTEGRATE_HELPER(kind)                               \
   template <> struct IntegrateHelper<kind> {                                   \
     template <class I>                                                         \
     static void call(const I & integrator, const Array<Real> & f,              \
                      Array<Real> & intf, UInt nb_degree_of_freedom,            \
-                     ElementType type, GhostType ghost_type,   \
+                     ElementType type, GhostType ghost_type,                   \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_HELPER)
 
 #undef AKANTU_SPECIALIZE_INTEGRATE_HELPER
 #undef INTEGRATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
     const Array<Real> & f, Array<Real> & intf, UInt nb_degree_of_freedom,
     ElementType type, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 #ifndef AKANTU_NDEBUG
 
   UInt nb_quadrature_points = getNbIntegrationPoints(type);
 
   AKANTU_DEBUG_ASSERT(f.size() == nb_element * nb_quadrature_points,
                       "The vector f(" << f.getID() << " size " << f.size()
                                       << ") has not the good size ("
                                       << nb_element << ").");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
                       "The vector intf("
                           << intf.getID()
                           << ") has not the good number of component.");
 #endif
 
   intf.resize(nb_element);
 
   fe_engine::details::IntegrateHelper<kind>::call(integrator, f, intf,
                                                   nb_degree_of_freedom, type,
                                                   ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct IntegrateScalarHelper {};
 
 #define INTEGRATE(type)                                                        \
   integral =                                                                   \
       integrator.template integrate<type>(f, ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER(kind)                        \
   template <> struct IntegrateScalarHelper<kind> {                             \
     template <class I>                                                         \
     static Real call(const I & integrator, const Array<Real> & f,              \
-                     ElementType type, GhostType ghost_type,   \
+                     ElementType type, GhostType ghost_type,                   \
                      const Array<UInt> & filter_elements) {                    \
       Real integral = 0.;                                                      \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
       return integral;                                                         \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER)
 
 #undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER
 #undef INTEGRATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 Real FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
-    const Array<Real> & f, ElementType type,
-    GhostType ghost_type, const Array<UInt> & filter_elements) const {
+    const Array<Real> & f, ElementType type, GhostType ghost_type,
+    const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
 #ifndef AKANTU_NDEBUG
   //   std::stringstream sstr; sstr << ghost_type;
   //   AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
   //                  "The vector " << nablauq.getID() << " is not taged " <<
   //                  ghost_type);
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 
   UInt nb_quadrature_points = getNbIntegrationPoints(type, ghost_type);
 
   AKANTU_DEBUG_ASSERT(
       f.size() == nb_element * nb_quadrature_points,
       "The vector f(" << f.getID() << ") has not the good size. (" << f.size()
                       << "!=" << nb_quadrature_points * nb_element << ")");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == 1,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
 #endif
 
   Real integral = fe_engine::details::IntegrateScalarHelper<kind>::call(
       integrator, f, type, ghost_type, filter_elements);
   AKANTU_DEBUG_OUT();
   return integral;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct IntegrateScalarOnOneElementHelper {};
 
 #define INTEGRATE(type)                                                        \
   res = integrator.template integrate<type>(f, index, ghost_type);
 
 #define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER(kind)         \
   template <> struct IntegrateScalarOnOneElementHelper<kind> {                 \
     template <class I>                                                         \
     static Real call(const I & integrator, const Vector<Real> & f,             \
-                     ElementType type, UInt index,                     \
-                     GhostType ghost_type) {                           \
+                     ElementType type, UInt index, GhostType ghost_type) {     \
       Real res = 0.;                                                           \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
       return res;                                                              \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(
         AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER)
 
 #undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER
 #undef INTEGRATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 Real FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate(
     const Vector<Real> & f, ElementType type, UInt index,
     GhostType ghost_type) const {
 
   Real res = fe_engine::details::IntegrateScalarOnOneElementHelper<kind>::call(
       integrator, f, type, index, ghost_type);
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct IntegrateOnIntegrationPointsHelper {};
 
 #define INTEGRATE(type)                                                        \
   integrator.template integrateOnIntegrationPoints<type>(                      \
       f, intf, nb_degree_of_freedom, ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER(kind)         \
   template <> struct IntegrateOnIntegrationPointsHelper<kind> {                \
     template <class I>                                                         \
     static void call(const I & integrator, const Array<Real> & f,              \
                      Array<Real> & intf, UInt nb_degree_of_freedom,            \
-                     ElementType type, GhostType ghost_type,   \
+                     ElementType type, GhostType ghost_type,                   \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind);                       \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(
         AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER)
 
 #undef AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER
 #undef INTEGRATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     integrateOnIntegrationPoints(const Array<Real> & f, Array<Real> & intf,
-                                 UInt nb_degree_of_freedom,
-                                 ElementType type,
+                                 UInt nb_degree_of_freedom, ElementType type,
                                  GhostType ghost_type,
                                  const Array<UInt> & filter_elements) const {
 
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
   UInt nb_quadrature_points = getNbIntegrationPoints(type);
 #ifndef AKANTU_NDEBUG
   //   std::stringstream sstr; sstr << ghost_type;
   //   AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
   //                  "The vector " << nablauq.getID() << " is not taged " <<
   //                  ghost_type);
 
   AKANTU_DEBUG_ASSERT(f.size() == nb_element * nb_quadrature_points,
                       "The vector f(" << f.getID() << " size " << f.size()
                                       << ") has not the good size ("
                                       << nb_element << ").");
   AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom,
                       "The vector f("
                           << f.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
                       "The vector intf("
                           << intf.getID()
                           << ") has not the good number of component.");
 #endif
 
   intf.resize(nb_element * nb_quadrature_points);
   fe_engine::details::IntegrateOnIntegrationPointsHelper<kind>::call(
       integrator, f, intf, nb_degree_of_freedom, type, ghost_type,
       filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct InterpolateOnIntegrationPointsHelper {
       template <class S>
       static void call(const S & /*unused*/, const Array<Real> & /*unused*/,
                        Array<Real> & /*unused*/, const UInt /*unused*/,
-                       ElementType /*unused*/,
-                       GhostType /*unused*/,
+                       ElementType /*unused*/, GhostType /*unused*/,
                        const Array<UInt> & /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define INTERPOLATE(type)                                                      \
   shape_functions.template interpolateOnIntegrationPoints<type>(               \
       u, uq, nb_degree_of_freedom, ghost_type, filter_elements);
 
 #define AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER(kind)       \
   template <> struct InterpolateOnIntegrationPointsHelper<kind> {              \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & u,         \
                      Array<Real> & uq, const UInt nb_degree_of_freedom,        \
-                     ElementType type, GhostType ghost_type,   \
+                     ElementType type, GhostType ghost_type,                   \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind);                     \
     }                                                                          \
   };
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER,
         AKANTU_FE_ENGINE_LIST_INTERPOLATE_ON_INTEGRATION_POINTS)
 
 #undef AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER
 #undef INTERPOLATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateOnIntegrationPoints(const Array<Real> & u, Array<Real> & uq,
                                    const UInt nb_degree_of_freedom,
-                                   ElementType type,
-                                   GhostType ghost_type,
+                                   ElementType type, GhostType ghost_type,
                                    const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_points =
       shape_functions.getIntegrationPoints(type, ghost_type).cols();
   UInt nb_element = mesh.getNbElement(type, ghost_type);
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 #ifndef AKANTU_NDEBUG
 
   AKANTU_DEBUG_ASSERT(u.size() == mesh.getNbNodes(),
                       "The vector u(" << u.getID()
                                       << ") has not the good size.");
   AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom,
                       "The vector u("
                           << u.getID()
                           << ") has not the good number of component.");
   AKANTU_DEBUG_ASSERT(uq.getNbComponent() == nb_degree_of_freedom,
                       "The vector uq("
                           << uq.getID()
                           << ") has not the good number of component.");
 #endif
 
   uq.resize(nb_element * nb_points);
 
   fe_engine::details::InterpolateOnIntegrationPointsHelper<kind>::call(
       shape_functions, u, uq, nb_degree_of_freedom, type, ghost_type,
       filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateOnIntegrationPoints(
         const Array<Real> & u, ElementTypeMapArray<Real> & uq,
         const ElementTypeMapArray<UInt> * filter_elements) const {
   AKANTU_DEBUG_IN();
 
   const Array<UInt> * filter = nullptr;
 
   for (auto ghost_type : ghost_types) {
     for (auto && type : uq.elementTypes(_all_dimensions, ghost_type, kind)) {
       UInt nb_quad_per_element = getNbIntegrationPoints(type, ghost_type);
 
       UInt nb_element = 0;
 
       if (filter_elements != nullptr) {
         filter = &((*filter_elements)(type, ghost_type));
         nb_element = filter->size();
       } else {
         filter = &empty_filter;
         nb_element = mesh.getNbElement(type, ghost_type);
       }
 
       UInt nb_tot_quad = nb_quad_per_element * nb_element;
 
       Array<Real> & quad = uq(type, ghost_type);
       quad.resize(nb_tot_quad);
 
       interpolateOnIntegrationPoints(u, quad, quad.getNbComponent(), type,
                                      ghost_type, *filter);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeBtDHelper {};
 
 #define COMPUTE_BTD(type)                                                      \
   shape_functions.template computeBtD<type>(Ds, BtDs, ghost_type,              \
                                             filter_elements);
 
 #define AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER(kind)                             \
   template <> struct ComputeBtDHelper<kind> {                                  \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & Ds,        \
-                     Array<Real> & BtDs, ElementType type,             \
-                     GhostType ghost_type,                             \
+                     Array<Real> & BtDs, ElementType type,                     \
+                     GhostType ghost_type,                                     \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_BTD, kind);                     \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_BtD_HELPER
 #undef COMPUTE_BTD
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeBtD(
     const Array<Real> & Ds, Array<Real> & BtDs, ElementType type,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   fe_engine::details::ComputeBtDHelper<kind>::call(
       shape_functions, Ds, BtDs, type, ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeBtDBHelper {};
 
 #define COMPUTE_BTDB(type)                                                     \
   shape_functions.template computeBtDB<type>(Ds, BtDBs, order_d, ghost_type,   \
                                              filter_elements);
 
 #define AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER(kind)                            \
   template <> struct ComputeBtDBHelper<kind> {                                 \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & Ds,        \
-                     Array<Real> & BtDBs, UInt order_d,                        \
-                     ElementType type, GhostType ghost_type,   \
+                     Array<Real> & BtDBs, UInt order_d, ElementType type,      \
+                     GhostType ghost_type,                                     \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_BTDB, kind);                    \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_BtDB_HELPER
 #undef COMPUTE_BTDB
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeBtDB(
-    const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
-    ElementType type, GhostType ghost_type,
-    const Array<UInt> & filter_elements) const {
+    const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d, ElementType type,
+    GhostType ghost_type, const Array<UInt> & filter_elements) const {
   fe_engine::details::ComputeBtDBHelper<kind>::call(
       shape_functions, Ds, BtDBs, order_d, type, ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeNtbNHelper {};
 
 #define COMPUTE_NtbN(type)                                                     \
   shape_functions.template computeNtbN<type>(bs, NtbNs, ghost_type,            \
                                              filter_elements);
 
 #define AKANTU_SPECIALIZE_COMPUTE_NtbN_HELPER(kind)                            \
   template <> struct ComputeNtbNHelper<kind> {                                 \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & bs,        \
                      Array<Real> & NtbNs, ElementType type,                    \
                      GhostType ghost_type,                                     \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_NtbN, kind);                    \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_NtbN_HELPER)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_NtbN_HELPER
 #undef COMPUTE_NtbN
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeNtbN(
     const Array<Real> & bs, Array<Real> & NtbNs, ElementType type,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   fe_engine::details::ComputeNtbNHelper<kind>::call(
       shape_functions, bs, NtbNs, type, ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeNtbHelper {};
 
 #define COMPUTE_Ntb(type)                                                      \
   shape_functions.template computeNtb<type>(bs, Ntbs, ghost_type,              \
                                             filter_elements);
 
 #define AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER(kind)                             \
   template <> struct ComputeNtbHelper<kind> {                                  \
     template <class S>                                                         \
     static void call(const S & shape_functions, const Array<Real> & bs,        \
-                     Array<Real> & Ntbs, ElementType type,             \
-                     GhostType ghost_type,                             \
+                     Array<Real> & Ntbs, ElementType type,                     \
+                     GhostType ghost_type,                                     \
                      const Array<UInt> & filter_elements) {                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_Ntb, kind);                     \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_Ntb_HELPER
 #undef COMPUTE_Ntb
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeNtb(
     const Array<Real> & bs, Array<Real> & Ntbs, ElementType type,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   fe_engine::details::ComputeNtbHelper<kind>::call(
       shape_functions, bs, Ntbs, type, ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeIntegrationPointsCoordinates(
         ElementTypeMapArray<Real> & quadrature_points_coordinates,
         const ElementTypeMapArray<UInt> * filter_elements) const {
 
   const Array<Real> & nodes_coordinates = mesh.getNodes();
 
   interpolateOnIntegrationPoints(
       nodes_coordinates, quadrature_points_coordinates, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeIntegrationPointsCoordinates(
         Array<Real> & quadrature_points_coordinates, ElementType type,
-        GhostType ghost_type,
-        const Array<UInt> & filter_elements) const {
+        GhostType ghost_type, const Array<UInt> & filter_elements) const {
 
   const Array<Real> & nodes_coordinates = mesh.getNodes();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
 
   interpolateOnIntegrationPoints(
       nodes_coordinates, quadrature_points_coordinates, spatial_dimension, type,
       ghost_type, filter_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     initElementalFieldInterpolationFromIntegrationPoints(
         const ElementTypeMapArray<Real> & interpolation_points_coordinates,
         ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
         ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
         const ElementTypeMapArray<UInt> * element_filter) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = this->mesh.getSpatialDimension();
 
   ElementTypeMapArray<Real> quadrature_points_coordinates(
       "quadrature_points_coordinates_for_interpolation", getID());
 
   quadrature_points_coordinates.initialize(*this,
                                            _nb_component = spatial_dimension);
 
   computeIntegrationPointsCoordinates(quadrature_points_coordinates,
                                       element_filter);
   shape_functions.initElementalFieldInterpolationFromIntegrationPoints(
       interpolation_points_coordinates,
       interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, quadrature_points_coordinates,
       element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateElementalFieldFromIntegrationPoints(
         const ElementTypeMapArray<Real> & field,
         const ElementTypeMapArray<Real> & interpolation_points_coordinates,
         ElementTypeMapArray<Real> & result, GhostType ghost_type,
         const ElementTypeMapArray<UInt> * element_filter) const {
 
   ElementTypeMapArray<Real> interpolation_points_coordinates_matrices(
       "interpolation_points_coordinates_matrices", id);
   ElementTypeMapArray<Real> quad_points_coordinates_inv_matrices(
       "quad_points_coordinates_inv_matrices", id);
 
   initElementalFieldInterpolationFromIntegrationPoints(
       interpolation_points_coordinates,
       interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, element_filter);
 
   interpolateElementalFieldFromIntegrationPoints(
       field, interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, result, ghost_type, element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     interpolateElementalFieldFromIntegrationPoints(
         const ElementTypeMapArray<Real> & field,
         const ElementTypeMapArray<Real> &
             interpolation_points_coordinates_matrices,
         const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
         ElementTypeMapArray<Real> & result, GhostType ghost_type,
         const ElementTypeMapArray<UInt> * element_filter) const {
 
   shape_functions.interpolateElementalFieldFromIntegrationPoints(
       field, interpolation_points_coordinates_matrices,
       quad_points_coordinates_inv_matrices, result, ghost_type, element_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct InterpolateHelper {
       template <class S>
       static void call(const S & /*unused*/, const Vector<Real> & /*unused*/,
                        UInt /*unused*/, const Matrix<Real> & /*unused*/,
-                       Vector<Real> & /*unused*/,
-                       ElementType /*unused*/,
+                       Vector<Real> & /*unused*/, ElementType /*unused*/,
                        GhostType /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define INTERPOLATE(type)                                                      \
   shape_functions.template interpolate<type>(                                  \
       real_coords, element, nodal_values, interpolated, ghost_type);
 
 #define AKANTU_SPECIALIZE_INTERPOLATE_HELPER(kind)                             \
   template <> struct InterpolateHelper<kind> {                                 \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, UInt element,           \
                      const Matrix<Real> & nodal_values,                        \
-                     Vector<Real> & interpolated, ElementType type,    \
-                     GhostType ghost_type) {                           \
+                     Vector<Real> & interpolated, ElementType type,            \
+                     GhostType ghost_type) {                                   \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind);                     \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INTERPOLATE_HELPER,
                                AKANTU_FE_ENGINE_LIST_INTERPOLATE)
 
 #undef AKANTU_SPECIALIZE_INTERPOLATE_HELPER
 #undef INTERPOLATE
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::interpolate(
     const Vector<Real> & real_coords, const Matrix<Real> & nodal_values,
     Vector<Real> & interpolated, const Element & element) const {
 
   AKANTU_DEBUG_IN();
 
   fe_engine::details::InterpolateHelper<kind>::call(
       shape_functions, real_coords, element.element, nodal_values, interpolated,
       element.type, element.ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   computeNormalsOnIntegrationPoints(mesh.getNodes(), ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                       GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   //  Real * coord = mesh.getNodes().storage();
   UInt spatial_dimension = mesh.getSpatialDimension();
 
   // allocate the normal arrays
   normals_on_integration_points.initialize(
       *this, _nb_component = spatial_dimension,
       _spatial_dimension = element_dimension, _ghost_type = ghost_type,
       _element_kind = kind);
 
   // loop over the type to build the normals
   for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) {
     auto & normals_on_quad = normals_on_integration_points(type, ghost_type);
     computeNormalsOnIntegrationPoints(field, normals_on_quad, type, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeNormalsOnIntegrationPoints {
       template <template <ElementKind, class> class I,
                 template <ElementKind> class S, ElementKind k, class IOF>
       static void call(const FEEngineTemplate<I, S, k, IOF> & /*unused*/,
                        const Array<Real> & /*unused*/, Array<Real> & /*unused*/,
-                       ElementType /*unused*/,
-                       GhostType /*unused*/) {
+                       ElementType /*unused*/, GhostType /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_NORMALS_ON_INTEGRATION_POINTS(type)                            \
   fem.template computeNormalsOnIntegrationPoints<type>(field, normal,          \
                                                        ghost_type);
 
 #define AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS(kind)          \
   template <> struct ComputeNormalsOnIntegrationPoints<kind> {                 \
     template <template <ElementKind, class> class I,                           \
               template <ElementKind> class S, ElementKind k, class IOF>        \
     static void call(const FEEngineTemplate<I, S, k, IOF> & fem,               \
                      const Array<Real> & field, Array<Real> & normal,          \
-                     ElementType type, GhostType ghost_type) { \
+                     ElementType type, GhostType ghost_type) {                 \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_INTEGRATION_POINTS,  \
                                        kind);                                  \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS,
         AKANTU_FE_ENGINE_LIST_COMPUTE_NORMALS_ON_INTEGRATION_POINTS)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS
 #undef COMPUTE_NORMALS_ON_INTEGRATION_POINTS
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const Array<Real> & field,
-                                      Array<Real> & normal,
-                                      ElementType type,
+                                      Array<Real> & normal, ElementType type,
                                       GhostType ghost_type) const {
   fe_engine::details::ComputeNormalsOnIntegrationPoints<kind>::call(
       *this, field, normal, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 template <ElementType type>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
     computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                       Array<Real> & normal,
                                       GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   if (type == _point_1) {
     computeNormalsOnIntegrationPointsPoint1(field, normal, ghost_type);
     return;
   }
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   UInt nb_points = getNbIntegrationPoints(type, ghost_type);
 
   UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
   normal.resize(nb_element * nb_points);
   Array<Real>::matrix_iterator normals_on_quad =
       normal.begin_reinterpret(spatial_dimension, nb_points, nb_element);
   Array<Real> f_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, field, f_el, type, ghost_type);
 
   const Matrix<Real> & quads =
       integrator.template getIntegrationPoints<type>(ghost_type);
 
   Array<Real>::matrix_iterator f_it =
       f_el.begin(spatial_dimension, nb_nodes_per_element);
 
   for (UInt elem = 0; elem < nb_element; ++elem) {
     ElementClass<type>::computeNormalsOnNaturalCoordinates(quads, *f_it,
                                                            *normals_on_quad);
     ++normals_on_quad;
     ++f_it;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 template <ElementKind kind> struct InverseMapHelper {
   template <class S>
-  static void
-  call(const S & /*shape_functions*/, const Vector<Real> & /*real_coords*/,
-       UInt /*element*/, ElementType /*type*/,
-       Vector<Real> & /*natural_coords*/, GhostType /*ghost_type*/) {
+  static void call(const S & /*shape_functions*/,
+                   const Vector<Real> & /*real_coords*/, UInt /*element*/,
+                   ElementType /*type*/, Vector<Real> & /*natural_coords*/,
+                   GhostType /*ghost_type*/) {
     AKANTU_TO_IMPLEMENT();
   }
 };
 
 #define INVERSE_MAP(type)                                                      \
   shape_functions.template inverseMap<type>(real_coords, element,              \
                                             natural_coords, ghost_type);
 
 #define AKANTU_SPECIALIZE_INVERSE_MAP_HELPER(kind)                             \
   template <> struct InverseMapHelper<kind> {                                  \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, UInt element,           \
-                     ElementType type, Vector<Real> & natural_coords,  \
-                     GhostType ghost_type) {                           \
+                     ElementType type, Vector<Real> & natural_coords,          \
+                     GhostType ghost_type) {                                   \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(INVERSE_MAP, kind);                     \
     }                                                                          \
   };
 
 AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INVERSE_MAP_HELPER,
                            AKANTU_FE_ENGINE_LIST_INVERSE_MAP)
 
 #undef AKANTU_SPECIALIZE_INVERSE_MAP_HELPER
 #undef INVERSE_MAP
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::inverseMap(
     const Vector<Real> & real_coords, UInt element, ElementType type,
     Vector<Real> & natural_coords, GhostType ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   InverseMapHelper<kind>::call(shape_functions, real_coords, element, type,
                                natural_coords, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ContainsHelper {
       template <class S>
       static void call(const S & /*unused*/, const Vector<Real> & /*unused*/,
                        UInt /*unused*/, ElementType /*unused*/,
                        GhostType /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define CONTAINS(type)                                                         \
   contain = shape_functions.template contains<type>(real_coords, element,      \
                                                     ghost_type);
 
 #define AKANTU_SPECIALIZE_CONTAINS_HELPER(kind)                                \
   template <> struct ContainsHelper<kind> {                                    \
     template <template <ElementKind> class S, ElementKind k>                   \
     static bool call(const S<k> & shape_functions,                             \
                      const Vector<Real> & real_coords, UInt element,           \
-                     ElementType type, GhostType ghost_type) { \
+                     ElementType type, GhostType ghost_type) {                 \
       bool contain = false;                                                    \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(CONTAINS, kind);                        \
       return contain;                                                          \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_CONTAINS_HELPER,
                                AKANTU_FE_ENGINE_LIST_CONTAINS)
 
 #undef AKANTU_SPECIALIZE_CONTAINS_HELPER
 #undef CONTAINS
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline bool FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::contains(
     const Vector<Real> & real_coords, UInt element, ElementType type,
     GhostType ghost_type) const {
   return fe_engine::details::ContainsHelper<kind>::call(
       shape_functions, real_coords, element, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeShapesHelper {
       template <class S>
       static void call(const S & /*unused*/, const Vector<Real> & /*unused*/,
                        UInt /*unused*/, const ElementType /*unused*/,
-                       Vector<Real> & /*unused*/,
-                       GhostType /*unused*/) {
+                       Vector<Real> & /*unused*/, GhostType /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_SHAPES(type)                                                   \
   shape_functions.template computeShapes<type>(real_coords, element, shapes,   \
                                                ghost_type);
 
 #define AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER(kind)                          \
   template <> struct ComputeShapesHelper<kind> {                               \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, UInt element,           \
                      const ElementType type, Vector<Real> & shapes,            \
-                     GhostType ghost_type) {                           \
+                     GhostType ghost_type) {                                   \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPES, kind);                  \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER,
                                AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER
 #undef COMPUTE_SHAPES
   } // namespace details
 } // namespace fe_engine
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapes(
     const Vector<Real> & real_coords, UInt element, ElementType type,
     Vector<Real> & shapes, GhostType ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   fe_engine::details::ComputeShapesHelper<kind>::call(
       shape_functions, real_coords, element, type, shapes, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct ComputeShapeDerivativesHelper {
       template <class S>
       static void call(__attribute__((unused)) const S & shape_functions,
                        __attribute__((unused)) const Vector<Real> & real_coords,
                        __attribute__((unused)) UInt element,
                        __attribute__((unused)) const ElementType type,
                        __attribute__((unused)) Matrix<Real> & shape_derivatives,
                        __attribute__((unused)) GhostType ghost_type) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define COMPUTE_SHAPE_DERIVATIVES(type)                                        \
   Matrix<Real> coords_mat(real_coords.storage(), shape_derivatives.rows(), 1); \
   Tensor3<Real> shapesd_tensor(shape_derivatives.storage(),                    \
                                shape_derivatives.rows(),                       \
                                shape_derivatives.cols(), 1);                   \
   shape_functions.template computeShapeDerivatives<type>(                      \
       coords_mat, element, shapesd_tensor, ghost_type);
 
 #define AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER(kind)               \
   template <> struct ComputeShapeDerivativesHelper<kind> {                     \
     template <class S>                                                         \
     static void call(const S & shape_functions,                                \
                      const Vector<Real> & real_coords, UInt element,           \
                      const ElementType type, Matrix<Real> & shape_derivatives, \
-                     GhostType ghost_type) {                           \
+                     GhostType ghost_type) {                                   \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPE_DERIVATIVES, kind);       \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(
         AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER,
         AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES_DERIVATIVES)
 
 #undef AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER
 #undef COMPUTE_SHAPE_DERIVATIVES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapeDerivatives(
     const Vector<Real> & real_coords, UInt element, ElementType type,
     Matrix<Real> & shape_derivatives, GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   fe_engine::details::ComputeShapeDerivativesHelper<kind>::call(
       shape_functions, real_coords, element, type, shape_derivatives,
       ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetNbIntegrationPointsHelper {};
 
 #define GET_NB_INTEGRATION_POINTS(type)                                        \
   nb_quad_points = integrator.template getNbIntegrationPoints<type>(ghost_type);
 
 #define AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER(kind)               \
   template <> struct GetNbIntegrationPointsHelper<kind> {                      \
     template <template <ElementKind, class> class I, ElementKind k, class IOF> \
     static UInt call(const I<k, IOF> & integrator, const ElementType type,     \
-                     GhostType ghost_type) {                           \
+                     GhostType ghost_type) {                                   \
       UInt nb_quad_points = 0;                                                 \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_NB_INTEGRATION_POINTS, kind);       \
       return nb_quad_points;                                                   \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER
 #undef GET_NB_INTEGRATION
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline UInt
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getNbIntegrationPoints(
     ElementType type, GhostType ghost_type) const {
   return fe_engine::details::GetNbIntegrationPointsHelper<kind>::call(
       integrator, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetShapesHelper {};
 
 #define GET_SHAPES(type) ret = &(shape_functions.getShapes(type, ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_SHAPES_HELPER(kind)                              \
   template <> struct GetShapesHelper<kind> {                                   \
     template <class S>                                                         \
     static const Array<Real> & call(const S & shape_functions,                 \
                                     const ElementType type,                    \
-                                    GhostType ghost_type) {            \
+                                    GhostType ghost_type) {                    \
       const Array<Real> * ret = NULL;                                          \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES, kind);                      \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_SHAPES_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_SHAPES_HELPER
 #undef GET_SHAPES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Array<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getShapes(
     ElementType type, GhostType ghost_type,
     __attribute__((unused)) UInt id) const {
   return fe_engine::details::GetShapesHelper<kind>::call(shape_functions, type,
                                                          ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetShapesDerivativesHelper {
       template <template <ElementKind> class S, ElementKind k>
-      static const Array<Real> &
-      call(const S<k> & /*unused*/, ElementType /*unused*/,
-           GhostType /*unused*/, UInt /*unused*/) {
+      static const Array<Real> & call(const S<k> & /*unused*/,
+                                      ElementType /*unused*/,
+                                      GhostType /*unused*/, UInt /*unused*/) {
         AKANTU_TO_IMPLEMENT();
       }
     };
 
 #define GET_SHAPES_DERIVATIVES(type)                                           \
   ret = &(shape_functions.getShapesDerivatives(type, ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER(kind)                  \
   template <> struct GetShapesDerivativesHelper<kind> {                        \
     template <template <ElementKind> class S, ElementKind k>                   \
     static const Array<Real> &                                                 \
     call(const S<k> & shape_functions, const ElementType type,                 \
-         GhostType ghost_type, __attribute__((unused)) UInt id) {      \
+         GhostType ghost_type, __attribute__((unused)) UInt id) {              \
       const Array<Real> * ret = NULL;                                          \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES_DERIVATIVES, kind);          \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER,
                                AKANTU_FE_ENGINE_LIST_GET_SHAPES_DERIVATIVES)
 
 #undef AKANTU_SPECIALIZE_GET_SHAPE_DERIVATIVES_HELPER
 #undef GET_SHAPES_DERIVATIVES
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Array<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getShapesDerivatives(
     ElementType type, GhostType ghost_type,
     __attribute__((unused)) UInt id) const {
   return fe_engine::details::GetShapesDerivativesHelper<kind>::call(
       shape_functions, type, ghost_type, id);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Helper class to be able to write a partial specialization on the element kind
  */
 namespace fe_engine {
   namespace details {
     template <ElementKind kind> struct GetIntegrationPointsHelper {};
 
 #define GET_INTEGRATION_POINTS(type)                                           \
   ret = &(integrator.template getIntegrationPoints<type>(ghost_type));
 
 #define AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER(kind)                  \
   template <> struct GetIntegrationPointsHelper<kind> {                        \
     template <template <ElementKind, class> class I, ElementKind k, class IOF> \
     static const Matrix<Real> & call(const I<k, IOF> & integrator,             \
                                      const ElementType type,                   \
-                                     GhostType ghost_type) {           \
+                                     GhostType ghost_type) {                   \
       const Matrix<Real> * ret = NULL;                                         \
       AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_INTEGRATION_POINTS, kind);          \
       return *ret;                                                             \
     }                                                                          \
   };
 
     AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER)
 
 #undef AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER
 #undef GET_INTEGRATION_POINTS
   } // namespace details
 } // namespace fe_engine
 
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline const Matrix<Real> &
 FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getIntegrationPoints(
     ElementType type, GhostType ghost_type) const {
   return fe_engine::details::GetIntegrationPointsHelper<kind>::call(
       integrator, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::printself(
     std::ostream & stream, int indent) const {
   std::string space(indent, AKANTU_INDENT);
 
   stream << space << "FEEngineTemplate [" << std::endl;
   stream << space << " + parent [" << std::endl;
   FEEngine::printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << " + shape functions [" << std::endl;
   shape_functions.printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << " + integrator [" << std::endl;
   integrator.printself(stream, indent + 3);
   stream << space << "   ]" << std::endl;
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsAdded(
     const Array<Element> & new_elements, const NewElementsEvent & /*unused*/) {
   integrator.onElementsAdded(new_elements);
   shape_functions.onElementsAdded(new_elements);
 }
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsRemoved(
     const Array<Element> & /*unused*/,
     const ElementTypeMapArray<UInt> & /*unused*/,
     const RemovedElementsEvent & /*unused*/) {}
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsChanged(
     const Array<Element> & /*unused*/, const Array<Element> & /*unused*/,
     const ElementTypeMapArray<UInt> & /*unused*/,
     const ChangedElementsEvent & /*unused*/) {}
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind, class IntegrationOrderFunctor>
 inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::
-    computeNormalsOnIntegrationPointsPoint1(
-        const Array<Real> & /*unused*/, Array<Real> & normal,
-        GhostType ghost_type) const {
+    computeNormalsOnIntegrationPointsPoint1(const Array<Real> & /*unused*/,
+                                            Array<Real> & normal,
+                                            GhostType ghost_type) const {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == 1,
                       "Mesh dimension must be 1 to compute normals on points!");
   const auto type = _point_1;
   auto spatial_dimension = mesh.getSpatialDimension();
   // UInt nb_nodes_per_element  = Mesh::getNbNodesPerElement(type);
   auto nb_points = getNbIntegrationPoints(type, ghost_type);
   const auto & connectivity = mesh.getConnectivity(type, ghost_type);
   auto nb_element = connectivity.size();
 
   normal.resize(nb_element * nb_points);
   auto normals_on_quad =
       normal.begin_reinterpret(spatial_dimension, nb_points, nb_element);
   const auto & segments = mesh.getElementToSubelement(type, ghost_type);
   const auto & coords = mesh.getNodes();
 
   const Mesh * mesh_segment;
   if (mesh.isMeshFacets()) {
     mesh_segment = &(mesh.getMeshParent());
   } else {
     mesh_segment = &mesh;
   }
 
   for (UInt elem = 0; elem < nb_element; ++elem) {
     UInt nb_segment = segments(elem).size();
     AKANTU_DEBUG_ASSERT(
         nb_segment > 0,
         "Impossible to compute a normal on a point connected to 0 segments");
 
     Real normal_value = 1;
     if (nb_segment == 1) {
       auto point = connectivity(elem);
       const auto segment = segments(elem)[0];
       const auto & segment_connectivity =
           mesh_segment->getConnectivity(segment.type, segment.ghost_type);
       Vector<UInt> segment_points = segment_connectivity.begin(
           Mesh::getNbNodesPerElement(segment.type))[segment.element];
       Real difference;
       if (segment_points(0) == point) {
         difference = coords(elem) - coords(segment_points(1));
       } else {
         difference = coords(elem) - coords(segment_points(0));
       }
 
       normal_value = difference / std::abs(difference);
     }
 
     for (UInt n(0); n < nb_points; ++n) {
       (*normals_on_quad)(0, n) = normal_value;
     }
     ++normals_on_quad;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
diff --git a/src/fe_engine/shape_cohesive.hh b/src/fe_engine/shape_cohesive.hh
index accd39f41..21818137c 100644
--- a/src/fe_engine/shape_cohesive.hh
+++ b/src/fe_engine/shape_cohesive.hh
@@ -1,190 +1,188 @@
 /**
  * @file   shape_cohesive.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Tue Feb 15 2011
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  shape functions for cohesive elements description
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 #include "shape_lagrange.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SHAPE_COHESIVE_HH_
 #define AKANTU_SHAPE_COHESIVE_HH_
 
 namespace akantu {
 
 struct CohesiveReduceFunctionMean {
   inline Real operator()(Real u_plus, Real u_minus) {
     return .5 * (u_plus + u_minus);
   }
 };
 
 struct CohesiveReduceFunctionOpening {
   inline Real operator()(Real u_plus, Real u_minus) {
     return (u_plus - u_minus);
   }
 };
 
 template <> class ShapeLagrange<_ek_cohesive> : public ShapeLagrangeBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   ShapeLagrange(const Mesh & mesh, UInt spatial_dimension,
                 const ID & id = "shape_cohesive");
 
   ~ShapeLagrange() override = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   inline void initShapeFunctions(const Array<Real> & nodes,
                                  const Matrix<Real> & integration_points,
-                                 ElementType type,
-                                 GhostType ghost_type);
+                                 ElementType type, GhostType ghost_type);
 
   /// extract the nodal values and store them per element
   template <ElementType type, class ReduceFunction>
   void extractNodalToElementField(
       const Array<Real> & nodal_f, Array<Real> & elemental_f,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// computes the shape functions derivatives for given interpolation points
   template <ElementType type>
   void computeShapeDerivativesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
       Array<Real> & shape_derivatives, GhostType ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   void computeShapeDerivativesOnIntegrationPoints(
       const Array<Real> & nodes, const Matrix<Real> & integration_points,
-      Array<Real> & shape_derivatives, ElementType type,
-      GhostType ghost_type,
+      Array<Real> & shape_derivatives, ElementType type, GhostType ghost_type,
       const Array<UInt> & filter_elements) const override;
 
   /// pre compute all shapes on the element integration points from natural
   /// coordinates
   template <ElementType type>
   void precomputeShapesOnIntegrationPoints(const Array<Real> & nodes,
                                            GhostType ghost_type);
 
   /// pre compute all shape derivatives on the element integration points from
   /// natural coordinates
   template <ElementType type>
   void precomputeShapeDerivativesOnIntegrationPoints(const Array<Real> & nodes,
                                                      GhostType ghost_type);
 
   /// interpolate nodal values on the integration points
   template <ElementType type, class ReduceFunction>
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   template <ElementType type>
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const {
     interpolateOnIntegrationPoints<type, CohesiveReduceFunctionMean>(
         u, uq, nb_degree_of_freedom, ghost_type, filter_elements);
   }
 
   /// compute the gradient of u on the integration points in the natural
   /// coordinates
   template <ElementType type>
   void gradientOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const {
     variationOnIntegrationPoints<type, CohesiveReduceFunctionMean>(
         u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements);
   }
 
   /* ------------------------------------------------------------------------ */
   template <ElementType type>
   void computeBtD(const Array<Real> & /*Ds*/, Array<Real> & /*BtDs*/,
                   GhostType /*ghost_type*/,
                   const Array<UInt> & /*filter_elements*/) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   template <ElementType type>
   void computeBtDB(const Array<Real> & /*Ds*/, Array<Real> & /*BtDBs*/,
-                   UInt /*order_d*/, GhostType/*ghost_type*/,
+                   UInt /*order_d*/, GhostType /*ghost_type*/,
                    const Array<UInt> & /*filter_elements*/) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// multiply a field by shape functions
   template <ElementType type>
   void
   computeNtb(const Array<Real> & /*bs*/, Array<Real> & /*Ntbs*/,
              GhostType /*ghost_type*/,
              const Array<UInt> & /*filter_elements*/ = empty_filter) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   template <ElementType type>
   void computeNtbN(const Array<Real> & /*bs*/, Array<Real> & /*NtbNs*/,
                    GhostType /*ghost_type*/,
                    const Array<UInt> & /*filter_elements*/) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   /* ------------------------------------------------------------------------ */
   /// compute the gradient of u on the integration points
   template <ElementType type, class ReduceFunction>
   void variationOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 
   /// compute the normals to the field u on integration points
   template <ElementType type, class ReduceFunction>
   void computeNormalsOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & normals_u,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const;
 };
 
 /// standard output stream operator
 template <class ShapeFunction>
 inline std::ostream & operator<<(std::ostream & stream,
                                  const ShapeCohesive<ShapeFunction> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "shape_cohesive_inline_impl.hh"
 
 #endif /* AKANTU_SHAPE_COHESIVE_HH_ */
diff --git a/src/fe_engine/shape_lagrange_inline_impl.hh b/src/fe_engine/shape_lagrange_inline_impl.hh
index 54133b3db..0cd7ee9c4 100644
--- a/src/fe_engine/shape_lagrange_inline_impl.hh
+++ b/src/fe_engine/shape_lagrange_inline_impl.hh
@@ -1,579 +1,579 @@
 /**
  * @file   shape_lagrange_inline_impl.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Oct 27 2010
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  ShapeLagrange inline implementation
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_iterators.hh"
 #include "aka_voigthelper.hh"
 #include "fe_engine.hh"
 #include "shape_lagrange.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_HH_
 #define AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_HH_
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 #define INIT_SHAPE_FUNCTIONS(type)                                             \
   setIntegrationPointsByType<type>(integration_points, ghost_type);            \
   precomputeShapesOnIntegrationPoints<type>(nodes, ghost_type);                \
   if (ElementClass<type>::getNaturalSpaceDimension() ==                        \
           mesh.getSpatialDimension() ||                                        \
       kind != _ek_regular)                                                     \
     precomputeShapeDerivativesOnIntegrationPoints<type>(nodes, ghost_type);
 
 template <ElementKind kind>
 inline void ShapeLagrange<kind>::initShapeFunctions(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     ElementType type, GhostType ghost_type) {
   AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
 }
 
 #undef INIT_SHAPE_FUNCTIONS
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 inline void ShapeLagrange<kind>::computeShapeDerivativesOnCPointsByElement(
     const Matrix<Real> & node_coords, const Matrix<Real> & natural_coords,
     Tensor3<Real> & shapesd) const {
   AKANTU_DEBUG_IN();
 
   // compute dnds
   Tensor3<Real> dnds(node_coords.rows(), node_coords.cols(),
                      natural_coords.cols());
   ElementClass<type>::computeDNDS(natural_coords, dnds);
   // compute jacobian
   Tensor3<Real> J(node_coords.rows(), natural_coords.rows(),
                   natural_coords.cols());
   ElementClass<type>::computeJMat(dnds, node_coords, J);
 
   // compute dndx
   ElementClass<type>::computeShapeDerivatives(J, dnds, shapesd);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::inverseMap(const Vector<Real> & real_coords,
                                      UInt elem, Vector<Real> & natural_coords,
                                      GhostType ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage();
   Matrix<Real> nodes_coord(spatial_dimension, nb_nodes_per_element);
 
   mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(),
                                      elem_val + elem * nb_nodes_per_element,
                                      nb_nodes_per_element, spatial_dimension);
 
   ElementClass<type>::inverseMap(real_coords, nodes_coord, natural_coords);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 bool ShapeLagrange<kind>::contains(const Vector<Real> & real_coords, UInt elem,
                                    GhostType ghost_type) const {
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   Vector<Real> natural_coords(spatial_dimension);
 
   inverseMap<type>(real_coords, elem, natural_coords, ghost_type);
   return ElementClass<type>::contains(natural_coords);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::interpolate(const Vector<Real> & real_coords,
                                       UInt elem,
                                       const Matrix<Real> & nodal_values,
                                       Vector<Real> & interpolated,
                                       GhostType ghost_type) const {
   UInt nb_shapes = ElementClass<type>::getShapeSize();
   Vector<Real> shapes(nb_shapes);
   computeShapes<type>(real_coords, elem, shapes, ghost_type);
   ElementClass<type>::interpolate(nodal_values, shapes, interpolated);
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeShapes(const Vector<Real> & real_coords,
                                         UInt elem, Vector<Real> & shapes,
                                         GhostType ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   Vector<Real> natural_coords(spatial_dimension);
 
   inverseMap<type>(real_coords, elem, natural_coords, ghost_type);
   ElementClass<type>::computeShapes(natural_coords, shapes);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeShapeDerivatives(
     const Matrix<Real> & real_coords, UInt elem, Tensor3<Real> & shapesd,
     GhostType ghost_type) const {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   UInt nb_points = real_coords.cols();
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == shapesd.size(0) &&
                           nb_nodes_per_element == shapesd.size(1),
                       "Shape size doesn't match");
   AKANTU_DEBUG_ASSERT(nb_points == shapesd.size(2),
                       "Number of points doesn't match shapes size");
 
   Matrix<Real> natural_coords(spatial_dimension, nb_points);
 
   // Creates the matrix of natural coordinates
   for (UInt i = 0; i < nb_points; i++) {
     Vector<Real> real_point = real_coords(i);
     Vector<Real> natural_point = natural_coords(i);
 
     inverseMap<type>(real_point, elem, natural_point, ghost_type);
   }
 
   UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage();
   Matrix<Real> nodes_coord(spatial_dimension, nb_nodes_per_element);
 
   mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(),
                                      elem_val + elem * nb_nodes_per_element,
                                      nb_nodes_per_element, spatial_dimension);
 
   computeShapeDerivativesOnCPointsByElement<type>(nodes_coord, natural_coords,
                                                   shapesd);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 ShapeLagrange<kind>::ShapeLagrange(const Mesh & mesh, UInt spatial_dimension,
                                    const ID & id)
     : ShapeLagrangeBase(mesh, spatial_dimension, kind, id) {}
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
     Array<Real> & shape_derivatives, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
 
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   UInt nb_points = integration_points.cols();
   UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
 
   UInt size_of_shapesd = ElementClass<type>::getShapeDerivativesSize();
   AKANTU_DEBUG_ASSERT(shape_derivatives.getNbComponent() == size_of_shapesd,
                       "The shapes_derivatives array does not have the correct "
                           << "number of component");
   shape_derivatives.resize(nb_element * nb_points);
 
   Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type,
                                        filter_elements);
 
   Real * shapesd_val = shape_derivatives.storage();
   Array<Real>::matrix_iterator x_it =
       x_el.begin(spatial_dimension, nb_nodes_per_element);
 
   if (filter_elements != empty_filter) {
     nb_element = filter_elements.size();
   }
 
   for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) {
     if (filter_elements != empty_filter) {
       shapesd_val = shape_derivatives.storage() +
                     filter_elements(elem) * size_of_shapesd * nb_points;
     }
 
     Matrix<Real> & X = *x_it;
     Tensor3<Real> B(shapesd_val, spatial_dimension, nb_nodes_per_element,
                     nb_points);
     computeShapeDerivativesOnCPointsByElement<type>(X, integration_points, B);
 
     if (filter_elements == empty_filter) {
       shapesd_val += size_of_shapesd * nb_points;
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 void ShapeLagrange<kind>::computeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, const Matrix<Real> & integration_points,
-    Array<Real> & shape_derivatives, ElementType type,
-    GhostType ghost_type, const Array<UInt> & filter_elements) const {
+    Array<Real> & shape_derivatives, ElementType type, GhostType ghost_type,
+    const Array<UInt> & filter_elements) const {
 #define AKANTU_COMPUTE_SHAPES(type)                                            \
   computeShapeDerivativesOnIntegrationPoints<type>(                            \
       nodes, integration_points, shape_derivatives, ghost_type,                \
       filter_elements);
 
   AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES);
 
 #undef AKANTU_COMPUTE_SHAPES
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::precomputeShapesOnIntegrationPoints(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   Matrix<Real> & natural_coords = integration_points(type, ghost_type);
   UInt size_of_shapes = ElementClass<type>::getShapeSize();
 
   Array<Real> & shapes_tmp =
       shapes.alloc(0, size_of_shapes, itp_type, ghost_type);
 
   this->computeShapesOnIntegrationPoints<type>(nodes, natural_coords,
                                                shapes_tmp, ghost_type);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::precomputeShapeDerivativesOnIntegrationPoints(
     const Array<Real> & nodes, GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   Matrix<Real> & natural_coords = integration_points(type, ghost_type);
   UInt size_of_shapesd = ElementClass<type>::getShapeDerivativesSize();
 
   Array<Real> & shapes_derivatives_tmp =
       shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type);
 
   this->computeShapeDerivativesOnIntegrationPoints<type>(
       nodes, natural_coords, shapes_derivatives_tmp, ghost_type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::interpolateOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
     const Array<Real> & shapes, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   Array<Real> u_el(0, nb_degree_of_freedom * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
                                        filter_elements);
 
   this->interpolateElementalFieldOnIntegrationPoints<type>(
       u_el, out_uq, ghost_type, shapes, filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::interpolateOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_degree_of_freedom,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   AKANTU_DEBUG_ASSERT(shapes.exists(itp_type, ghost_type),
                       "No shapes for the type "
                           << shapes.printType(itp_type, ghost_type));
 
   this->interpolateOnIntegrationPoints<type>(in_u, out_uq, nb_degree_of_freedom,
                                              shapes(itp_type, ghost_type),
                                              ghost_type, filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::gradientOnIntegrationPoints(
     const Array<Real> & in_u, Array<Real> & out_nablauq,
     UInt nb_degree_of_freedom, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
   AKANTU_DEBUG_ASSERT(
       shapes_derivatives.exists(itp_type, ghost_type),
       "No shapes derivatives for the type "
           << shapes_derivatives.printType(itp_type, ghost_type));
 
   UInt nb_nodes_per_element =
       ElementClass<type>::getNbNodesPerInterpolationElement();
 
   Array<Real> u_el(0, nb_degree_of_freedom * nb_nodes_per_element);
   FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
                                        filter_elements);
 
   this->gradientElementalFieldOnIntegrationPoints<type>(
       u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type),
       filter_elements);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeBtD(
     const Array<Real> & Ds, Array<Real> & BtDs, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   auto itp_type = ElementClassProperty<type>::interpolation_type;
   const auto & shapes_derivatives =
       this->shapes_derivatives(itp_type, ghost_type);
 
   auto spatial_dimension = mesh.getSpatialDimension();
   auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
 
   Array<Real> shapes_derivatives_filtered(0,
                                           shapes_derivatives.getNbComponent());
   auto && view =
       make_view(shapes_derivatives, spatial_dimension, nb_nodes_per_element);
   auto B_it = view.begin();
   auto B_end = view.end();
 
   if (filter_elements != empty_filter) {
     FEEngine::filterElementalData(this->mesh, shapes_derivatives,
                                   shapes_derivatives_filtered, type, ghost_type,
                                   filter_elements);
     auto && view = make_view(shapes_derivatives_filtered, spatial_dimension,
                              nb_nodes_per_element);
     B_it = view.begin();
     B_end = view.end();
   }
 
   for (auto && values :
        zip(range(B_it, B_end),
            make_view(Ds, Ds.getNbComponent() / spatial_dimension,
                      spatial_dimension),
            make_view(BtDs, BtDs.getNbComponent() / nb_nodes_per_element,
                      nb_nodes_per_element))) {
     const auto & B = std::get<0>(values);
     const auto & D = std::get<1>(values);
     auto & Bt_D = std::get<2>(values);
     // transposed due to the storage layout of B
     Bt_D.template mul<false, false>(D, B);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeBtDB(
     const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
     GhostType ghost_type, const Array<UInt> & filter_elements) const {
   auto itp_type = ElementClassProperty<type>::interpolation_type;
   const auto & shapes_derivatives =
       this->shapes_derivatives(itp_type, ghost_type);
 
   constexpr auto dim = ElementClass<type>::getSpatialDimension();
   auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
 
   Array<Real> shapes_derivatives_filtered(0,
                                           shapes_derivatives.getNbComponent());
   auto && view = make_view(shapes_derivatives, dim, nb_nodes_per_element);
   auto B_it = view.begin();
   auto B_end = view.end();
 
   if (filter_elements != empty_filter) {
     FEEngine::filterElementalData(this->mesh, shapes_derivatives,
                                   shapes_derivatives_filtered, type, ghost_type,
                                   filter_elements);
     auto && view =
         make_view(shapes_derivatives_filtered, dim, nb_nodes_per_element);
     B_it = view.begin();
     B_end = view.end();
   }
 
   if (order_d == 4) {
     UInt tangent_size = VoigtHelper<dim>::size;
     Matrix<Real> B(tangent_size, dim * nb_nodes_per_element);
     Matrix<Real> Bt_D(dim * nb_nodes_per_element, tangent_size);
 
     for (auto && values :
          zip(range(B_it, B_end), make_view(Ds, tangent_size, tangent_size),
              make_view(BtDBs, dim * nb_nodes_per_element,
                        dim * nb_nodes_per_element))) {
       const auto & Bfull = std::get<0>(values);
       const auto & D = std::get<1>(values);
       auto & Bt_D_B = std::get<2>(values);
 
       VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(Bfull, B,
                                                          nb_nodes_per_element);
       Bt_D.template mul<true, false>(B, D);
       Bt_D_B.template mul<false, false>(Bt_D, B);
     }
   } else if (order_d == 2) {
     Matrix<Real> Bt_D(nb_nodes_per_element, dim);
     for (auto && values :
          zip(range(B_it, B_end), make_view(Ds, dim, dim),
              make_view(BtDBs, nb_nodes_per_element, nb_nodes_per_element))) {
       const auto & B = std::get<0>(values);
       const auto & D = std::get<1>(values);
       auto & Bt_D_B = std::get<2>(values);
       Bt_D.template mul<true, false>(B, D);
       Bt_D_B.template mul<false, false>(Bt_D, B);
     }
   }
 }
 
 template <>
 template <>
 inline void ShapeLagrange<_ek_regular>::computeBtDB<_point_1>(
     const Array<Real> & /*Ds*/, Array<Real> & /*BtDBs*/, UInt /*order_d*/,
     GhostType /*ghost_type*/, const Array<UInt> & /*filter_elements*/) const {
   AKANTU_TO_IMPLEMENT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeNtbN(
-    const Array<Real> & bs, Array<Real> & NtbNs, UInt order_d,
-    GhostType ghost_type, const Array<UInt> & filter_elements) const {
+    const Array<Real> & bs, Array<Real> & NtbNs, GhostType ghost_type,
+    const Array<UInt> & filter_elements) const {
 
   auto itp_type = ElementClassProperty<type>::interpolation_type;
-  UInt size_of_shapes = ElementClass<type>::getShapeSize();
-  UInt nb_degree_of_freedom = bs.getNbComponent();
+  auto size_of_shapes = ElementClass<type>::getShapeSize();
+  auto nb_degree_of_freedom = bs.getNbComponent();
 
   auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);
   Array<Real> shapes_filtered(0, size_of_shapes);
 
   auto && view = make_view(shapes(itp_type, ghost_type), 1, size_of_shapes);
   auto N_it = view.begin();
   auto N_end = view.end();
 
   if (filter_elements != empty_filter) {
     FEEngine::filterElementalData(this->mesh, shapes(itp_type, ghost_type),
                                   shapes_filtered, type, ghost_type,
                                   filter_elements);
     auto && view = make_view(shapes_filtered, 1, size_of_shapes);
     N_it = view.begin();
     N_end = view.end();
   }
 
   Matrix<Real> Nt_b(nb_nodes_per_element, nb_degree_of_freedom);
   for (auto && values :
        zip(range(N_it, N_end), make_view(bs, nb_degree_of_freedom, 1),
            make_view(NtbNs, nb_nodes_per_element, nb_nodes_per_element))) {
     const auto & N = std::get<0>(values);
     const auto & b = std::get<1>(values);
     auto & Nt_b_N = std::get<2>(values);
     Nt_b.template mul<true, false>(N, b);
     Nt_b_N.template mul<false, false>(Nt_b, N);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <ElementKind kind>
 template <ElementType type>
 void ShapeLagrange<kind>::computeNtb(
     const Array<Real> & bs, Array<Real> & Ntbs, GhostType ghost_type,
     const Array<UInt> & filter_elements) const {
   AKANTU_DEBUG_IN();
 
   Ntbs.resize(bs.size());
 
-  UInt size_of_shapes = ElementClass<type>::getShapeSize();
-  InterpolationType itp_type = ElementClassProperty<type>::interpolation_type;
-  UInt nb_degree_of_freedom = bs.getNbComponent();
+  auto size_of_shapes = ElementClass<type>::getShapeSize();
+  auto itp_type = ElementClassProperty<type>::interpolation_type;
+  auto nb_degree_of_freedom = bs.getNbComponent();
 
   Array<Real> shapes_filtered(0, size_of_shapes);
   auto && view = make_view(shapes(itp_type, ghost_type), 1, size_of_shapes);
   auto N_it = view.begin();
   auto N_end = view.end();
 
   if (filter_elements != empty_filter) {
     FEEngine::filterElementalData(this->mesh, shapes(itp_type, ghost_type),
                                   shapes_filtered, type, ghost_type,
                                   filter_elements);
     auto && view = make_view(shapes_filtered, 1, size_of_shapes);
     N_it = view.begin();
     N_end = view.end();
   }
 
   for (auto && values :
        zip(make_view(bs, nb_degree_of_freedom, 1), range(N_it, N_end),
            make_view(Ntbs, nb_degree_of_freedom, size_of_shapes))) {
     const auto & b = std::get<0>(values);
     const auto & N = std::get<1>(values);
     auto & Ntb = std::get<2>(values);
 
     Ntb.template mul<false, false>(b, N);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
 
 #endif /* AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_HH_ */
diff --git a/src/model/phase_field/phasefield.cc b/src/model/phase_field/phasefield.cc
index ec947afec..8ece183ec 100644
--- a/src/model/phase_field/phasefield.cc
+++ b/src/model/phase_field/phasefield.cc
@@ -1,301 +1,298 @@
 /**
  * @file   pahsefield.cc
  *
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  *
  * @date creation: Mon Mar 2 2020
  * @date last modification: Mon Mar 2 2020
  *
  * @brief  Implementation of the common part of the phasefield 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 "phasefield.hh"
 #include "phase_field_model.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 PhaseField::PhaseField(PhaseFieldModel & model, const ID & id)
-  : Parsable(ParserType::_phasefield, id), id(id),
-    fem(model.getFEEngine()), name(""), model(model),
-    spatial_dimension(this->model.getSpatialDimension()),
-    element_filter("element_filter", id),
-    damage("damage", *this), phi("phi", *this),
-    strain("strain", *this), driving_force("driving_force", *this),
-    damage_energy("damage_energy", *this),
-    damage_energy_density("damage_energy_density", *this) {
+    : Parsable(ParserType::_phasefield, id), id(id), fem(model.getFEEngine()),
+      model(model), spatial_dimension(this->model.getSpatialDimension()),
+      element_filter("element_filter", id), damage("damage", *this),
+      phi("phi", *this), strain("strain", *this),
+      driving_force("driving_force", *this),
+      damage_energy("damage_energy", *this),
+      damage_energy_density("damage_energy_density", *this) {
 
   AKANTU_DEBUG_IN();
 
   /// for each connectivity types allocate the element filer array of the
   /// material
   element_filter.initialize(model.getMesh(),
                             _spatial_dimension = spatial_dimension,
                             _element_kind = _ek_regular);
   this->initialize();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
-  PhaseField::PhaseField(PhaseFieldModel & model, UInt dim, const Mesh & mesh,
-			 FEEngine & fe_engine, const ID & id)
-    : Parsable(ParserType::_phasefield, id), id(id),
-    fem(fe_engine), name(""), model(model),
-    spatial_dimension(this->model.getSpatialDimension()),
-    element_filter("element_filter", id),
-    damage("damage", *this, dim, fe_engine, this->element_filter),
-    phi("phi", *this, dim, fe_engine, this->element_filter),
-    strain("strain", *this, dim, fe_engine, this->element_filter),
-    driving_force("driving_force", *this, dim, fe_engine, this->element_filter),
-    damage_energy("damage_energy", *this, dim, fe_engine, this->element_filter),
-    damage_energy_density("damage_energy_density", *this, dim, fe_engine,
-			  this->element_filter){
+PhaseField::PhaseField(PhaseFieldModel & model, UInt dim, const Mesh & mesh,
+                       FEEngine & fe_engine, const ID & id)
+    : Parsable(ParserType::_phasefield, id), id(id), fem(fe_engine),
+      model(model), spatial_dimension(this->model.getSpatialDimension()),
+      element_filter("element_filter", id),
+      damage("damage", *this, dim, fe_engine, this->element_filter),
+      phi("phi", *this, dim, fe_engine, this->element_filter),
+      strain("strain", *this, dim, fe_engine, this->element_filter),
+      driving_force("driving_force", *this, dim, fe_engine,
+                    this->element_filter),
+      damage_energy("damage_energy", *this, dim, fe_engine,
+                    this->element_filter),
+      damage_energy_density("damage_energy_density", *this, dim, fe_engine,
+                            this->element_filter) {
 
   AKANTU_DEBUG_IN();
 
   /// for each connectivity types allocate the element filer array of the
   /// material
-  element_filter.initialize(mesh,
-                            _spatial_dimension = spatial_dimension,
+  element_filter.initialize(mesh, _spatial_dimension = spatial_dimension,
                             _element_kind = _ek_regular);
   this->initialize();
 
   AKANTU_DEBUG_OUT();
 }
 
-/* -------------------------------------------------------------------------- */  
+/* -------------------------------------------------------------------------- */
 PhaseField::~PhaseField() = default;
 
 /* -------------------------------------------------------------------------- */
 void PhaseField::initialize() {
   registerParam("name", name, std::string(), _pat_parsable | _pat_readable);
   registerParam("l0", l0, Real(0.), _pat_parsable | _pat_readable,
-		"length scale parameter");
+                "length scale parameter");
   registerParam("gc", g_c, _pat_parsable | _pat_readable,
-		"critical local fracture energy density");
+                "critical local fracture energy density");
   registerParam("E", E, _pat_parsable | _pat_readable, "Young's modulus");
   registerParam("nu", nu, _pat_parsable | _pat_readable, "Poisson ratio");
 
   damage.initialize(1);
-  
+
   phi.initialize(1);
   driving_force.initialize(1);
 
   strain.initialize(spatial_dimension * spatial_dimension);
   damage_energy_density.initialize(1);
   damage_energy.initialize(spatial_dimension * spatial_dimension);
 }
 
 /* -------------------------------------------------------------------------- */
 void PhaseField::initPhaseField() {
   AKANTU_DEBUG_IN();
 
   this->phi.initializeHistory();
 
   this->resizeInternals();
 
   updateInternalParameters();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void PhaseField::resizeInternals() {
   AKANTU_DEBUG_IN();
   for (auto it = internal_vectors_real.begin();
-       it != internal_vectors_real.end(); ++it)
+       it != internal_vectors_real.end(); ++it) {
     it->second->resize();
+  }
 
   for (auto it = internal_vectors_uint.begin();
-       it != internal_vectors_uint.end(); ++it)
+       it != internal_vectors_uint.end(); ++it) {
     it->second->resize();
+  }
 
   for (auto it = internal_vectors_bool.begin();
-       it != internal_vectors_bool.end(); ++it)
+       it != internal_vectors_bool.end(); ++it) {
     it->second->resize();
+  }
+
   AKANTU_DEBUG_OUT();
 }
-  
 
 /* -------------------------------------------------------------------------- */
 void PhaseField::updateInternalParameters() {
   this->lambda = this->nu * this->E / ((1 + this->nu) * (1 - 2 * this->nu));
   this->mu = this->E / (2 * (1 + this->nu));
 }
 
 /* -------------------------------------------------------------------------- */
 void PhaseField::computeAllDrivingForces(GhostType ghost_type) {
 
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = model.getSpatialDimension();
 
   for (const auto & type :
        element_filter.elementTypes(spatial_dimension, ghost_type)) {
-    Array<UInt> & elem_filter = element_filter(type, ghost_type);
+    auto & elem_filter = element_filter(type, ghost_type);
 
-    if (elem_filter.size() == 0)
+    if (elem_filter.empty()) {
       continue;
-    
+    }
+
     computeDrivingForce(type, ghost_type);
   }
 
   AKANTU_DEBUG_OUT();
-  
 }
-  
+
 /* -------------------------------------------------------------------------- */
 void PhaseField::assembleInternalForces(GhostType ghost_type) {
 
   AKANTU_DEBUG_IN();
 
   Array<Real> & internal_force = model.getInternalForce();
 
   for (auto type : element_filter.elementTypes(_ghost_type = ghost_type)) {
-  
-    Array<UInt> & elem_filter = element_filter(type, ghost_type);
-    if (elem_filter.size() == 0)
+    auto & elem_filter = element_filter(type, ghost_type);
+    if (elem_filter.empty()) {
       continue;
- 
-    UInt nb_element = elem_filter.size();
-    UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
-    UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type);
+    }
+
+    auto nb_element = elem_filter.size();
+    auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
+    auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type);
 
     Array<Real> nt_driving_force(nb_quadrature_points, nb_nodes_per_element);
     fem.computeNtb(driving_force(type, ghost_type), nt_driving_force, type,
-		   ghost_type, elem_filter);
-    
+                   ghost_type, elem_filter);
+
     Array<Real> int_nt_driving_force(nb_element, nb_nodes_per_element);
-    fem.integrate(nt_driving_force, int_nt_driving_force,
-		  nb_nodes_per_element, type, ghost_type, elem_filter);
+    fem.integrate(nt_driving_force, int_nt_driving_force, nb_nodes_per_element,
+                  type, ghost_type, elem_filter);
 
     model.getDOFManager().assembleElementalArrayLocalArray(
-	 int_nt_driving_force, internal_force, type, ghost_type, 1,
-	 elem_filter);
+        int_nt_driving_force, internal_force, type, ghost_type, 1, elem_filter);
   }
 
   AKANTU_DEBUG_OUT();
 }
- 
+
 /* -------------------------------------------------------------------------- */
 void PhaseField::assembleStiffnessMatrix(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_INFO("Assemble the new stiffness matrix");
 
   for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) {
-
-    Array<UInt> & elem_filter = element_filter(type, ghost_type);
-    if (elem_filter.size() == 0) {
+    auto & elem_filter = element_filter(type, ghost_type);
+    if (elem_filter.empty()) {
       AKANTU_DEBUG_OUT();
       return;
     }
-    
+
     auto nb_element = elem_filter.size();
     auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
     auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type);
 
-
     auto nt_b_n = std::make_unique<Array<Real>>(
         nb_element * nb_quadrature_points,
         nb_nodes_per_element * nb_nodes_per_element, "N^t*b*N");
 
     auto bt_d_b = std::make_unique<Array<Real>>(
         nb_element * nb_quadrature_points,
         nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B");
 
     // damage_energy_density_on_qpoints = gc/l0 + phi = scalar
     auto & damage_energy_density_vect = damage_energy_density(type, ghost_type);
 
     // damage_energy_on_qpoints = gc*l0 = scalar
     auto & damage_energy_vect = damage_energy(type, ghost_type);
 
     fem.computeBtDB(damage_energy_vect, *bt_d_b, 2, type, ghost_type,
                     elem_filter);
 
     fem.computeNtbN(damage_energy_density_vect, *nt_b_n, type, ghost_type,
                     elem_filter);
 
     /// compute @f$ K_{\grad d} = \int_e \mathbf{N}^t * \mathbf{w} *
     /// \mathbf{N}@f$
     auto K_n = std::make_unique<Array<Real>>(
         nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_n");
 
     fem.integrate(*nt_b_n, *K_n, nb_nodes_per_element * nb_nodes_per_element,
                   type, ghost_type, elem_filter);
 
     model.getDOFManager().assembleElementalMatricesToMatrix(
-		  "K", "damage", *K_n, type, _not_ghost, _symmetric, elem_filter);
+        "K", "damage", *K_n, type, _not_ghost, _symmetric, elem_filter);
 
     /// compute @f$ K_{\grad d} = \int_e \mathbf{B}^t * \mathbf{W} *
     /// \mathbf{B}@f$
     auto K_b = std::make_unique<Array<Real>>(
         nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_b");
 
     fem.integrate(*bt_d_b, *K_b, nb_nodes_per_element * nb_nodes_per_element,
                   type, ghost_type, elem_filter);
 
     model.getDOFManager().assembleElementalMatricesToMatrix(
-		 "K", "damage", *K_b, type, _not_ghost, _symmetric, elem_filter);
+        "K", "damage", *K_b, type, _not_ghost, _symmetric, elem_filter);
   }
 
   AKANTU_DEBUG_OUT();
-}  
+}
 
 /* -------------------------------------------------------------------------- */
 void PhaseField::beforeSolveStep() {
   this->savePreviousState();
   this->computeAllDrivingForces(_not_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
-void PhaseField::afterSolveStep() {
+void PhaseField::afterSolveStep() {}
 
-}
-  
 /* -------------------------------------------------------------------------- */
 void PhaseField::savePreviousState() {
   AKANTU_DEBUG_IN();
 
-  for (auto pair : internal_vectors_real)
-    if (pair.second->hasHistory())
+  for (auto pair : internal_vectors_real) {
+    if (pair.second->hasHistory()) {
       pair.second->saveCurrentValues();
+    }
+  }
 
   AKANTU_DEBUG_OUT();
 }
-  
-  
+
 /* -------------------------------------------------------------------------- */
 void PhaseField::printself(std::ostream & stream, int indent) const {
- std::string space(indent, AKANTU_INDENT);
+  std::string space(indent, AKANTU_INDENT);
   std::string type = getID().substr(getID().find_last_of(':') + 1);
 
   stream << space << "PhaseField Material " << type << " [" << std::endl;
   Parsable::printself(stream, indent);
   stream << space << "]" << std::endl;
 }
 
-  
-}
+} // namespace akantu