diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh
index 3fd0b247b..e8985a5a9 100644
--- a/src/model/solid_mechanics/material.hh
+++ b/src/model/solid_mechanics/material.hh
@@ -1,718 +1,718 @@
 /**
  * @file   material.hh
  *
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  Mother class for all materials
  *
  *
  * 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_factory.hh"
 #include "aka_voigthelper.hh"
 #include "data_accessor.hh"
 #include "integration_point.hh"
 #include "parsable.hh"
 #include "parser.hh"
 /* -------------------------------------------------------------------------- */
 #include "internal_field.hh"
 #include "random_internal_field.hh"
 /* -------------------------------------------------------------------------- */
 #include "mesh_events.hh"
 #include "solid_mechanics_model_event_handler.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_MATERIAL_HH_
 #define AKANTU_MATERIAL_HH_
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 class Model;
 class SolidMechanicsModel;
 class Material;
 } // namespace akantu
 
 namespace akantu {
 
 using MaterialFactory =
     Factory<Material, ID, UInt, const ID &, SolidMechanicsModel &, const ID &>;
 
 /**
  * Interface of all materials
  * Prerequisites for a new material
  * - inherit from this class
  * - implement the following methods:
  * \code
  *  virtual Real getStableTimeStep(Real h, const Element & element =
  * ElementNull);
  *
  *  virtual void computeStress(ElementType el_type,
  *                             GhostType ghost_type = _not_ghost);
  *
  *  virtual void computeTangentStiffness(ElementType el_type,
  *                                       Array<Real> & tangent_matrix,
  *                                       GhostType ghost_type = _not_ghost);
  * \endcode
  *
  */
 class Material : public DataAccessor<Element>,
                  public Parsable,
                  public MeshEventHandler,
                  protected SolidMechanicsModelEventHandler {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   Material(const Material & mat) = delete;
   Material & operator=(const Material & mat) = delete;
 
   /// Initialize material with defaults
   Material(SolidMechanicsModel & model, const ID & id = "");
 
   /// Initialize material with custom mesh & fe_engine
   Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh,
            FEEngine & fe_engine, const ID & id = "");
 
   /// Destructor
   ~Material() override;
 
 protected:
   void initialize();
 
   /* ------------------------------------------------------------------------ */
   /* Function that materials can/should reimplement                           */
   /* ------------------------------------------------------------------------ */
 protected:
   /// constitutive law
   virtual void computeStress(ElementType /* el_type */,
                              GhostType /* ghost_type */ = _not_ghost) {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// compute the tangent stiffness matrix
   virtual void computeTangentModuli(ElementType /*el_type*/,
                                     Array<Real> & /*tangent_matrix*/,
                                     GhostType /*ghost_type*/ = _not_ghost) {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// compute the potential energy
   virtual void computePotentialEnergy(ElementType el_type);
 
   /// compute the potential energy for an element
   virtual void
   computePotentialEnergyByElement(ElementType /*type*/, UInt /*index*/,
                                   Vector<Real> & /*epot_on_quad_points*/) {
     AKANTU_TO_IMPLEMENT();
   }
 
   virtual void updateEnergies(ElementType /*el_type*/) {}
 
   virtual void updateEnergiesAfterDamage(ElementType /*el_type*/) {}
 
   /// set the material to steady state (to be implemented for materials that
   /// need it)
   virtual void setToSteadyState(ElementType /*el_type*/,
                                 GhostType /*ghost_type*/ = _not_ghost) {}
 
   /// function called to update the internal parameters when the modifiable
   /// parameters are modified
   virtual void updateInternalParameters() {}
 
 public:
   /// extrapolate internal values
   virtual void extrapolateInternal(const ID & id, const Element & element,
                                    const Matrix<Real> & points,
                                    Matrix<Real> & extrapolated);
 
   /// compute the p-wave speed in the material
   virtual Real getPushWaveSpeed(const Element & /*element*/) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// compute the s-wave speed in the material
   virtual Real getShearWaveSpeed(const Element & /*element*/) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// get a material celerity to compute the stable time step (default: is the
   /// push wave speed)
   virtual Real getCelerity(const Element & element) const {
     return getPushWaveSpeed(element);
   }
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   template <typename T> void registerInternal(InternalField<T> & /*vect*/) {
     AKANTU_TO_IMPLEMENT();
   }
 
   template <typename T> void unregisterInternal(InternalField<T> & /*vect*/) {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// initialize the material computed parameter
   virtual void initMaterial();
 
   /// compute the residual for this material
   //  virtual void updateResidual(GhostType ghost_type = _not_ghost);
 
   /// assemble the residual for this material
   virtual void assembleInternalForces(GhostType ghost_type);
 
   /// save the stress in the previous_stress if needed
   virtual void savePreviousState();
 
   /// restore the stress from previous_stress if needed
   virtual void restorePreviousState();
 
   /// compute the stresses for this material
   virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
   // virtual void
   // computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost);
   virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost);
 
   /// set material to steady state
   void setToSteadyState(GhostType ghost_type = _not_ghost);
 
   /// compute the stiffness matrix
   virtual void assembleStiffnessMatrix(GhostType ghost_type);
 
   /// add an element to the local mesh filter
   inline UInt addElement(ElementType type, UInt element, GhostType ghost_type);
   inline UInt addElement(const Element & element);
 
   /// add many elements at once
   void addElements(const Array<Element> & elements_to_add);
 
   /// remove many element at once
   void removeElements(const Array<Element> & elements_to_remove);
 
   /// function to print the contain of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   /**
    * interpolate stress on given positions for each element by means
    * of a geometrical interpolation on quadrature points
    */
   void interpolateStress(ElementTypeMapArray<Real> & result,
                          GhostType ghost_type = _not_ghost);
 
   /**
    * interpolate stress on given positions for each element by means
    * of a geometrical interpolation on quadrature points and store the
    * results per facet
    */
   void interpolateStressOnFacets(ElementTypeMapArray<Real> & result,
                                  ElementTypeMapArray<Real> & by_elem_result,
                                  GhostType ghost_type = _not_ghost);
 
   /**
    * function to initialize the elemental field interpolation
    * function by inverting the quadrature points' coordinates
    */
   void initElementalFieldInterpolation(
       const ElementTypeMapArray<Real> & interpolation_points_coordinates);
 
   /* ------------------------------------------------------------------------ */
   /* Common part                                                              */
   /* ------------------------------------------------------------------------ */
 protected:
   /* ------------------------------------------------------------------------ */
   static inline UInt getTangentStiffnessVoigtSize(UInt dim);
 
   /// compute the potential energy by element
   void computePotentialEnergyByElements();
 
   /// resize the intenals arrays
   virtual void resizeInternals();
 
   /* ------------------------------------------------------------------------ */
   /* Finite deformation functions                                             */
   /* This functions area implementing what is described in the paper of Bathe */
   /* et al, in IJNME, Finite Element Formulations For Large Deformation       */
   /* Dynamic Analysis, Vol 9, 353-386, 1975                                   */
   /* ------------------------------------------------------------------------ */
 protected:
   /// assemble the residual
   template <UInt dim> void assembleInternalForces(GhostType ghost_type);
 
   template <UInt dim>
   void computeAllStressesFromTangentModuli(ElementType type,
                                            GhostType ghost_type);
 
   template <UInt dim>
   void assembleStiffnessMatrix(ElementType type, GhostType ghost_type);
 
   /// assembling in finite deformation
   template <UInt dim>
   void assembleStiffnessMatrixNL(ElementType type, GhostType ghost_type);
 
   template <UInt dim>
   void assembleStiffnessMatrixL2(ElementType type, GhostType ghost_type);
 
   /* ------------------------------------------------------------------------ */
   /* Conversion functions                                                     */
   /* ------------------------------------------------------------------------ */
 public:
   /// Size of the Stress matrix for the case of finite deformation see: Bathe et
   /// al, IJNME, Vol 9, 353-386, 1975
   static inline UInt getCauchyStressMatrixSize(UInt dim);
 
   /// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386,
   /// 1975
   template <UInt dim>
   static inline void setCauchyStressMatrix(const Matrix<Real> & S_t,
                                            Matrix<Real> & sigma);
 
   /// write the stress tensor in the Voigt notation.
   template <UInt dim>
   static inline decltype(auto) stressToVoigt(const Matrix<Real> & stress) {
     return VoigtHelper<dim>::matrixToVoigt(stress);
   }
 
   /// write the strain tensor in the Voigt notation.
   template <UInt dim>
   static inline decltype(auto) strainToVoigt(const Matrix<Real> & strain) {
     return VoigtHelper<dim>::matrixToVoigtWithFactors(strain);
   }
 
   /// write a voigt vector to stress
   template <UInt dim>
   static inline void voigtToStress(const Vector<Real> & voigt,
                                    Matrix<Real> & stress) {
     VoigtHelper<dim>::voigtToMatrix(voigt, stress);
   }
 
   /// Computation of Cauchy stress tensor in the case of finite deformation from
   /// the 2nd Piola-Kirchhoff for a given element type
   template <UInt dim>
   void StoCauchy(ElementType el_type, GhostType ghost_type = _not_ghost);
 
   /// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation
   /// gradient
   template <UInt dim>
   inline void StoCauchy(const Matrix<Real> & F, const Matrix<Real> & S,
                         Matrix<Real> & sigma, const Real & C33 = 1.0) const;
 
   template <UInt dim>
   static inline void gradUToF(const Matrix<Real> & grad_u, Matrix<Real> & F);
 
   template <UInt dim>
   static inline decltype(auto) gradUToF(const Matrix<Real> & grad_u);
 
   static inline void rightCauchy(const Matrix<Real> & F, Matrix<Real> & C);
   static inline void leftCauchy(const Matrix<Real> & F, Matrix<Real> & B);
 
   template <UInt dim>
   static inline void gradUToEpsilon(const Matrix<Real> & grad_u,
                                     Matrix<Real> & epsilon);
   template <UInt dim>
   static inline decltype(auto) gradUToEpsilon(const Matrix<Real> & grad_u);
 
   template <UInt dim>
   static inline void gradUToE(const Matrix<Real> & grad_u,
                               Matrix<Real> & epsilon);
 
   template <UInt dim>
   static inline decltype(auto) gradUToE(const Matrix<Real> & grad_u);
 
   static inline Real stressToVonMises(const Matrix<Real> & stress);
 
 protected:
   /// converts global element to local element
   inline Element convertToLocalElement(const Element & global_element) const;
   /// converts local element to global element
   inline Element convertToGlobalElement(const Element & local_element) const;
 
   /// converts global quadrature point to local quadrature point
   inline IntegrationPoint
   convertToLocalPoint(const IntegrationPoint & global_point) const;
   /// converts local quadrature point to global quadrature point
   inline IntegrationPoint
   convertToGlobalPoint(const IntegrationPoint & local_point) const;
 
   /* ------------------------------------------------------------------------ */
   /* DataAccessor inherited members                                           */
   /* ------------------------------------------------------------------------ */
 public:
   inline UInt getNbData(const Array<Element> & elements,
                         const SynchronizationTag & tag) const override;
 
   inline void packData(CommunicationBuffer & buffer,
                        const Array<Element> & elements,
                        const SynchronizationTag & tag) const override;
 
   inline void unpackData(CommunicationBuffer & buffer,
                          const Array<Element> & elements,
                          const SynchronizationTag & tag) override;
 
   template <typename T>
   inline void packElementDataHelper(const ElementTypeMapArray<T> & data_to_pack,
                                     CommunicationBuffer & buffer,
                                     const Array<Element> & elements,
                                     const ID & fem_id = ID()) const;
 
   template <typename T>
   inline void unpackElementDataHelper(ElementTypeMapArray<T> & data_to_unpack,
                                       CommunicationBuffer & buffer,
                                       const Array<Element> & elements,
                                       const ID & fem_id = ID());
 
   /* ------------------------------------------------------------------------ */
   /* MeshEventHandler inherited members                                       */
   /* ------------------------------------------------------------------------ */
 public:
   /* ------------------------------------------------------------------------ */
   void onNodesAdded(const Array<UInt> & /*unused*/,
                     const NewNodesEvent & /*unused*/) override{};
   void onNodesRemoved(const Array<UInt> & /*unused*/,
                       const Array<UInt> & /*unused*/,
                       const RemovedNodesEvent & /*unused*/) override{};
   void onElementsAdded(const Array<Element> & element_list,
                        const NewElementsEvent & event) override;
   void onElementsRemoved(const Array<Element> & element_list,
                          const ElementTypeMapArray<UInt> & new_numbering,
                          const RemovedElementsEvent & event) override;
   void onElementsChanged(const Array<Element> & /*unused*/,
                          const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const ChangedElementsEvent & /*unused*/) override{};
 
   /* ------------------------------------------------------------------------ */
   /* SolidMechanicsModelEventHandler inherited members                        */
   /* ------------------------------------------------------------------------ */
 public:
   virtual void beforeSolveStep();
   virtual void afterSolveStep(bool converged = true);
 
   void onDamageIteration() override;
   void onDamageUpdate() override;
   void onDump() override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   AKANTU_GET_MACRO(Name, name, const std::string &);
 
   AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &)
 
   AKANTU_GET_MACRO(ID, id, const ID &);
   AKANTU_GET_MACRO(Rho, rho, Real);
   AKANTU_SET_MACRO(Rho, rho, Real);
 
   AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
 
   /// tells if the material can compute energy energy_id
   inline bool hasEnergy(const ID & energy_id);
 
-  /// get the lists of possible energies
-  inline decltype(auto) getEnergyLists() { return (list_of_energies); }
+  /// get the list of possible energies
+  inline decltype(auto) getEnergiesList() { return (list_of_energies); }
 
   /// return the potential energy for the subset of elements contained by the
   /// material
   Real getPotentialEnergy();
   /// return the potential energy for the provided element
   Real getPotentialEnergy(ElementType & type, UInt index);
 
   /// return the energy (identified by id) for the subset of elements contained
   /// by the material
   virtual Real getEnergy(const ID & energy_id);
   /// return the energy (identified by id) for the provided element
   virtual Real getEnergy(const ID & energy_id, ElementType type, UInt index);
 
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy,
                                          Real);
   AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray<Real> &);
   AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray<Real> &);
   AKANTU_GET_MACRO(ElementFilter, element_filter,
                    const ElementTypeMapArray<UInt> &);
   AKANTU_GET_MACRO(FEEngine, fem, FEEngine &);
 
   bool isNonLocal() const { return is_non_local; }
 
   template <typename T>
   const Array<T> & getArray(const ID & id, ElementType type,
                             GhostType ghost_type = _not_ghost) const;
   template <typename T>
   Array<T> & getArray(const ID & id, ElementType type,
                       GhostType ghost_type = _not_ghost);
 
   template <typename T>
   const InternalField<T> & getInternal(const ID & id) const;
   template <typename T> InternalField<T> & getInternal(const ID & id);
 
   template <typename T>
   inline bool isInternal(const ID & id, ElementKind element_kind) const;
 
   template <typename T>
   ElementTypeMap<UInt> getInternalDataPerElem(const ID & id,
                                               ElementKind element_kind) const;
 
   bool isFiniteDeformation() const { return finite_deformation; }
   bool isInelasticDeformation() const { return inelastic_deformation; }
 
   template <typename T> inline void setParam(const ID & param, T value);
   inline const Parameter & getParam(const ID & param) const;
 
   template <typename T>
   void flattenInternal(const std::string & field_id,
                        ElementTypeMapArray<T> & internal_flat,
                        GhostType ghost_type = _not_ghost,
                        ElementKind element_kind = _ek_not_defined) const;
 
   /// apply a constant eigengrad_u everywhere in the material
   virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
                                GhostType /*ghost_type*/ = _not_ghost);
 
   bool hasMatrixChanged(const ID & id) {
     if (id == "K") {
       return hasStiffnessMatrixChanged() or finite_deformation;
     }
 
     return true;
   }
 
   MatrixType getMatrixType(const ID & id) {
     if (id == "K") {
       return getTangentType();
     }
 
     if (id == "M") {
       return _symmetric;
     }
 
     return _mt_not_defined;
   }
 
   /// specify if the matrix need to be recomputed for this material
   virtual bool hasStiffnessMatrixChanged() { return true; }
 
   /// specify the type of matrix, if not overloaded the material is not valid
   /// for static or implicit computations
   virtual MatrixType getTangentType() { return _mt_not_defined; }
 
   /// static method to reteive the material factory
   static MaterialFactory & getFactory();
 
 protected:
   bool isInit() const { return is_init; }
 
   /// register the fact that material can compute a given energy
   void registerEnergy(const ID & energy_id);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// boolean to know if the material has been initialized
   bool is_init{false};
 
   std::map<ID, InternalField<Real> *> internal_vectors_real;
   std::map<ID, InternalField<UInt> *> internal_vectors_uint;
   std::map<ID, InternalField<bool> *> internal_vectors_bool;
 
 protected:
   ID id;
 
   /// Link to the fem object in the model
   FEEngine & fem;
 
   /// Finite deformation
   bool finite_deformation{false};
 
   /// Finite deformation
   bool inelastic_deformation{false};
 
   /// material name
   std::string name;
 
   /// The model to witch the material belong
   SolidMechanicsModel & model;
 
   /// density : rho
   Real rho{0.};
 
   /// spatial dimension
   UInt spatial_dimension;
 
   /// list of element handled by the material
   ElementTypeMapArray<UInt> element_filter;
 
   /// stresses arrays ordered by element types
   InternalField<Real> stress;
 
   /// eigengrad_u arrays ordered by element types
   InternalField<Real> eigengradu;
 
   /// grad_u arrays ordered by element types
   InternalField<Real> gradu;
 
   /// Green Lagrange strain (Finite deformation)
   InternalField<Real> green_strain;
 
   /// Second Piola-Kirchhoff stress tensor arrays ordered by element types
   /// (Finite deformation)
   InternalField<Real> piola_kirchhoff_2;
 
   /// potential energy by element
   InternalField<Real> potential_energy;
 
   /// tell if using in non local mode or not
   bool is_non_local{false};
 
   /// tell if the material need the previous stress state
   bool use_previous_stress{false};
 
   /// tell if the material need the previous strain state
   bool use_previous_gradu{false};
 
   /// elemental field interpolation coordinates
   InternalField<Real> interpolation_inverse_coordinates;
 
   /// elemental field interpolation points
   InternalField<Real> interpolation_points_matrices;
 
   /// vector that contains the names of all the internals that need to
   /// be transferred when material interfaces move
   std::vector<ID> internals_to_transfer;
 
 private:
   /// eigen_grad_u for the parser
   Matrix<Real> eigen_grad_u;
 
   /// list of usable energy
   std::set<ID> list_of_energies;
 };
 
 /// standard output stream operator
 inline std::ostream & operator<<(std::ostream & stream,
                                  const Material & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "material_inline_impl.hh"
 
 #include "internal_field_tmpl.hh"
 #include "random_internal_field_tmpl.hh"
 
 /* -------------------------------------------------------------------------- */
 /* Auto loop                                                                  */
 /* -------------------------------------------------------------------------- */
 /// This can be used to automatically write the loop on quadrature points in
 /// functions such as computeStress. This macro in addition to write the loop
 /// provides two tensors (matrices) sigma and grad_u
 #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type)       \
   auto && grad_u_view =                                                        \
       make_view(this->gradu(el_type, ghost_type), this->spatial_dimension,     \
                 this->spatial_dimension);                                      \
                                                                                \
   auto stress_view =                                                           \
       make_view(this->stress(el_type, ghost_type), this->spatial_dimension,    \
                 this->spatial_dimension);                                      \
                                                                                \
   if (this->isFiniteDeformation()) {                                           \
     stress_view = make_view(this->piola_kirchhoff_2(el_type, ghost_type),      \
                             this->spatial_dimension, this->spatial_dimension); \
   }                                                                            \
                                                                                \
   for (auto && data : zip(grad_u_view, stress_view)) {                         \
     [[gnu::unused]] Matrix<Real> & grad_u = std::get<0>(data);                 \
     [[gnu::unused]] Matrix<Real> & sigma = std::get<1>(data)
 
 #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END }
 
 /// This can be used to automatically write the loop on quadrature points in
 /// functions such as computeTangentModuli. This macro in addition to write the
 /// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix
 /// where the elemental tangent moduli should be stored in Voigt Notation
 #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat)              \
   auto && grad_u_view =                                                        \
       make_view(this->gradu(el_type, ghost_type), this->spatial_dimension,     \
                 this->spatial_dimension);                                      \
                                                                                \
   auto && stress_view =                                                        \
       make_view(this->stress(el_type, ghost_type), this->spatial_dimension,    \
                 this->spatial_dimension);                                      \
                                                                                \
   auto tangent_size =                                                          \
       this->getTangentStiffnessVoigtSize(this->spatial_dimension);             \
                                                                                \
   auto && tangent_view = make_view(tangent_mat, tangent_size, tangent_size);   \
                                                                                \
   for (auto && data : zip(grad_u_view, stress_view, tangent_view)) {           \
     [[gnu::unused]] Matrix<Real> & grad_u = std::get<0>(data);                 \
     [[gnu::unused]] Matrix<Real> & sigma = std::get<1>(data);                  \
     Matrix<Real> & tangent = std::get<2>(data);
 
 #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END }
 
 /* -------------------------------------------------------------------------- */
 
 #define INSTANTIATE_MATERIAL_ONLY(mat_name)                                    \
   template class mat_name<1>; /* NOLINT */                                     \
   template class mat_name<2>; /* NOLINT */                                     \
   template class mat_name<3>  /* NOLINT */
 
 #define MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name)                       \
   [](UInt dim, const ID &, SolidMechanicsModel & model,                        \
      const ID & id) /* NOLINT */                                               \
       -> std::unique_ptr<                                                      \
           Material> { /* NOLINT */                                             \
                       switch (dim) {                                           \
                       case 1:                                                  \
                         return std::make_unique<mat_name<1>>(/* NOLINT */      \
                                                              model, id);       \
                       case 2:                                                  \
                         return std::make_unique<mat_name<2>>(/* NOLINT */      \
                                                              model, id);       \
                       case 3:                                                  \
                         return std::make_unique<mat_name<3>>(/* NOLINT */      \
                                                              model, id);       \
                       default:                                                 \
                         AKANTU_EXCEPTION(                                      \
                             "The dimension "                                   \
                             << dim                                             \
                             << "is not a valid dimension for the material "    \
                             << #id);                                           \
                       }                                                        \
   }
 
 #define INSTANTIATE_MATERIAL(id, mat_name)                                     \
   INSTANTIATE_MATERIAL_ONLY(mat_name);                                         \
   static bool material_is_alocated_##id [[gnu::unused]] =                      \
       MaterialFactory::getInstance().registerAllocator(                        \
           #id, MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name))
 
 #endif /* AKANTU_MATERIAL_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index b04015d0b..133d392e3 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1264 +1,1277 @@
 /**
  * @file   solid_mechanics_model.cc
  *
  * @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Clement Roux <clement.roux@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Tue Jul 27 2010
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  Implementation of the SolidMechanicsModel class
  *
  *
  * 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 "solid_mechanics_model.hh"
 #include "integrator_gauss.hh"
 #include "shape_lagrange.hh"
 #include "solid_mechanics_model_tmpl.hh"
 
 #include "communicator.hh"
 #include "element_synchronizer.hh"
 #include "sparse_matrix.hh"
 #include "synchronizer_registry.hh"
 
 #include "dumpable_inline_impl.hh"
 #ifdef AKANTU_USE_IOHELPER
 #include "dumper_iohelper_paraview.hh"
 #endif
 
 #include "material_non_local.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /**
  * A solid mechanics model need a mesh  and a dimension to be created. the model
  * by it  self can not  do a lot,  the good init  functions should be  called in
  * order to configure the model depending on what we want to do.
  *
  * @param  mesh mesh  representing  the model  we  want to  simulate
  * @param dim spatial  dimension of the problem, if dim =  0 (default value) the
  * dimension of the problem is assumed to be the on of the mesh
  * @param id an id to identify the model
  * @param model_type this is an internal parameter for inheritance purposes
  */
 
 SolidMechanicsModel::SolidMechanicsModel(
     Mesh & mesh, UInt dim, const ID & id,
     std::shared_ptr<DOFManager> dof_manager, const ModelType model_type)
     : Model(mesh, model_type, std::move(dof_manager), dim, id),
       material_index("material index", id),
       material_local_numbering("material local numbering", id) {
   AKANTU_DEBUG_IN();
 
   this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh,
                                                Model::spatial_dimension);
 
 #if defined(AKANTU_USE_IOHELPER)
   this->mesh.registerDumper<DumperParaview>("solid_mechanics_model", id, true);
   this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost,
                          _ek_regular);
 #endif
 
   material_selector = std::make_shared<DefaultMaterialSelector>(material_index);
 
   this->registerDataAccessor(*this);
 
   if (this->mesh.isDistributed()) {
     auto & synchronizer = this->mesh.getElementSynchronizer();
     this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id);
     this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass);
     this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress);
     this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 SolidMechanicsModel::~SolidMechanicsModel() = default;
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) {
   Model::setTimeStep(time_step, solver_id);
 
 #if defined(AKANTU_USE_IOHELPER)
   this->mesh.getDumper().setTimeStep(time_step);
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 /* Initialization                                                             */
 /* -------------------------------------------------------------------------- */
 /**
  * This function groups  many of the initialization in on  function. For most of
  * basics  case the  function should  be  enough. The  functions initialize  the
  * model, the internal  vectors, set them to 0, and  depending on the parameters
  * it also initialize the explicit or implicit solver.
  *
  * @param options
  * \parblock
  * contains the different options to initialize the model
  * \li \c analysis_method specify the type of solver to use
  * \endparblock
  */
 void SolidMechanicsModel::initFullImpl(const ModelOptions & options) {
   material_index.initialize(mesh, _element_kind = _ek_not_defined,
                             _default_value = UInt(-1), _with_nb_element = true);
   material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined,
                                       _with_nb_element = true);
 
   Model::initFullImpl(options);
 
   // initialize the materials
   if (not this->parser.getLastParsedFile().empty()) {
     this->instantiateMaterials();
     this->initMaterials();
   }
 
   this->initBC(*this, *displacement, *displacement_increment, *external_force);
 }
 
 /* -------------------------------------------------------------------------- */
 TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const {
   return TimeStepSolverType::_dynamic_lumped;
 }
 
 /* -------------------------------------------------------------------------- */
 ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions(
     const TimeStepSolverType & type) const {
   ModelSolverOptions options;
 
   switch (type) {
   case TimeStepSolverType::_dynamic_lumped: {
     options.non_linear_solver_type = NonLinearSolverType::_lumped;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_central_difference;
     options.solution_type["displacement"] = IntegrationScheme::_acceleration;
     break;
   }
   case TimeStepSolverType::_static: {
     options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_pseudo_time;
     options.solution_type["displacement"] = IntegrationScheme::_not_defined;
     break;
   }
   case TimeStepSolverType::_dynamic: {
     if (this->method == _explicit_consistent_mass) {
       options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
       options.integration_scheme_type["displacement"] =
           IntegrationSchemeType::_central_difference;
       options.solution_type["displacement"] = IntegrationScheme::_acceleration;
     } else {
       options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
       options.integration_scheme_type["displacement"] =
           IntegrationSchemeType::_trapezoidal_rule_2;
       options.solution_type["displacement"] = IntegrationScheme::_displacement;
     }
     break;
   }
   default:
     AKANTU_EXCEPTION(type << " is not a valid time step solver type");
   }
 
   return options;
 }
 
 /* -------------------------------------------------------------------------- */
 std::tuple<ID, TimeStepSolverType>
 SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
   switch (method) {
   case _explicit_lumped_mass: {
     return std::make_tuple("explicit_lumped",
                            TimeStepSolverType::_dynamic_lumped);
   }
   case _explicit_consistent_mass: {
     return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
   }
   case _static: {
     return std::make_tuple("static", TimeStepSolverType::_static);
   }
   case _implicit_dynamic: {
     return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
   }
   default:
     return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type,
                                      NonLinearSolverType /*unused*/) {
   auto & dof_manager = this->getDOFManager();
 
   /* ------------------------------------------------------------------------ */
   // for alloc type of solvers
   this->allocNodalField(this->displacement, spatial_dimension, "displacement");
   this->allocNodalField(this->previous_displacement, spatial_dimension,
                         "previous_displacement");
   this->allocNodalField(this->displacement_increment, spatial_dimension,
                         "displacement_increment");
   this->allocNodalField(this->internal_force, spatial_dimension,
                         "internal_force");
   this->allocNodalField(this->external_force, spatial_dimension,
                         "external_force");
   this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs");
   this->allocNodalField(this->current_position, spatial_dimension,
                         "current_position");
 
   // initialize the current positions
   this->current_position->copy(this->mesh.getNodes());
 
   /* ------------------------------------------------------------------------ */
   if (!dof_manager.hasDOFs("displacement")) {
     dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal);
     dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
     dof_manager.registerDOFsIncrement("displacement",
                                       *this->displacement_increment);
     dof_manager.registerDOFsPrevious("displacement",
                                      *this->previous_displacement);
   }
 
   /* ------------------------------------------------------------------------ */
   // for dynamic
   if (time_step_solver_type == TimeStepSolverType::_dynamic ||
       time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
     this->allocNodalField(this->velocity, spatial_dimension, "velocity");
     this->allocNodalField(this->acceleration, spatial_dimension,
                           "acceleration");
 
     if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
       dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
       dof_manager.registerDOFsDerivative("displacement", 2,
                                          *this->acceleration);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Initialize the model,basically it  pre-compute the shapes, shapes derivatives
  * and jacobian
  */
 void SolidMechanicsModel::initModel() {
   /// \todo add  the current position  as a parameter to  initShapeFunctions for
   /// large deformation
   getFEEngine().initShapeFunctions(_not_ghost);
   getFEEngine().initShapeFunctions(_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleResidual() {
   AKANTU_DEBUG_IN();
 
   /* ------------------------------------------------------------------------ */
   // computes the internal forces
   this->assembleInternalForces();
 
   /* ------------------------------------------------------------------------ */
   this->getDOFManager().assembleToResidual("displacement",
                                            *this->external_force, 1);
   this->getDOFManager().assembleToResidual("displacement",
                                            *this->internal_force, 1);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleResidual(const ID & residual_part) {
   AKANTU_DEBUG_IN();
 
   if ("external" == residual_part) {
     this->getDOFManager().assembleToResidual("displacement",
                                              *this->external_force, 1);
     AKANTU_DEBUG_OUT();
     return;
   }
 
   if ("internal" == residual_part) {
     this->assembleInternalForces();
     this->getDOFManager().assembleToResidual("displacement",
                                              *this->internal_force, 1);
     AKANTU_DEBUG_OUT();
     return;
   }
 
   AKANTU_CUSTOM_EXCEPTION(
       debug::SolverCallbackResidualPartUnknown(residual_part));
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) {
   // \TODO check the materials to know what is the correct answer
   if (matrix_id == "C") {
     return _mt_not_defined;
   }
 
   if (matrix_id == "K") {
     auto matrix_type = _unsymmetric;
 
     for (auto & material : materials) {
       matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id));
     }
   }
   return _symmetric;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) {
   if (matrix_id == "K") {
     this->assembleStiffnessMatrix();
   } else if (matrix_id == "M") {
     this->assembleMass();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) {
   if (matrix_id == "M") {
     this->assembleMassLumped();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::beforeSolveStep() {
   for (auto & material : materials) {
     material->beforeSolveStep();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::afterSolveStep(bool converged) {
   for (auto & material : materials) {
     material->afterSolveStep(converged);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::predictor() { ++displacement_release; }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::corrector() { ++displacement_release; }
 
 /* -------------------------------------------------------------------------- */
 /**
  * This function computes the internal forces as \f$F_{int} = \int_{\Omega} N
  * \sigma d\Omega@\f$
  */
 void SolidMechanicsModel::assembleInternalForces() {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_INFO("Assemble the internal forces");
 
   this->internal_force->zero();
 
   // compute the stresses of local elements
   AKANTU_DEBUG_INFO("Compute local stresses");
   for (auto & material : materials) {
     material->computeAllStresses(_not_ghost);
   }
 
   /* ------------------------------------------------------------------------ */
   /* Computation of the non local part */
   if (this->non_local_manager) {
     this->non_local_manager->computeAllNonLocalStresses();
   }
 
   // communicate the stresses
   AKANTU_DEBUG_INFO("Send data for residual assembly");
   this->asynchronousSynchronize(SynchronizationTag::_smm_stress);
 
   // assemble the forces due to local stresses
   AKANTU_DEBUG_INFO("Assemble residual for local elements");
   for (auto & material : materials) {
     material->assembleInternalForces(_not_ghost);
   }
 
   // finalize communications
   AKANTU_DEBUG_INFO("Wait distant stresses");
   this->waitEndSynchronize(SynchronizationTag::_smm_stress);
 
   // assemble the stresses due to ghost elements
   AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
   for (auto & material : materials) {
     material->assembleInternalForces(_ghost);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleStiffnessMatrix(bool need_to_reassemble) {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
 
   // Check if materials need to recompute the matrix
   for (auto & material : materials) {
     need_to_reassemble |= material->hasMatrixChanged("K");
   }
 
   if (need_to_reassemble) {
     this->getDOFManager().getMatrix("K").zero();
 
     // call compute stiffness matrix on each local elements
     for (auto & material : materials) {
       material->assembleStiffnessMatrix(_not_ghost);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateCurrentPosition() {
   if (this->current_position_release == this->displacement_release) {
     return;
   }
 
   this->current_position->copy(this->mesh.getNodes());
 
   auto cpos_it = this->current_position->begin(Model::spatial_dimension);
   auto cpos_end = this->current_position->end(Model::spatial_dimension);
   auto disp_it = this->displacement->begin(Model::spatial_dimension);
 
   for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) {
     *cpos_it += *disp_it;
   }
 
   this->current_position_release = this->displacement_release;
 }
 
 /* -------------------------------------------------------------------------- */
 const Array<Real> & SolidMechanicsModel::getCurrentPosition() {
   this->updateCurrentPosition();
   return *this->current_position;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateDataForNonLocalCriterion(
     ElementTypeMapReal & criterion) {
   const ID field_name = criterion.getName();
   for (auto & material : materials) {
     if (!material->isInternal<Real>(field_name, _ek_regular)) {
       continue;
     }
 
     for (auto ghost_type : ghost_types) {
       material->flattenInternal(field_name, criterion, ghost_type, _ek_regular);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 /* Information                                                                */
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getStableTimeStep() {
   AKANTU_DEBUG_IN();
 
   Real min_dt = getStableTimeStep(_not_ghost);
 
   /// reduction min over all processors
   mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min);
 
   AKANTU_DEBUG_OUT();
   return min_dt;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getStableTimeStep(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Real min_dt = std::numeric_limits<Real>::max();
 
   this->updateCurrentPosition();
 
   Element elem;
   elem.ghost_type = ghost_type;
 
   for (auto type :
        mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) {
     elem.type = type;
     UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
     UInt nb_element = mesh.getNbElement(type);
 
     auto mat_indexes = material_index(type, ghost_type).begin();
     auto mat_loc_num = material_local_numbering(type, ghost_type).begin();
 
     Array<Real> X(0, nb_nodes_per_element * Model::spatial_dimension);
     FEEngine::extractNodalToElementField(mesh, *current_position, X, type,
                                          _not_ghost);
 
     auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element);
 
     for (UInt el = 0; el < nb_element;
          ++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
       elem.element = *mat_loc_num;
       Real el_h = getFEEngine().getElementInradius(*X_el, type);
       Real el_c = this->materials[*mat_indexes]->getCelerity(elem);
       Real el_dt = el_h / el_c;
 
       min_dt = std::min(min_dt, el_dt);
     }
   }
 
   AKANTU_DEBUG_OUT();
   return min_dt;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getKineticEnergy() {
   AKANTU_DEBUG_IN();
 
   Real ekin = 0.;
   UInt nb_nodes = mesh.getNbNodes();
 
   if (this->getDOFManager().hasLumpedMatrix("M")) {
     auto m_it = this->mass->begin(Model::spatial_dimension);
     auto m_end = this->mass->end(Model::spatial_dimension);
     auto v_it = this->velocity->begin(Model::spatial_dimension);
 
     for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) {
       const auto & v = *v_it;
       const auto & m = *m_it;
 
       Real mv2 = 0.;
       auto is_local_node = mesh.isLocalOrMasterNode(n);
       // bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
       auto count_node = is_local_node; // && is_not_pbc_slave_node;
       if (count_node) {
         for (UInt i = 0; i < Model::spatial_dimension; ++i) {
           if (m(i) > std::numeric_limits<Real>::epsilon()) {
             mv2 += v(i) * v(i) * m(i);
           }
         }
       }
 
       ekin += mv2;
     }
   } else if (this->getDOFManager().hasMatrix("M")) {
     Array<Real> Mv(nb_nodes, Model::spatial_dimension);
     this->getDOFManager().assembleMatMulVectToArray("displacement", "M",
                                                     *this->velocity, Mv);
 
     for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension),
                             make_view(*this->velocity, spatial_dimension))) {
       ekin += std::get<2>(data).dot(std::get<1>(data)) *
               static_cast<Real>(mesh.isLocalOrMasterNode(std::get<0>(data)));
     }
   } else {
     AKANTU_ERROR("No function called to assemble the mass matrix.");
   }
 
   mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum);
 
   AKANTU_DEBUG_OUT();
   return ekin * .5;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getKineticEnergy(ElementType type, UInt index) {
   AKANTU_DEBUG_IN();
 
   UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
 
   Array<Real> vel_on_quad(nb_quadrature_points, Model::spatial_dimension);
   Array<UInt> filter_element(1, 1, index);
 
   getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad,
                                                Model::spatial_dimension, type,
                                                _not_ghost, filter_element);
 
   auto vit = vel_on_quad.begin(Model::spatial_dimension);
   auto vend = vel_on_quad.end(Model::spatial_dimension);
 
   Vector<Real> rho_v2(nb_quadrature_points);
 
   Real rho = materials[material_index(type)(index)]->getRho();
 
   for (UInt q = 0; vit != vend; ++vit, ++q) {
     rho_v2(q) = rho * vit->dot(*vit);
   }
 
   AKANTU_DEBUG_OUT();
 
   return .5 * getFEEngine().integrate(rho_v2, type, index);
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getExternalWork() {
   AKANTU_DEBUG_IN();
 
   auto ext_force_it = external_force->begin(Model::spatial_dimension);
   auto int_force_it = internal_force->begin(Model::spatial_dimension);
   auto boun_it = blocked_dofs->begin(Model::spatial_dimension);
 
   decltype(ext_force_it) incr_or_velo_it;
   if (this->method == _static) {
     incr_or_velo_it =
         this->displacement_increment->begin(Model::spatial_dimension);
   } else {
     incr_or_velo_it = this->velocity->begin(Model::spatial_dimension);
   }
 
   Real work = 0.;
 
   UInt nb_nodes = this->mesh.getNbNodes();
 
   for (UInt n = 0; n < nb_nodes;
        ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) {
     const auto & int_force = *int_force_it;
     const auto & ext_force = *ext_force_it;
     const auto & boun = *boun_it;
     const auto & incr_or_velo = *incr_or_velo_it;
 
     bool is_local_node = this->mesh.isLocalOrMasterNode(n);
     // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n);
     bool count_node = is_local_node; // && is_not_pbc_slave_node;
 
     if (count_node) {
       for (UInt i = 0; i < Model::spatial_dimension; ++i) {
         if (boun(i)) {
           work -= int_force(i) * incr_or_velo(i);
         } else {
           work += ext_force(i) * incr_or_velo(i);
         }
       }
     }
   }
 
   mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum);
 
   if (this->method != _static) {
     work *= this->getTimeStep();
   }
 
   AKANTU_DEBUG_OUT();
   return work;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getEnergy(const ID & energy_id) {
   AKANTU_DEBUG_IN();
 
   if (energy_id == "kinetic") {
     return getKineticEnergy();
   }
 
   if (energy_id == "external work") {
     return getExternalWork();
   }
 
   std::set<ID> list_of_energies;
   for (auto && material : materials) {
-    for(auto && energy : material->getEnergyLists()) {
+    for (auto && energy : material->getEnergiesList()) {
       list_of_energies.insert(energy);
     }
   }
 
-  if(list_of_energies.find(energy_id) != list_of_energies.end()) {
+  if (list_of_energies.find(energy_id) != list_of_energies.end()) {
+    list_of_energies.insert("kinetic");
+    list_of_energies.insert("external work");
 
+    std::string str_list;
+    for (auto it = list_of_energies.begin(); it != list_of_energies.end();
+         ++it) {
+      if (it != list_of_energies.begin()) {
+        str_list += ", "
+      }
+      str_list += "\"" + *it + "\"";
+    }
+
+    AKANTU_EXCEPTION("The energy \""
+                     << energy_id << "\" does not exists, valid energies are ["
+                     << str_list << "]");
   }
 
   Real energy = 0.;
   for (auto & material : materials) {
     energy += material->getEnergy(energy_id);
   }
 
   /// reduction sum over all processors
   mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
 
   AKANTU_DEBUG_OUT();
   return energy;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
                                     ElementType type, UInt index) {
   AKANTU_DEBUG_IN();
 
   if (energy_id == "kinetic") {
     return getKineticEnergy(type, index);
   }
 
   UInt mat_index = this->material_index(type, _not_ghost)(index);
   UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
   Real energy =
       this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
 
   AKANTU_DEBUG_OUT();
   return energy;
 }
 
 /* -------------------------------------------------------------------------- */
-Real SolidMechanicsModel::getEnergy(const ID & energy_id,
-                                    const ID & group_id) {
+Real SolidMechanicsModel::getEnergy(const ID & energy_id, const ID & group_id) {
   auto && group = mesh.getElementGroup(group_id);
   auto energy = 0.;
-  for(auto && type : group.elementTypes()) {
-    for(auto el : group.getElementsIterable(type)) {
+  for (auto && type : group.elementTypes()) {
+    for (auto el : group.getElementsIterable(type)) {
       energy += getEnergy(energy_id, el);
     }
   }
 
   /// reduction sum over all processors
   mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
 
   return energy;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
                                           const NewElementsEvent & event) {
   AKANTU_DEBUG_IN();
 
   this->material_index.initialize(mesh, _element_kind = _ek_not_defined,
                                   _with_nb_element = true,
                                   _default_value = UInt(-1));
   this->material_local_numbering.initialize(
       mesh, _element_kind = _ek_not_defined, _with_nb_element = true,
       _default_value = UInt(-1));
 
   ElementTypeMapArray<UInt> filter("new_element_filter", this->getID());
 
   for (const auto & elem : element_list) {
     if (mesh.getSpatialDimension(elem.type) != spatial_dimension) {
       continue;
     }
 
     if (!filter.exists(elem.type, elem.ghost_type)) {
       filter.alloc(0, 1, elem.type, elem.ghost_type);
     }
     filter(elem.type, elem.ghost_type).push_back(elem.element);
   }
 
   // this fails in parallel if the event is sent on facet between constructor
   // and initFull \todo: to debug...
   this->assignMaterialToElements(&filter);
 
   for (auto & material : materials) {
     material->onElementsAdded(element_list, event);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onElementsRemoved(
     const Array<Element> & element_list,
     const ElementTypeMapArray<UInt> & new_numbering,
     const RemovedElementsEvent & event) {
   for (auto & material : materials) {
     material->onElementsRemoved(element_list, new_numbering, event);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
                                        const NewNodesEvent & event) {
   AKANTU_DEBUG_IN();
   UInt nb_nodes = mesh.getNbNodes();
 
   if (displacement) {
     displacement->resize(nb_nodes, 0.);
     ++displacement_release;
   }
   if (mass) {
     mass->resize(nb_nodes, 0.);
   }
   if (velocity) {
     velocity->resize(nb_nodes, 0.);
   }
   if (acceleration) {
     acceleration->resize(nb_nodes, 0.);
   }
   if (external_force) {
     external_force->resize(nb_nodes, 0.);
   }
   if (internal_force) {
     internal_force->resize(nb_nodes, 0.);
   }
   if (blocked_dofs) {
     blocked_dofs->resize(nb_nodes, false);
   }
   if (current_position) {
     current_position->resize(nb_nodes, 0.);
   }
 
   if (previous_displacement) {
     previous_displacement->resize(nb_nodes, 0.);
   }
   if (displacement_increment) {
     displacement_increment->resize(nb_nodes, 0.);
   }
 
   for (auto & material : materials) {
     material->onNodesAdded(nodes_list, event);
   }
 
   need_to_reassemble_lumped_mass = true;
   need_to_reassemble_mass = true;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onNodesRemoved(const Array<UInt> & /*element_list*/,
                                          const Array<UInt> & new_numbering,
                                          const RemovedNodesEvent & /*event*/) {
   if (displacement) {
     mesh.removeNodesFromArray(*displacement, new_numbering);
     ++displacement_release;
   }
   if (mass) {
     mesh.removeNodesFromArray(*mass, new_numbering);
   }
   if (velocity) {
     mesh.removeNodesFromArray(*velocity, new_numbering);
   }
   if (acceleration) {
     mesh.removeNodesFromArray(*acceleration, new_numbering);
   }
   if (internal_force) {
     mesh.removeNodesFromArray(*internal_force, new_numbering);
   }
   if (external_force) {
     mesh.removeNodesFromArray(*external_force, new_numbering);
   }
   if (blocked_dofs) {
     mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
   }
 
   // if (increment_acceleration)
   //   mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
   if (displacement_increment) {
     mesh.removeNodesFromArray(*displacement_increment, new_numbering);
   }
 
   if (previous_displacement) {
     mesh.removeNodesFromArray(*previous_displacement, new_numbering);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
   std::string space(indent, AKANTU_INDENT);
 
   stream << space << "Solid Mechanics Model [" << std::endl;
   stream << space << " + id                : " << id << std::endl;
   stream << space << " + spatial dimension : " << Model::spatial_dimension
          << std::endl;
 
   stream << space << " + fem [" << std::endl;
   getFEEngine().printself(stream, indent + 2);
   stream << space << " ]" << std::endl;
 
   stream << space << " + nodals information [" << std::endl;
   displacement->printself(stream, indent + 2);
   if (velocity) {
     velocity->printself(stream, indent + 2);
   }
   if (acceleration) {
     acceleration->printself(stream, indent + 2);
   }
   if (mass) {
     mass->printself(stream, indent + 2);
   }
   external_force->printself(stream, indent + 2);
   internal_force->printself(stream, indent + 2);
   blocked_dofs->printself(stream, indent + 2);
   stream << space << " ]" << std::endl;
 
   stream << space << " + material information [" << std::endl;
   material_index.printself(stream, indent + 2);
   stream << space << " ]" << std::endl;
 
   stream << space << " + materials [" << std::endl;
   for (const auto & material : materials) {
     material->printself(stream, indent + 2);
   }
   stream << space << " ]" << std::endl;
 
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::initializeNonLocal() {
   this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods(
     GhostType ghost_type) {
   for (auto & mat : materials) {
     MaterialNonLocalInterface * mat_non_local;
     if ((mat_non_local =
              dynamic_cast<MaterialNonLocalInterface *>(mat.get())) == nullptr) {
       continue;
     }
 
     ElementTypeMapArray<Real> quadrature_points_coordinates(
         "quadrature_points_coordinates_tmp_nl", this->id);
     quadrature_points_coordinates.initialize(this->getFEEngine(),
                                              _nb_component = spatial_dimension,
                                              _ghost_type = ghost_type);
 
     for (const auto & type : quadrature_points_coordinates.elementTypes(
              Model::spatial_dimension, ghost_type)) {
       this->getFEEngine().computeIntegrationPointsCoordinates(
           quadrature_points_coordinates(type, ghost_type), type, ghost_type);
     }
 
     mat_non_local->initMaterialNonLocal();
 
     mat_non_local->insertIntegrationPointsInNeighborhoods(
         ghost_type, quadrature_points_coordinates);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::computeNonLocalStresses(GhostType ghost_type) {
   for (auto & mat : materials) {
     if (not aka::is_of_type<MaterialNonLocalInterface>(*mat)) {
       continue;
     }
 
     auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
     mat_non_local.computeNonLocalStresses(ghost_type);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateLocalInternal(
     ElementTypeMapReal & internal_flat, GhostType ghost_type,
     ElementKind kind) {
   const ID field_name = internal_flat.getName();
   for (auto & material : materials) {
     if (material->isInternal<Real>(field_name, kind)) {
       material->flattenInternal(field_name, internal_flat, ghost_type, kind);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateNonLocalInternal(
     ElementTypeMapReal & internal_flat, GhostType ghost_type,
     ElementKind kind) {
 
   const ID field_name = internal_flat.getName();
 
   for (auto & mat : materials) {
     if (not aka::is_of_type<MaterialNonLocalInterface>(*mat)) {
       continue;
     }
 
     auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
     mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type,
                                           kind);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) {
   return getFEEngineClassBoundary<MyFEEngineType>(name);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::splitElementByMaterial(
     const Array<Element> & elements,
     std::vector<Array<Element>> & elements_per_mat) const {
   for (const auto & el : elements) {
     Element mat_el = el;
     mat_el.element = this->material_local_numbering(el);
     elements_per_mat[this->material_index(el)].push_back(mat_el);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 UInt SolidMechanicsModel::getNbData(const Array<Element> & elements,
                                     const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   UInt size = 0;
   UInt nb_nodes_per_element = 0;
 
   for (const Element & el : elements) {
     nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
   }
 
   switch (tag) {
   case SynchronizationTag::_material_id: {
     size += elements.size() * sizeof(UInt);
     break;
   }
   case SynchronizationTag::_smm_mass: {
     size += nb_nodes_per_element * sizeof(Real) *
             Model::spatial_dimension; // mass vector
     break;
   }
   case SynchronizationTag::_smm_for_gradu: {
     size += nb_nodes_per_element * Model::spatial_dimension *
             sizeof(Real); // displacement
     break;
   }
   case SynchronizationTag::_smm_boundary: {
     // force, displacement, boundary
     size += nb_nodes_per_element * Model::spatial_dimension *
             (2 * sizeof(Real) + sizeof(bool));
     break;
   }
   case SynchronizationTag::_for_dump: {
     // displacement, velocity, acceleration, residual, force
     size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5;
     break;
   }
   default: {
   }
   }
 
   if (tag != SynchronizationTag::_material_id) {
     splitByMaterial(elements, [&](auto && mat, auto && elements) {
       size += mat.getNbData(elements, tag);
     });
   }
 
   AKANTU_DEBUG_OUT();
   return size;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
                                    const Array<Element> & elements,
                                    const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case SynchronizationTag::_material_id: {
-    packElementalDataHelper(
-        material_index, buffer, elements, false, getFEEngine());
+    packElementalDataHelper(material_index, buffer, elements, false,
+                            getFEEngine());
     break;
   }
   case SynchronizationTag::_smm_mass: {
     packNodalDataHelper(*mass, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_smm_for_gradu: {
     packNodalDataHelper(*displacement, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_for_dump: {
     packNodalDataHelper(*displacement, buffer, elements, mesh);
     packNodalDataHelper(*velocity, buffer, elements, mesh);
     packNodalDataHelper(*acceleration, buffer, elements, mesh);
     packNodalDataHelper(*internal_force, buffer, elements, mesh);
     packNodalDataHelper(*external_force, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_smm_boundary: {
     packNodalDataHelper(*external_force, buffer, elements, mesh);
     packNodalDataHelper(*velocity, buffer, elements, mesh);
     packNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
     break;
   }
   default: {
   }
   }
 
   if (tag != SynchronizationTag::_material_id) {
     splitByMaterial(elements, [&](auto && mat, auto && elements) {
       mat.packData(buffer, elements, tag);
     });
   }
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
                                      const Array<Element> & elements,
                                      const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case SynchronizationTag::_material_id: {
     for (auto && element : elements) {
       UInt recv_mat_index;
       buffer >> recv_mat_index;
       UInt & mat_index = material_index(element);
       if (mat_index != UInt(-1)) {
         continue;
       }
 
       // add ghosts element to the correct material
       mat_index = recv_mat_index;
       UInt index = materials[mat_index]->addElement(element);
       material_local_numbering(element) = index;
     }
     break;
   }
   case SynchronizationTag::_smm_mass: {
     unpackNodalDataHelper(*mass, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_smm_for_gradu: {
     unpackNodalDataHelper(*displacement, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_for_dump: {
     unpackNodalDataHelper(*displacement, buffer, elements, mesh);
     unpackNodalDataHelper(*velocity, buffer, elements, mesh);
     unpackNodalDataHelper(*acceleration, buffer, elements, mesh);
     unpackNodalDataHelper(*internal_force, buffer, elements, mesh);
     unpackNodalDataHelper(*external_force, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_smm_boundary: {
     unpackNodalDataHelper(*external_force, buffer, elements, mesh);
     unpackNodalDataHelper(*velocity, buffer, elements, mesh);
     unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
     break;
   }
   default: {
   }
   }
 
   if (tag != SynchronizationTag::_material_id) {
     splitByMaterial(elements, [&](auto && mat, auto && elements) {
       mat.unpackData(buffer, elements, tag);
     });
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 UInt SolidMechanicsModel::getNbData(const Array<UInt> & dofs,
                                     const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   UInt size = 0;
   //  UInt nb_nodes = mesh.getNbNodes();
 
   switch (tag) {
   case SynchronizationTag::_smm_uv: {
     size += sizeof(Real) * Model::spatial_dimension * 2;
     break;
   }
   case SynchronizationTag::_smm_res: /* FALLTHRU */
   case SynchronizationTag::_smm_mass: {
     size += sizeof(Real) * Model::spatial_dimension;
     break;
   }
   case SynchronizationTag::_for_dump: {
     size += sizeof(Real) * Model::spatial_dimension * 5;
     break;
   }
   default: {
     AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
   return size * dofs.size();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
                                    const Array<UInt> & dofs,
                                    const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case SynchronizationTag::_smm_uv: {
     packDOFDataHelper(*displacement, buffer, dofs);
     packDOFDataHelper(*velocity, buffer, dofs);
     break;
   }
   case SynchronizationTag::_smm_res: {
     packDOFDataHelper(*internal_force, buffer, dofs);
     break;
   }
   case SynchronizationTag::_smm_mass: {
     packDOFDataHelper(*mass, buffer, dofs);
     break;
   }
   case SynchronizationTag::_for_dump: {
     packDOFDataHelper(*displacement, buffer, dofs);
     packDOFDataHelper(*velocity, buffer, dofs);
     packDOFDataHelper(*acceleration, buffer, dofs);
     packDOFDataHelper(*internal_force, buffer, dofs);
     packDOFDataHelper(*external_force, buffer, dofs);
     break;
   }
   default: {
     AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
                                      const Array<UInt> & dofs,
                                      const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case SynchronizationTag::_smm_uv: {
     unpackDOFDataHelper(*displacement, buffer, dofs);
     unpackDOFDataHelper(*velocity, buffer, dofs);
     break;
   }
   case SynchronizationTag::_smm_res: {
     unpackDOFDataHelper(*internal_force, buffer, dofs);
     break;
   }
   case SynchronizationTag::_smm_mass: {
     unpackDOFDataHelper(*mass, buffer, dofs);
     break;
   }
   case SynchronizationTag::_for_dump: {
     unpackDOFDataHelper(*displacement, buffer, dofs);
     unpackDOFDataHelper(*velocity, buffer, dofs);
     unpackDOFDataHelper(*acceleration, buffer, dofs);
     unpackDOFDataHelper(*internal_force, buffer, dofs);
     unpackDOFDataHelper(*external_force, buffer, dofs);
     break;
   }
   default: {
     AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu