diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index c04cec9af..f8ad9811e 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -1,101 +1,100 @@
 
 #===============================================================================
 # @file   CMakeLists.txt
 #
 # @author Nicolas Richart <nicolas.richart@epfl.ch>
 #
 # @date creation: Fri Dec 12 2014
 # @date last modification: Mon Jan 18 2016
 #
 # @brief  CMake file for the python wrapping of akantu
 #
 # @section LICENSE
 #
 # Copyright (©) 2015 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/>.
 #
 #===============================================================================
 
 set(PYAKANTU_SRCS
   py_aka_common.cc
   py_aka_error.cc
   py_akantu.cc
   py_boundary_conditions.cc
   py_fe_engine.cc
   py_group_manager.cc
   py_mesh.cc
   py_model.cc
   py_parser.cc
   )
 
 package_is_activated(iohelper _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_dumpable.cc
     )
 endif()
 
 
 package_is_activated(solid_mechanics _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_solid_mechanics_model.cc
     py_material.cc
     )
 endif()
 
 package_is_activated(cohesive_element _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_solid_mechanics_model_cohesive.cc
     )
 endif()
 
 package_is_activated(heat_transfer _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_heat_transfer_model.cc
     )
 endif()
 
 package_is_activated(contact_mechanics _is_activated)
 if(_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_contact_mechanics_model.cc
     )
 endif()
 
 package_is_activated(model_couplers _is_activated)
 if(_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_model_couplers.cc
     )
 endif()
 
-
 pybind11_add_module(py11_akantu ${PYAKANTU_SRCS})
 target_link_libraries(py11_akantu PUBLIC akantu)
 set_target_properties(py11_akantu PROPERTIES DEBUG_POSTFIX "")
 
 set(_python_install_dir
   ${CMAKE_INSTALL_LIBDIR}/python${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}/site-packages)
 
 install(TARGETS py11_akantu
   LIBRARY DESTINATION ${_python_install_dir})
 
 install(DIRECTORY akantu
   DESTINATION ${_python_install_dir}
   FILES_MATCHING PATTERN "*.py")
diff --git a/python/py_material.cc b/python/py_material.cc
index e7418c522..a6f163006 100644
--- a/python/py_material.cc
+++ b/python/py_material.cc
@@ -1,257 +1,246 @@
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <material_selector.hh>
 #include <solid_mechanics_model.hh>
 #if defined(AKANTU_COHESIVE_ELEMENT)
 #include <solid_mechanics_model_cohesive.hh>
 #endif
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 #if not defined(PYBIND11_OVERRIDE)
 #define PYBIND11_OVERRIDE PYBIND11_OVERLOAD
 #define PYBIND11_OVERRIDE_PURE PYBIND11_OVERLOAD_PURE
 #endif
 
 namespace akantu {
 
 template <typename _Material> class PyMaterial : public _Material {
 
 public:
   /* Inherit the constructors */
   using _Material::_Material;
 
   ~PyMaterial() override = default;
 
   void initMaterial() override {
-      PYBIND11_OVERRIDE(void, _Material, initMaterial, ); // NOLINT
+    PYBIND11_OVERRIDE(void, _Material, initMaterial, ); // NOLINT
   };
   void computeStress(ElementType el_type,
                      GhostType ghost_type = _not_ghost) override {
     PYBIND11_OVERRIDE_PURE(void, _Material, computeStress, el_type, ghost_type);
   }
-  void computeTangentModuli(ElementType el_type,
-                            Array<Real> & tangent_matrix,
+  void computeTangentModuli(ElementType el_type, Array<Real> & tangent_matrix,
                             GhostType ghost_type = _not_ghost) override {
     PYBIND11_OVERRIDE(void, _Material, computeTangentModuli, el_type,
                       tangent_matrix, ghost_type);
   }
 
   void computePotentialEnergy(ElementType el_type) override {
     PYBIND11_OVERRIDE(void, _Material, computePotentialEnergy, el_type);
   }
 
   Real getPushWaveSpeed(const Element & element) const override {
     PYBIND11_OVERRIDE(Real, _Material, getPushWaveSpeed, element);
   }
 
   Real getShearWaveSpeed(const Element & element) const override {
     PYBIND11_OVERRIDE(Real, _Material, getShearWaveSpeed, element);
   }
 
   void registerInternal(const std::string & name, UInt nb_component) {
     this->internals[name] = std::make_shared<InternalField<Real>>(name, *this);
     AKANTU_DEBUG_INFO("alloc internal " << name << " "
                                         << &this->internals[name]);
 
     this->internals[name]->initialize(nb_component);
   }
 
   auto & getInternals() { return this->internals; }
 
 protected:
   std::map<std::string, std::shared_ptr<InternalField<Real>>> internals;
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void register_element_type_map_array(py::module & mod,
                                      const std::string & name) {
   py::class_<ElementTypeMapArray<T>, std::shared_ptr<ElementTypeMapArray<T>>>(
       mod, ("ElementTypeMapArray" + name).c_str())
       .def(
           "__call__",
           [](ElementTypeMapArray<T> & self, ElementType type,
              GhostType ghost_type) -> decltype(auto) {
             return self(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost,
           py::return_value_policy::reference)
       .def(
           "elementTypes",
           [](ElementTypeMapArray<T> & self, UInt _dim, GhostType _ghost_type,
              ElementKind _kind) -> decltype(auto) {
             auto types = self.elementTypes(_dim, _ghost_type, _kind);
             std::vector<ElementType> _types;
             for (auto && t : types) {
               _types.push_back(t);
             }
             return _types;
           },
           py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost,
           py::arg("kind") = _ek_regular);
 
   py::class_<InternalField<T>, ElementTypeMapArray<T>,
              std::shared_ptr<InternalField<T>>>(
       mod, ("InternalField" + name).c_str());
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename _Material>
 void define_material(py::module & mod, const std::string & name) {
   py::class_<_Material, PyMaterial<_Material>, Parsable>(
       mod, name.c_str(), py::multiple_inheritance())
       .def(py::init<SolidMechanicsModel &, const ID &>())
       .def(
           "getGradU",
           [](Material & self, ElementType el_type,
              GhostType ghost_type = _not_ghost) -> decltype(auto) {
             return self.getGradU(el_type, ghost_type);
           },
           py::arg("el_type"), py::arg("ghost_type") = _not_ghost,
           py::return_value_policy::reference)
       .def(
           "getStress",
           [](Material & self, ElementType el_type,
              GhostType ghost_type = _not_ghost) -> decltype(auto) {
             return self.getStress(el_type, ghost_type);
           },
           py::arg("el_type"), py::arg("ghost_type") = _not_ghost,
           py::return_value_policy::reference)
       .def(
           "getPotentialEnergy",
           [](Material & self, ElementType el_type) -> decltype(auto) {
             return self.getPotentialEnergy(el_type);
           },
           py::return_value_policy::reference)
       .def("initMaterial", &Material::initMaterial)
       .def("getModel", &Material::getModel)
       .def("registerInternal",
            [](Material & self, const std::string & name, UInt nb_component) {
              return dynamic_cast<PyMaterial<Material> &>(self).registerInternal(
                  name, nb_component);
            })
       .def(
           "getInternalFieldReal",
           [](Material & self, const ID & id, const ElementType & type,
              const GhostType & ghost_type) -> Array<Real> & {
             return self.getArray<Real>(id, type, ghost_type);
           },
           py::arg("id"), py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def(
           "getInternalFieldUInt",
           [](Material & self, const ID & id, const ElementType & type,
              const GhostType & ghost_type) -> Array<UInt> & {
             return self.getArray<UInt>(id, type, ghost_type);
           },
           py::arg("id"), py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def(
           "getElementFilter",
           [](Material & self, const ElementType & type,
              const GhostType & ghost_type) -> const Array<UInt> & {
             return self.getElementFilter()(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void register_material(py::module & mod) {
   py::class_<MaterialFactory>(mod, "MaterialFactory")
       .def_static(
           "getInstance",
           []() -> MaterialFactory & { return Material::getFactory(); },
           py::return_value_policy::reference)
       .def("registerAllocator",
            [](MaterialFactory & self, const std::string id, py::function func) {
              self.registerAllocator(
                  id,
                  [func, id](UInt dim, const ID & /*unused*/,
                             SolidMechanicsModel & model,
                             const ID & option) -> std::unique_ptr<Material> {
                    py::object obj = func(dim, id, model, option);
                    auto & ptr = py::cast<Material &>(obj);
 
                    obj.release();
                    return std::unique_ptr<Material>(&ptr);
                  });
            });
 
   register_element_type_map_array<Real>(mod, "Real");
   register_element_type_map_array<UInt>(mod, "UInt");
 
   define_material<Material>(mod, "Material");
+}
+
+/* -------------------------------------------------------------------------- */
+template <typename T>
+void register_data_material_selector(py::module & mod,
+                                     const std::string & name) {
+  py::class_<ElementDataMaterialSelector<T>, MaterialSelector,
+             std::shared_ptr<ElementDataMaterialSelector<T>>>(
+      mod, ("ElementDataMaterialSelector" + name).c_str());
+
+  py::class_<MeshDataMaterialSelector<T>, ElementDataMaterialSelector<T>,
+             std::shared_ptr<MeshDataMaterialSelector<T>>>(
+      mod, ("MeshDataMaterialSelector" + name).c_str())
+      .def(py::init<const std::string &, SolidMechanicsModel &, UInt>(),
+           py::arg("name"), py::arg("model"), py::arg("first_index") = 1);
+}
 
+/* -------------------------------------------------------------------------- */
+void register_material_selector(py::module & mod) {
   py::class_<MaterialSelector, std::shared_ptr<MaterialSelector>>(
       mod, "MaterialSelector")
       .def("setFallback",
            [](MaterialSelector & self, UInt f) { self.setFallback(f); })
       .def("setFallback",
            [](MaterialSelector & self,
               const std::shared_ptr<MaterialSelector> & fallback_selector) {
              self.setFallback(fallback_selector);
            })
       .def("setFallback",
            [](MaterialSelector & self, MaterialSelector & fallback_selector) {
              self.setFallback(fallback_selector);
            });
 
-  py::class_<MeshDataMaterialSelector<std::string>, MaterialSelector,
-             std::shared_ptr<MeshDataMaterialSelector<std::string>>>(
-      mod, "MeshDataMaterialSelectorString")
-      .def(py::init<const std::string &, const SolidMechanicsModel &, UInt>(),
-           py::arg("name"), py::arg("model"), py::arg("first_index") = 1);
-
 #if defined(AKANTU_COHESIVE_ELEMENT)
   py::class_<DefaultMaterialCohesiveSelector, MaterialSelector,
              std::shared_ptr<DefaultMaterialCohesiveSelector>>(
       mod, "DefaultMaterialCohesiveSelector")
       .def(py::init<const SolidMechanicsModelCohesive &>());
 
   py::class_<MeshDataMaterialCohesiveSelector, MaterialSelector,
              std::shared_ptr<MeshDataMaterialCohesiveSelector>>(
       mod, "MeshDataMaterialCohesiveSelector")
       .def(py::init<const SolidMechanicsModelCohesive &>());
 
   py::class_<MaterialCohesiveRulesSelector, MaterialSelector,
              std::shared_ptr<MaterialCohesiveRulesSelector>>(
       mod, "MaterialCohesiveRulesSelector")
       .def(py::init<const SolidMechanicsModelCohesive &,
                     const MaterialCohesiveRules &, const ID &>(),
            py::arg("model"), py::arg("rules"),
            py::arg("mesh_data_id") = "physical_names");
 #endif
-}
-
-/* -------------------------------------------------------------------------- */
-template <typename T>
-void register_data_material_selector(py::module & mod,
-                                          const std::string & name) {
-  py::class_<ElementDataMaterialSelector<T>, MaterialSelector,
-             std::shared_ptr<ElementDataMaterialSelector<T>>>(
-      mod, ("ElementDataMaterialSelector" + name).c_str());
-
-  py::class_<MeshDataMaterialSelector<T>, ElementDataMaterialSelector<T>,
-             std::shared_ptr<MeshDataMaterialSelector<T>>>(
-      mod, ("MeshDataMaterialSelector" + name).c_str())
-      .def(py::init<const std::string &, SolidMechanicsModel &, UInt>(),
-           py::arg("name"), py::arg("model"), py::arg("first_index") = 1);
-}
-
-/* -------------------------------------------------------------------------- */
-void register_material_selector(py::module & mod) {
-  py::class_<MaterialSelector, std::shared_ptr<MaterialSelector>>(
-      mod, "MaterialSelector");
 
   register_data_material_selector<std::string>(mod, "String");
 }
 
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu
-
diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh
index 468464f1e..bdbb0a45d 100644
--- a/src/common/aka_array_tmpl.hh
+++ b/src/common/aka_array_tmpl.hh
@@ -1,1362 +1,1364 @@
 /**
  * @file   aka_array_tmpl.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Jul 15 2010
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Inline functions of the classes Array<T> and ArrayBase
  *
  *
  * 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 /* Inline Functions Array<T>                                                  */
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh" // NOLINT
 #include "aka_static_memory.hh"
 /* -------------------------------------------------------------------------- */
 #include <memory>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_AKA_ARRAY_TMPL_HH_
 #define AKANTU_AKA_ARRAY_TMPL_HH_
 
 namespace akantu {
 
 namespace debug {
   struct ArrayException : public Exception {};
 } // namespace debug
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(UInt size,
                                                     UInt nb_component,
                                                     const ID & id)
     : ArrayBase(id) {
   allocate(size, nb_component);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(UInt size,
                                                     UInt nb_component,
                                                     const_reference value,
                                                     const ID & id)
     : ArrayBase(id) {
   allocate(size, nb_component, value);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(const ArrayDataLayer & vect,
                                                     const ID & id)
     : ArrayBase(vect, id) {
   this->data_storage = vect.data_storage;
   this->size_ = vect.size_;
   this->nb_component = vect.nb_component;
   this->values = this->data_storage.data();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 ArrayDataLayer<T, allocation_trait>::ArrayDataLayer(
     const std::vector<value_type> & vect) {
   this->data_storage = vect;
   this->size_ = vect.size();
   this->nb_component = 1;
   this->values = this->data_storage.data();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 ArrayDataLayer<T, allocation_trait> &
 ArrayDataLayer<T, allocation_trait>::operator=(const ArrayDataLayer & other) {
   if (this != &other) {
     this->data_storage = other.data_storage;
     this->nb_component = other.nb_component;
     this->size_ = other.size_;
     this->values = this->data_storage.data();
   }
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 void ArrayDataLayer<T, allocation_trait>::allocate(UInt new_size,
                                                    UInt nb_component) {
   this->nb_component = nb_component;
   this->resize(new_size);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 void ArrayDataLayer<T, allocation_trait>::allocate(UInt new_size,
                                                    UInt nb_component,
                                                    const T & val) {
   this->nb_component = nb_component;
   this->resize(new_size, val);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 void ArrayDataLayer<T, allocation_trait>::resize(UInt new_size) {
   this->data_storage.resize(new_size * this->nb_component);
   this->values = this->data_storage.data();
   this->size_ = new_size;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 void ArrayDataLayer<T, allocation_trait>::resize(UInt new_size,
                                                  const T & value) {
   this->data_storage.resize(new_size * this->nb_component, value);
   this->values = this->data_storage.data();
   this->size_ = new_size;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 void ArrayDataLayer<T, allocation_trait>::reserve(UInt size, UInt new_size) {
   if (new_size != UInt(-1)) {
     this->data_storage.resize(new_size * this->nb_component);
   }
 
   this->data_storage.reserve(size * this->nb_component);
   this->values = this->data_storage.data();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * append a tuple to the array with the value value for all components
  * @param value the new last tuple or the array will contain nb_component copies
  * of value
  */
 template <typename T, ArrayAllocationType allocation_trait>
 inline void ArrayDataLayer<T, allocation_trait>::push_back(const T & value) {
   this->data_storage.push_back(value);
   this->values = this->data_storage.data();
   this->size_ += 1;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * append a matrix or a vector to the array
  * @param new_elem a reference to a Matrix<T> or Vector<T> */
 template <typename T, ArrayAllocationType allocation_trait>
 template <template <typename> class C, typename>
 inline void
 ArrayDataLayer<T, allocation_trait>::push_back(const C<T> & new_elem) {
   AKANTU_DEBUG_ASSERT(
       nb_component == new_elem.size(),
       "The vector("
           << new_elem.size()
           << ") as not a size compatible with the Array (nb_component="
           << nb_component << ").");
   for (UInt i = 0; i < new_elem.size(); ++i) {
     this->data_storage.push_back(new_elem[i]);
   }
   this->values = this->data_storage.data();
   this->size_ += 1;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 inline UInt ArrayDataLayer<T, allocation_trait>::getAllocatedSize() const {
   return this->data_storage.capacity() / this->nb_component;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T, ArrayAllocationType allocation_trait>
 inline UInt ArrayDataLayer<T, allocation_trait>::getMemorySize() const {
   return this->data_storage.capacity() * sizeof(T);
 }
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 class ArrayDataLayer<T, ArrayAllocationType::_pod> : public ArrayBase {
 public:
   using value_type = T;
   using reference = value_type &;
   using pointer_type = value_type *;
   using const_reference = const value_type &;
 
 public:
   ~ArrayDataLayer() override { deallocate(); }
 
   /// Allocation of a new vector
   ArrayDataLayer(UInt size = 0, UInt nb_component = 1, const ID & id = "")
       : ArrayBase(id) {
     allocate(size, nb_component);
   }
 
   /// Allocation of a new vector with a default value
   ArrayDataLayer(UInt size, UInt nb_component, const_reference value,
                  const ID & id = "")
       : ArrayBase(id) {
     allocate(size, nb_component, value);
   }
 
   /// Copy constructor (deep copy)
   ArrayDataLayer(const ArrayDataLayer & vect, const ID & id = "")
       : ArrayBase(vect, id) {
     allocate(vect.size(), vect.getNbComponent());
     std::copy_n(vect.storage(), this->size_ * this->nb_component, values);
   }
 
   /// Copy constructor (deep copy)
   explicit ArrayDataLayer(const std::vector<value_type> & vect) {
     allocate(vect.size(), 1);
     std::copy_n(vect.data(), this->size_ * this->nb_component, values);
   }
 
   // copy operator
   inline ArrayDataLayer & operator=(const ArrayDataLayer & other) {
     if (this != &other) {
       allocate(other.size(), other.getNbComponent());
       std::copy_n(other.storage(), this->size_ * this->nb_component, values);
     }
     return *this;
   }
 
   // move constructor
   inline ArrayDataLayer(ArrayDataLayer && other) noexcept = default;
 
   // move assign
   inline ArrayDataLayer & operator=(ArrayDataLayer && other) noexcept = default;
 
 protected:
   // deallocate the memory
   virtual void deallocate() {
     // NOLINTNEXTLINE(cppcoreguidelines-owning-memory,
     // cppcoreguidelines-no-malloc)
     free(this->values);
   }
 
   // allocate the memory
   virtual inline void allocate(UInt size, UInt nb_component) {
     if (size != 0) { // malloc can return a non NULL pointer in case size is 0
       this->values = static_cast<T *>(                   // NOLINT
           std::malloc(nb_component * size * sizeof(T))); // NOLINT
     }
 
     if (this->values == nullptr and size != 0) {
       throw std::bad_alloc();
     }
     this->nb_component = nb_component;
     this->allocated_size = this->size_ = size;
   }
 
   // allocate and initialize the memory
   virtual inline void allocate(UInt size, UInt nb_component, const T & value) {
     allocate(size, nb_component);
     std::fill_n(values, size * nb_component, value);
   }
 
 public:
   /// append a tuple of size nb_component containing value
   inline void push_back(const_reference value) {
     resize(this->size_ + 1, value);
   }
 
   /// append a Vector or a Matrix
   template <template <typename> class C,
             typename = std::enable_if_t<aka::is_tensor<C<T>>::value or
                                         aka::is_tensor_proxy<C<T>>::value>>
   inline void push_back(const C<T> & new_elem) {
     AKANTU_DEBUG_ASSERT(
         nb_component == new_elem.size(),
         "The vector("
             << new_elem.size()
             << ") as not a size compatible with the Array (nb_component="
             << nb_component << ").");
     this->resize(this->size_ + 1);
     std::copy_n(new_elem.storage(), new_elem.size(),
                 values + this->nb_component * (this->size_ - 1));
   }
 
   /// changes the allocated size but not the size
   virtual void reserve(UInt size, UInt new_size = UInt(-1)) {
     UInt tmp_size = this->size_;
     if (new_size != UInt(-1)) {
       tmp_size = new_size;
     }
     this->resize(size);
     this->size_ = std::min(this->size_, tmp_size);
   }
 
   /// change the size of the Array
   virtual void resize(UInt size) {
     if (size * this->nb_component == 0) {
       free(values); // NOLINT: cppcoreguidelines-no-malloc
       values = nullptr;
       this->allocated_size = 0;
     } else {
       if (this->values == nullptr) {
         this->allocate(size, this->nb_component);
         return;
       }
 
       Int diff = size - allocated_size;
       UInt size_to_allocate = (std::abs(diff) > AKANTU_MIN_ALLOCATION)
                                   ? size
                                   : (diff > 0)
                                         ? allocated_size + AKANTU_MIN_ALLOCATION
                                         : allocated_size;
 
       if (size_to_allocate ==
           allocated_size) { // otherwhy the reserve + push_back might fail...
         this->size_ = size;
         return;
       }
 
       auto * tmp_ptr = reinterpret_cast<T *>( // NOLINT
           realloc(this->values,
                   size_to_allocate * this->nb_component * sizeof(T)));
 
       if (tmp_ptr == nullptr) {
         StaticMemory::getStaticMemory().printself(std::cerr);
         
         throw std::bad_alloc();
       }
 
+      this->values = tmp_ptr;
+      this->allocated_size = size_to_allocate;
     }
 
     this->size_ = size;
   }
 
   /// change the size of the Array and initialize the values
   virtual void resize(UInt size, const T & val) {
     UInt tmp_size = this->size_;
     this->resize(size);
     if (size > tmp_size) {
       // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
       std::fill_n(values + this->nb_component * tmp_size,
                   (size - tmp_size) * this->nb_component, val);
     }
   }
 
   /// get the amount of space allocated in bytes
   inline UInt getMemorySize() const final {
     return this->allocated_size * this->nb_component * sizeof(T);
   }
 
   /// Get the real size allocated in memory
   inline UInt getAllocatedSize() const { return this->allocated_size; }
 
   /// give the address of the memory allocated for this vector
   T * storage() const { return values; };
 
 protected:
   /// allocation type agnostic  data access
   T * values{nullptr};
 
   UInt allocated_size{0};
 };
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 inline auto Array<T, is_scal>::operator()(UInt i, UInt j) -> reference {
   AKANTU_DEBUG_ASSERT(this->size_ > 0,
                       "The array \"" << this->id << "\" is empty");
   AKANTU_DEBUG_ASSERT((i < this->size_) && (j < this->nb_component),
                       "The value at position ["
                           << i << "," << j << "] is out of range in array \""
                           << this->id << "\"");
   return this->values[i * this->nb_component + j];
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 inline auto Array<T, is_scal>::operator()(UInt i, UInt j) const
     -> const_reference {
   AKANTU_DEBUG_ASSERT(this->size_ > 0,
                       "The array \"" << this->id << "\" is empty");
   AKANTU_DEBUG_ASSERT((i < this->size_) && (j < this->nb_component),
                       "The value at position ["
                           << i << "," << j << "] is out of range in array \""
                           << this->id << "\"");
   // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
   return this->values[i * this->nb_component + j];
 }
 
 template <class T, bool is_scal>
 inline auto Array<T, is_scal>::operator[](UInt i) -> reference {
   AKANTU_DEBUG_ASSERT(this->size_ > 0,
                       "The array \"" << this->id << "\" is empty");
   AKANTU_DEBUG_ASSERT((i < this->size_ * this->nb_component),
                       "The value at position ["
                           << i << "] is out of range in array \"" << this->id
                           << "\"");
   return this->values[i];
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 inline auto Array<T, is_scal>::operator[](UInt i) const -> const_reference {
   AKANTU_DEBUG_ASSERT(this->size_ > 0,
                       "The array \"" << this->id << "\" is empty");
   AKANTU_DEBUG_ASSERT((i < this->size_ * this->nb_component),
                       "The value at position ["
                           << i << "] is out of range in array \"" << this->id
                           << "\"");
   return this->values[i];
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * erase an element. If the erased element is not the last of the array, the
  * last element is moved into the hole in order to maintain contiguity. This
  * may invalidate existing iterators (For instance an iterator obtained by
  * Array::end() is no longer correct) and will change the order of the
  * elements.
  * @param i index of element to erase
  */
 template <class T, bool is_scal> inline void Array<T, is_scal>::erase(UInt i) {
   AKANTU_DEBUG_IN();
   AKANTU_DEBUG_ASSERT((this->size_ > 0), "The array is empty");
   AKANTU_DEBUG_ASSERT((i < this->size_), "The element at position ["
                                              << i << "] is out of range (" << i
                                              << ">=" << this->size_ << ")");
 
   if (i != (this->size_ - 1)) {
     for (UInt j = 0; j < this->nb_component; ++j) {
       // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
       this->values[i * this->nb_component + j] =
           this->values[(this->size_ - 1) * this->nb_component + j];
     }
   }
 
   this->resize(this->size_ - 1);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Subtract another array entry by entry from this array in place. Both arrays
  * must
  * have the same size and nb_component. If the arrays have different shapes,
  * code compiled in debug mode will throw an expeption and optimised code
  * will behave in an unpredicted manner
  * @param other array to subtract from this
  * @return reference to modified this
  */
 template <class T, bool is_scal>
 Array<T, is_scal> &
 Array<T, is_scal>::operator-=(const Array<T, is_scal> & other) {
   AKANTU_DEBUG_ASSERT((this->size_ == other.size_) &&
                           (this->nb_component == other.nb_component),
                       "The too array don't have the same sizes");
 
   T * a = this->values;
   T * b = other.storage();
   for (UInt i = 0; i < this->size_ * this->nb_component; ++i) {
     *a -= *b;
     ++a;
     ++b;
   }
 
   return *this;
 }
 
 /* --------------------------------------------------------------------------
  */
 /**
  * Add another array entry by entry to this array in
  * place. Both arrays must have the same size and
  * nb_component. If the arrays have different shapes, code
  * compiled in debug mode will throw an expeption and
  * optimised code will behave in an unpredicted manner
  * @param other array to add to this
  * @return reference to modified this
  */
 template <class T, bool is_scal>
 Array<T, is_scal> &
 Array<T, is_scal>::operator+=(const Array<T, is_scal> & other) {
   AKANTU_DEBUG_ASSERT((this->size_ == other.size()) &&
                           (this->nb_component == other.nb_component),
                       "The too array don't have the same sizes");
 
   T * a = this->values;
   T * b = other.storage();
   for (UInt i = 0; i < this->size_ * this->nb_component; ++i) {
     *a++ += *b++;
   }
 
   return *this;
 }
 
 /* --------------------------------------------------------------------------
  */
 /**
  * Multiply all entries of this array by a scalar in place
  * @param alpha scalar multiplicant
  * @return reference to modified this
  */
 
 template <class T, bool is_scal>
 Array<T, is_scal> & Array<T, is_scal>::operator*=(const T & alpha) {
   T * a = this->values;
   for (UInt i = 0; i < this->size_ * this->nb_component; ++i) {
     *a++ *= alpha;
   }
 
   return *this;
 }
 
 /* --------------------------------------------------------------------------
  */
 /**
  * Compare this array element by element to another.
  * @param other array to compare to
  * @return true it all element are equal and arrays have
  * the same shape, else false
  */
 template <class T, bool is_scal>
 bool Array<T, is_scal>::operator==(const Array<T, is_scal> & other) const {
   bool equal = this->nb_component == other.nb_component &&
                this->size_ == other.size_ && this->id == other.id;
   if (not equal) {
     return false;
   }
   if (this->values == other.storage()) {
     return true;
   }
 
   // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
   return std::equal(this->values,
                     this->values + this->size_ * this->nb_component,
                     other.storage());
 }
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 bool Array<T, is_scal>::operator!=(const Array<T, is_scal> & other) const {
   return !operator==(other);
 }
 
 /* --------------------------------------------------------------------------
  */
 /**
  * set all tuples of the array to a given vector or matrix
  * @param vm Matrix or Vector to fill the array with
  */
 template <class T, bool is_scal>
 template <template <typename> class C, typename>
 inline void Array<T, is_scal>::set(const C<T> & vm) {
   AKANTU_DEBUG_ASSERT(this->nb_component == vm.size(),
                       "The size of the object does not "
                       "match the number of components");
   for (T * it = this->values;
        it < this->values + this->nb_component * this->size_;
        it += this->nb_component) {
     std::copy_n(vm.storage(), this->nb_component, it);
   }
 }
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 void Array<T, is_scal>::append(const Array<T> & other) {
   AKANTU_DEBUG_ASSERT(this->nb_component == other.nb_component,
                       "Cannot append an array with a "
                       "different number of component");
   UInt old_size = this->size_;
   this->resize(this->size_ + other.size());
 
   T * tmp = this->values + this->nb_component * old_size;
   std::copy_n(other.storage(), other.size() * this->nb_component, tmp);
 }
 
 /* --------------------------------------------------------------------------
  */
 /* Functions Array<T, is_scal> */
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(UInt size, UInt nb_component, const ID & id)
     : parent(size, nb_component, id) {}
 
 template <>
 inline Array<std::string, false>::Array(UInt size, UInt nb_component,
                                         const ID & id)
     : parent(size, nb_component, "", id) {}
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(UInt size, UInt nb_component, const_reference value,
                          const ID & id)
     : parent(size, nb_component, value, id) {}
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(const Array & vect, const ID & id)
     : parent(vect, id) {}
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 Array<T, is_scal> &
 Array<T, is_scal>::operator=(const Array<T, is_scal> & other) {
   AKANTU_DEBUG_WARNING("You are copying the array "
                        << this->id << " are you sure it is on purpose");
 
   if (&other == this) {
     return *this;
   }
 
   parent::operator=(other);
   return *this;
 }
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(const std::vector<T> & vect) : parent(vect) {}
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal> Array<T, is_scal>::~Array() = default;
 
 /* --------------------------------------------------------------------------
  */
 /**
  * search elem in the array, return  the position of the
  * first occurrence or -1 if not found
  *  @param elem the element to look for
  *  @return index of the first occurrence of elem or -1 if
  * elem is not present
  */
 template <class T, bool is_scal>
 UInt Array<T, is_scal>::find(const_reference elem) const {
   AKANTU_DEBUG_IN();
 
   auto begin = this->begin();
   auto end = this->end();
   auto it = std::find(begin, end, elem);
 
   AKANTU_DEBUG_OUT();
   return (it != end) ? it - begin : UInt(-1);
 }
 
 /* --------------------------------------------------------------------------
  */
 // template <class T, bool is_scal> UInt Array<T,
 // is_scal>::find(T elem[]) const
 // {
 //   AKANTU_DEBUG_IN();
 //   T * it = this->values;
 //   UInt i = 0;
 //   for (; i < this->size_; ++i) {
 //     if (*it == elem[0]) {
 //       T * cit = it;
 //       UInt c = 0;
 //       for (; (c < this->nb_component) && (*cit ==
 //       elem[c]); ++c, ++cit)
 //         ;
 //       if (c == this->nb_component) {
 //         AKANTU_DEBUG_OUT();
 //         return i;
 //       }
 //     }
 //     it += this->nb_component;
 //   }
 //   return UInt(-1);
 // }
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 template <template <typename> class C, typename>
 inline UInt Array<T, is_scal>::find(const C<T> & elem) {
   AKANTU_DEBUG_ASSERT(elem.size() == this->nb_component,
                       "Cannot find an element with a wrong size ("
                           << elem.size() << ") != " << this->nb_component);
   return this->find(*elem.storage());
 }
 
 /* --------------------------------------------------------------------------
  */
 /**
  * copy the content of another array. This overwrites the
  * current content.
  * @param other Array to copy into this array. It has to
  * have the same nb_component as this. If compiled in
  * debug mode, an incorrect other will result in an
  * exception being thrown. Optimised code may result in
  * unpredicted behaviour.
  * @param no_sanity_check turns off all checkes
  */
 template <class T, bool is_scal>
 void Array<T, is_scal>::copy(const Array<T, is_scal> & other,
                              bool no_sanity_check) {
   AKANTU_DEBUG_IN();
 
   if (not no_sanity_check and (other.nb_component != this->nb_component)) {
     AKANTU_ERROR("The two arrays do not have the same "
                  "number of components");
   }
 
   this->resize((other.size_ * other.nb_component) / this->nb_component);
 
   std::copy_n(other.storage(), this->size_ * this->nb_component, this->values);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* --------------------------------------------------------------------------
  */
 template <bool is_scal> class ArrayPrintHelper {
 public:
   template <typename T>
   static void print_content(const Array<T> & vect, std::ostream & stream,
                             int indent) {
     std::string space(indent, AKANTU_INDENT);
 
     stream << space << " + values         : {";
     for (UInt i = 0; i < vect.size(); ++i) {
       stream << "{";
       for (UInt j = 0; j < vect.getNbComponent(); ++j) {
         stream << vect(i, j);
         if (j != vect.getNbComponent() - 1) {
           stream << ", ";
         }
       }
       stream << "}";
       if (i != vect.size() - 1) {
         stream << ", ";
       }
     }
     stream << "}" << std::endl;
   }
 };
 
 template <> class ArrayPrintHelper<false> {
 public:
   template <typename T>
   static void print_content(__attribute__((unused)) const Array<T> & vect,
                             __attribute__((unused)) std::ostream & stream,
                             __attribute__((unused)) int indent) {}
 };
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 void Array<T, is_scal>::printself(std::ostream & stream, int indent) const {
   std::string space(indent, AKANTU_INDENT);
 
   std::streamsize prec = stream.precision();
   std::ios_base::fmtflags ff = stream.flags();
 
   stream.setf(std::ios_base::showbase);
   stream.precision(2);
 
   stream << space << "Array<" << debug::demangle(typeid(T).name()) << "> ["
          << std::endl;
   stream << space << " + id             : " << this->id << std::endl;
   stream << space << " + size           : " << this->size_ << std::endl;
   stream << space << " + nb_component   : " << this->nb_component << std::endl;
   stream << space << " + allocated size : " << this->getAllocatedSize()
          << std::endl;
   stream << space
          << " + memory size    : " << printMemorySize<T>(this->getMemorySize())
          << std::endl;
   if (not AKANTU_DEBUG_LEVEL_IS_TEST()) {
     stream << space << " + address        : " << std::hex << this->values
            << std::dec << std::endl;
   }
 
   stream.precision(prec);
   stream.flags(ff);
 
   if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) {
     ArrayPrintHelper<is_scal or std::is_enum<T>::value>::print_content(
         *this, stream, indent);
   }
 
   stream << space << "]" << std::endl;
 }
 
 /* --------------------------------------------------------------------------
  */
 /* Inline Functions ArrayBase */
 /* --------------------------------------------------------------------------
  */
 
 // inline bool ArrayBase::empty() { return (this->size_ ==
 // 0); }
 
 /* --------------------------------------------------------------------------
  */
 /* Iterators */
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 template <class R, class daughter, class IR, bool is_tensor>
 class Array<T, is_scal>::iterator_internal {
 public:
   using value_type = R;
   using pointer = R *;
   using reference = R &;
   using const_reference = const R &;
   using internal_value_type = IR;
   using internal_pointer = IR *;
   using difference_type = std::ptrdiff_t;
   using iterator_category = std::random_access_iterator_tag;
   static_assert(not is_tensor, "Cannot handle tensors");
 
 public:
   iterator_internal(pointer data = nullptr) : ret(data), initial(data){};
   iterator_internal(const iterator_internal & it) = default;
   iterator_internal(iterator_internal && it) noexcept = default;
 
   virtual ~iterator_internal() = default;
 
   inline iterator_internal & operator=(const iterator_internal & it) = default;
   inline iterator_internal &
   operator=(iterator_internal && it) noexcept = default;
 
   UInt getCurrentIndex() { return (this->ret - this->initial); };
 
   inline reference operator*() { return *ret; };
   inline const_reference operator*() const { return *ret; };
   inline pointer operator->() { return ret; };
   inline daughter & operator++() {
     ++ret;
     return static_cast<daughter &>(*this);
   };
   inline daughter & operator--() {
     --ret;
     return static_cast<daughter &>(*this);
   };
 
   inline daughter & operator+=(const UInt n) {
     ret += n;
     return static_cast<daughter &>(*this);
   }
   inline daughter & operator-=(const UInt n) {
     ret -= n;
     return static_cast<daughter &>(*this);
   }
 
   inline reference operator[](const UInt n) { return ret[n]; }
 
   inline bool operator==(const iterator_internal & other) const {
     return ret == other.ret;
   }
   inline bool operator!=(const iterator_internal & other) const {
     return ret != other.ret;
   }
   inline bool operator<(const iterator_internal & other) const {
     return ret < other.ret;
   }
   inline bool operator<=(const iterator_internal & other) const {
     return ret <= other.ret;
   }
   inline bool operator>(const iterator_internal & other) const {
     return ret > other.ret;
   }
   inline bool operator>=(const iterator_internal & other) const {
     return ret >= other.ret;
   }
 
   inline daughter operator-(difference_type n) { return daughter(ret - n); }
   inline daughter operator+(difference_type n) { return daughter(ret + n); }
 
   inline difference_type operator-(const iterator_internal & b) {
     return ret - b.ret;
   }
 
   inline pointer data() const { return ret; }
 
 protected:
   pointer ret{nullptr};
   pointer initial{nullptr};
 };
 
 /* --------------------------------------------------------------------------
  */
 /**
  * Specialization for scalar types
  */
 template <class T, bool is_scal>
 template <class R, class daughter, class IR>
 class Array<T, is_scal>::iterator_internal<R, daughter, IR, true> {
 public:
   using value_type = R;
   using pointer = R *;
-  using pointer_type = typename Array<T, is_scal>::pointer_type;
   using reference = R &;
   using proxy = typename R::proxy;
   using const_proxy = const typename R::proxy;
   using const_reference = const R &;
   using internal_value_type = IR;
   using internal_pointer = IR *;
   using difference_type = std::ptrdiff_t;
   using iterator_category = std::random_access_iterator_tag;
+  using pointer_type = typename Array<T, is_scal>::pointer_type;
 
 public:
   iterator_internal() = default;
 
   iterator_internal(pointer_type data, UInt _offset)
       : _offset(_offset), initial(data), ret(nullptr), ret_ptr(data) {
     AKANTU_ERROR("The constructor should never be called "
                  "it is just an ugly trick...");
   }
 
   iterator_internal(std::unique_ptr<internal_value_type> && wrapped)
       : _offset(wrapped->size()), initial(wrapped->storage()),
         ret(std::move(wrapped)), ret_ptr(ret->storage()) {}
 
   iterator_internal(const iterator_internal & it) {
     if (this != &it) {
       this->_offset = it._offset;
       this->initial = it.initial;
       this->ret_ptr = it.ret_ptr;
       this->ret = std::make_unique<internal_value_type>(*it.ret, false);
     }
   }
 
   iterator_internal(iterator_internal && it) noexcept = default;
 
   virtual ~iterator_internal() = default;
 
   inline iterator_internal & operator=(const iterator_internal & it) {
     if (this != &it) {
       this->_offset = it._offset;
       this->initial = it.initial;
       this->ret_ptr = it.ret_ptr;
       if (this->ret) {
         this->ret->shallowCopy(*it.ret);
       } else {
         this->ret = std::make_unique<internal_value_type>(*it.ret, false);
       }
     }
     return *this;
   }
 
   inline iterator_internal &
   operator=(iterator_internal && it) noexcept = default;
 
   UInt getCurrentIndex() {
     return (this->ret_ptr - this->initial) / this->_offset;
   };
 
   inline reference operator*() {
     ret->values = ret_ptr;
     return *ret;
   };
   inline const_reference operator*() const {
     ret->values = ret_ptr;
     return *ret;
   };
   inline pointer operator->() {
     ret->values = ret_ptr;
     return ret.get();
   };
   inline daughter & operator++() {
     ret_ptr += _offset;
     return static_cast<daughter &>(*this);
   };
   inline daughter & operator--() {
     ret_ptr -= _offset;
     return static_cast<daughter &>(*this);
   };
 
   inline daughter & operator+=(const UInt n) {
     ret_ptr += _offset * n;
     return static_cast<daughter &>(*this);
   }
   inline daughter & operator-=(const UInt n) {
     ret_ptr -= _offset * n;
     return static_cast<daughter &>(*this);
   }
 
   inline proxy operator[](const UInt n) {
     ret->values = ret_ptr + n * _offset;
     return proxy(*ret);
   }
   inline const_proxy operator[](const UInt n) const { // NOLINT
     ret->values = ret_ptr + n * _offset;
     return const_proxy(*ret);
   }
 
   inline bool operator==(const iterator_internal & other) const {
     return this->ret_ptr == other.ret_ptr;
   }
   inline bool operator!=(const iterator_internal & other) const {
     return this->ret_ptr != other.ret_ptr;
   }
   inline bool operator<(const iterator_internal & other) const {
     return this->ret_ptr < other.ret_ptr;
   }
   inline bool operator<=(const iterator_internal & other) const {
     return this->ret_ptr <= other.ret_ptr;
   }
   inline bool operator>(const iterator_internal & other) const {
     return this->ret_ptr > other.ret_ptr;
   }
   inline bool operator>=(const iterator_internal & other) const {
     return this->ret_ptr >= other.ret_ptr;
   }
 
   inline daughter operator+(difference_type n) {
     daughter tmp(static_cast<daughter &>(*this));
     tmp += n;
     return tmp;
   }
   inline daughter operator-(difference_type n) {
     daughter tmp(static_cast<daughter &>(*this));
     tmp -= n;
     return tmp;
   }
 
   inline difference_type operator-(const iterator_internal & b) {
     return (this->ret_ptr - b.ret_ptr) / _offset;
   }
 
   inline pointer_type data() const { return ret_ptr; }
   inline difference_type offset() const { return _offset; }
 
 protected:
   UInt _offset{0};
   pointer_type initial{nullptr};
   std::unique_ptr<internal_value_type> ret{nullptr};
   pointer_type ret_ptr{nullptr};
 };
 
 /* -------------------------------------------------------------------------- */
 /* Iterators */
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 template <typename R>
 class Array<T, is_scal>::const_iterator
     : public iterator_internal<const R, Array<T, is_scal>::const_iterator<R>,
                                R> {
 public:
   using parent = iterator_internal<const R, const_iterator, R>;
   using value_type = typename parent::value_type;
   using pointer = typename parent::pointer;
   using reference = typename parent::reference;
   using difference_type = typename parent::difference_type;
   using iterator_category = typename parent::iterator_category;
 
 public:
   ~const_iterator() override = default;
 
   const_iterator() = default;
   const_iterator(const const_iterator & it) = default;
   const_iterator(const_iterator && it) noexcept = default;
 
   const_iterator & operator=(const const_iterator & it) = default;
   const_iterator & operator=(const_iterator && it) noexcept = default;
 
   template <typename P,
             typename = std::enable_if_t<not aka::is_tensor<P>::value>>
   const_iterator(P * data) : parent(data) {}
 
   template <typename UP_P, typename = std::enable_if_t<aka::is_tensor<
                                typename UP_P::element_type>::value>>
   const_iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
 };
 
 /* -------------------------------------------------------------------------- */
 template <class T, class R, bool is_tensor_ = aka::is_tensor<R>::value>
 struct ConstConverterIteratorHelper {
   using const_iterator = typename Array<T>::template const_iterator<R>;
   using iterator = typename Array<T>::template iterator<R>;
 
   static inline const_iterator convert(const iterator & it) {
     return const_iterator(std::unique_ptr<R>(new R(*it, false)));
   }
 };
 
 template <class T, class R> struct ConstConverterIteratorHelper<T, R, false> {
   using const_iterator = typename Array<T>::template const_iterator<R>;
   using iterator = typename Array<T>::template iterator<R>;
   static inline const_iterator convert(const iterator & it) {
     return const_iterator(it.data());
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 template <typename R>
 class Array<T, is_scal>::iterator
     : public iterator_internal<R, Array<T, is_scal>::iterator<R>> {
 public:
   using parent = iterator_internal<R, iterator>;
   using value_type = typename parent::value_type;
   using pointer = typename parent::pointer;
   using reference = typename parent::reference;
   using difference_type = typename parent::difference_type;
   using iterator_category = typename parent::iterator_category;
 
 public:
   ~iterator() override = default;
   iterator() = default;
   iterator(const iterator & it) = default;
   iterator(iterator && it) noexcept = default;
   iterator & operator=(const iterator & it) = default;
   iterator & operator=(iterator && it) noexcept = default;
 
   template <typename P,
             typename = std::enable_if_t<not aka::is_tensor<P>::value>>
   iterator(P * data) : parent(data) {}
 
   template <typename UP_P, typename = std::enable_if_t<aka::is_tensor<
                                typename UP_P::element_type>::value>>
   iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
 
   operator const_iterator<R>() {
     return ConstConverterIteratorHelper<T, R>::convert(*this);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 /* Begin/End functions implementation */
 /* -------------------------------------------------------------------------- */
 namespace detail {
   template <class Tuple, size_t... Is>
   constexpr auto take_front_impl(Tuple && t,
                                  std::index_sequence<Is...> /*idxs*/) {
     return std::make_tuple(std::get<Is>(std::forward<Tuple>(t))...);
   }
 
   template <size_t N, class Tuple> constexpr auto take_front(Tuple && t) {
     return take_front_impl(std::forward<Tuple>(t),
                            std::make_index_sequence<N>{});
   }
 
   template <typename... V> constexpr auto product_all(V &&... v) {
     std::common_type_t<int, V...> result = 1;
     (void)std::initializer_list<int>{(result *= v, 0)...};
     return result;
   }
 
   template <typename... T> std::string to_string_all(T &&... t) {
     if (sizeof...(T) == 0) {
       return "";
     }
 
     std::stringstream ss;
     bool noComma = true;
     ss << "(";
     (void)std::initializer_list<bool>{
         (ss << (noComma ? "" : ", ") << t, noComma = false)...};
     ss << ")";
     return ss.str();
   }
 
   template <std::size_t N> struct InstantiationHelper {
     template <typename type, typename T, typename... Ns>
     static auto instantiate(T && data, Ns... ns) {
       return std::make_unique<type>(data, ns...);
     }
   };
 
   template <> struct InstantiationHelper<0> {
     template <typename type, typename T> static auto instantiate(T && data) {
       return data;
     }
   };
 
   template <typename Arr, typename T, typename... Ns>
   decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) {
     using type = IteratorHelper_t<sizeof...(Ns) - 1, T>;
     using array_type = std::decay_t<Arr>;
     using iterator =
         std::conditional_t<std::is_const<std::remove_reference_t<Arr>>::value,
                            typename array_type::template const_iterator<type>,
                            typename array_type::template iterator<type>>;
     static_assert(sizeof...(Ns), "You should provide a least one size");
 
     if (array.getNbComponent() * array.size() !=
         product_all(std::forward<Ns>(ns)...)) {
       AKANTU_CUSTOM_EXCEPTION_INFO(
           debug::ArrayException(),
           "The iterator on "
               << debug::demangle(typeid(Arr).name())
               << to_string_all(array.size(), array.getNbComponent())
               << "is not compatible with the type "
               << debug::demangle(typeid(type).name()) << to_string_all(ns...));
     }
 
     auto && wrapped = aka::apply(
         [&](auto... n) {
           return InstantiationHelper<sizeof...(n)>::template instantiate<type>(
               data, n...);
         },
         take_front<sizeof...(Ns) - 1>(std::make_tuple(ns...)));
 
     return iterator(std::move(wrapped));
   }
 } // namespace detail
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) {
   return detail::get_iterator(*this, this->values, std::forward<Ns>(ns)...,
                               this->size_);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) {
   return detail::get_iterator(*this,
                               this->values + this->nb_component * this->size_,
                               std::forward<Ns>(ns)..., this->size_);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) const {
   return detail::get_iterator(*this, this->values, std::forward<Ns>(ns)...,
                               this->size_);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) const {
   return detail::get_iterator(*this,
                               this->values + this->nb_component * this->size_,
                               std::forward<Ns>(ns)..., this->size_);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) {
   return detail::get_iterator(*this, this->values, std::forward<Ns>(ns)...);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) {
   return detail::get_iterator(
       *this, this->values + detail::product_all(std::forward<Ns>(ns)...),
       std::forward<Ns>(ns)...);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) const {
   return detail::get_iterator(*this, this->values, std::forward<Ns>(ns)...);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) const {
   return detail::get_iterator(
       *this, this->values + detail::product_all(std::forward<Ns>(ns)...),
       std::forward<Ns>(ns)...);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Views */
 /* -------------------------------------------------------------------------- */
 namespace detail {
   template <typename Array, typename... Ns> class ArrayView {
     using tuple = std::tuple<Ns...>;
 
   public:
     ~ArrayView() = default;
     ArrayView(Array && array, Ns... ns) noexcept
         : array(array), sizes(std::move(ns)...) {}
 
     ArrayView(const ArrayView & array_view) = default;
     ArrayView & operator=(const ArrayView & array_view) = default;
 
     ArrayView(ArrayView && array_view) noexcept = default;
     ArrayView & operator=(ArrayView && array_view) noexcept = default;
 
     decltype(auto) begin() {
       return aka::apply(
           [&](auto &&... ns) { return array.get().begin_reinterpret(ns...); },
           sizes);
     }
 
     decltype(auto) begin() const {
       return aka::apply(
           [&](auto &&... ns) { return array.get().begin_reinterpret(ns...); },
           sizes);
     }
 
     decltype(auto) end() {
       return aka::apply(
           [&](auto &&... ns) { return array.get().end_reinterpret(ns...); },
           sizes);
     }
 
     decltype(auto) end() const {
       return aka::apply(
           [&](auto &&... ns) { return array.get().end_reinterpret(ns...); },
           sizes);
     }
 
     decltype(auto) size() const {
       return std::get<std::tuple_size<tuple>::value - 1>(sizes);
     }
 
     decltype(auto) dims() const { return std::tuple_size<tuple>::value - 1; }
 
   private:
     std::reference_wrapper<std::remove_reference_t<Array>> array;
     tuple sizes;
   };
 } // namespace detail
 
 /* -------------------------------------------------------------------------- */
 template <typename Array, typename... Ns>
 decltype(auto) make_view(Array && array, const Ns... ns) {
   static_assert(aka::conjunction<std::is_integral<std::decay_t<Ns>>...>::value,
                 "Ns should be integral types");
   AKANTU_DEBUG_ASSERT((detail::product_all(ns...) != 0),
                       "You must specify non zero dimensions");
   auto size = std::forward<decltype(array)>(array).size() *
               std::forward<decltype(array)>(array).getNbComponent() /
               detail::product_all(ns...);
 
   return detail::ArrayView<Array, std::common_type_t<size_t, Ns>...,
                            std::common_type_t<size_t, decltype(size)>>(
       std::forward<Array>(array), std::move(ns)..., size);
 }
 
 /* --------------------------------------------------------------------------
  */
 template <class T, bool is_scal>
 template <typename R>
 inline typename Array<T, is_scal>::template iterator<R>
 Array<T, is_scal>::erase(const iterator<R> & it) {
   T * curr = it.data();
   UInt pos = (curr - this->values) / this->nb_component;
   erase(pos);
   iterator<R> rit = it;
   return --rit;
 }
 
 } // namespace akantu
 
 #endif /* AKANTU_AKA_ARRAY_TMPL_HH_ */
diff --git a/src/mesh_utils/cohesive_element_inserter_helper.cc b/src/mesh_utils/cohesive_element_inserter_helper.cc
index b50fea3d3..22f2777fa 100644
--- a/src/mesh_utils/cohesive_element_inserter_helper.cc
+++ b/src/mesh_utils/cohesive_element_inserter_helper.cc
@@ -1,947 +1,948 @@
 /**
  * @file   cohesive_element_inserter_helper.cc
  *
  * @author Nicolas Richart
  *
  * @date creation  jeu sep 03 2020
  *
  * @brief A Documented file.
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2011 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 "cohesive_element_inserter_helper.hh"
 /* -------------------------------------------------------------------------- */
 #include "data_accessor.hh"
 #include "element_synchronizer.hh"
 #include "fe_engine.hh"
 #include "mesh_accessor.hh"
 /* -------------------------------------------------------------------------- */
 #include <queue>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 CohesiveElementInserterHelper::CohesiveElementInserterHelper(
     Mesh & mesh, const ElementTypeMapArray<bool> & facet_insertion)
     : doubled_nodes(0, 2, "doubled_nodes"), mesh(mesh),
       mesh_facets(mesh.getMeshFacets()) {
 
   auto spatial_dimension = mesh_facets.getSpatialDimension();
 
   for (auto gt : ghost_types) {
     for (auto type : mesh_facets.elementTypes(_ghost_type = gt)) {
       nb_new_facets(type, gt) = mesh_facets.getNbElement(type, gt);
     }
   }
 
   std::array<Int, 2> nb_facet_to_insert{0, 0};
   // creates a vector of the facets to insert
   std::vector<Element> potential_facets_to_double;
   for (auto gt_facet : ghost_types) {
     for (auto type_facet :
          mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
       const auto & f_insertion = facet_insertion(type_facet, gt_facet);
       auto & counter = nb_facet_to_insert[gt_facet == _not_ghost ? 0 : 1];
 
       for (auto && data : enumerate(f_insertion)) {
         if (std::get<1>(data)) {
           UInt el = std::get<0>(data);
           potential_facets_to_double.push_back({type_facet, el, gt_facet});
           ++counter;
         }
       }
     }
   }
 
   // Defines a global order of insertion oof cohesive element to ensure node
   // doubling appens in the smae order, this is necessary for the global node
   // numbering
   if (mesh_facets.isDistributed()) {
     const auto & comm = mesh_facets.getCommunicator();
 
     ElementTypeMapArray<Int> global_orderings;
 
     global_orderings.initialize(mesh_facets,
                                 _spatial_dimension = spatial_dimension - 1,
                                 _with_nb_element = true, _default_value = -1);
 
     auto starting_index = nb_facet_to_insert[0];
     comm.exclusiveScan(starting_index);
 
     // define the local numbers for facet to insert
     for (auto gt_facet : ghost_types) {
       for (auto type_facet :
            mesh_facets.elementTypes(spatial_dimension - 1, gt_facet)) {
         for (auto data : zip(facet_insertion(type_facet, gt_facet),
                              global_orderings(type_facet, gt_facet))) {
           if (std::get<0>(data)) {
             std::get<1>(data) = starting_index;
             ++starting_index;
           }
         }
       }
     }
 
     // retreives the oorder number from neighoring processors
     auto && synchronizer = mesh_facets.getElementSynchronizer();
     SimpleElementDataAccessor<Int> data_accessor(
         global_orderings, SynchronizationTag::_ce_insertion_order);
     synchronizer.synchronizeOnce(data_accessor,
                                  SynchronizationTag::_ce_insertion_order);
 
     // sort the facets to double based on this global ordering
     std::sort(potential_facets_to_double.begin(),
               potential_facets_to_double.end(),
               [&global_orderings](auto && el1, auto && el2) {
                 return global_orderings(el1) < global_orderings(el2);
               });
   }
 
   for (auto d : arange(spatial_dimension)) {
     facets_to_double_by_dim[d] = std::make_unique<Array<Element>>(
         0, 2, "facets_to_double_" + std::to_string(d));
   }
 
   auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
 
   MeshAccessor mesh_accessor(mesh_facets);
   auto & elements_to_subelements = mesh_accessor.getElementToSubelement();
 
   for (auto && facet_to_double : potential_facets_to_double) {
     auto gt_facet = facet_to_double.ghost_type;
     auto type_facet = facet_to_double.type;
     auto & elements_to_facets = elements_to_subelements(type_facet, gt_facet);
 
     auto & elements_to_facet = elements_to_facets(facet_to_double.element);
     if (elements_to_facet[1].type == _not_defined
 #if defined(AKANTU_COHESIVE_ELEMENT)
         || elements_to_facet[1].kind() == _ek_cohesive
 #endif
     ) {
       AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary");
       continue;
     }
 
     auto new_facet = nb_new_facets(type_facet, gt_facet)++;
 
     facets_to_double.push_back(Vector<Element>{
         facet_to_double, Element{type_facet, new_facet, gt_facet}});
 
     /// update facet_to_element vector
     auto & element_to_update = elements_to_facet[1];
     auto el = element_to_update.element;
 
-    auto & facets_to_elements = mesh_facets.getSubelementToElement(
+    const auto & facets_to_elements = mesh_facets.getSubelementToElement(
         element_to_update.type, element_to_update.ghost_type);
     auto facets_to_element = Vector<Element>(
         make_view(facets_to_elements, facets_to_elements.getNbComponent())
             .begin()[el]);
     auto facet_to_update = std::find(facets_to_element.begin(),
                                      facets_to_element.end(), facet_to_double);
 
     AKANTU_DEBUG_ASSERT(facet_to_update != facets_to_element.end(),
                         "Facet not found");
 
     facet_to_update->element = new_facet;
 
     /// update elements connected to facet
     const auto & first_facet_list = elements_to_facet;
     elements_to_facets.push_back(first_facet_list);
 
     /// set new and original facets as boundary facets
     elements_to_facets(new_facet)[0] = elements_to_facets(new_facet)[1];
     elements_to_facets(new_facet)[1] = ElementNull;
 
     elements_to_facets(facet_to_double.element)[1] = ElementNull;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto
 CohesiveElementInserterHelper::hasElement(const Vector<UInt> & nodes_element,
                                           const Vector<UInt> & nodes) -> bool {
   // if one of the nodes of "nodes" is not in "nodes_element" it stops
   auto it = std::mismatch(nodes.begin(), nodes.end(), nodes_element.begin(),
                           [&](auto && node, auto && /*node2*/) -> bool {
                             auto it = std::find(nodes_element.begin(),
                                                 nodes_element.end(), node);
                             return (it != nodes_element.end());
                           });
 
   // true if all "nodes" where found in "nodes_element"
   return (it.first == nodes.end());
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto CohesiveElementInserterHelper::removeElementsInVector(
     const std::vector<Element> & elem_to_remove,
     std::vector<Element> & elem_list) -> bool {
   if (elem_list.size() <= elem_to_remove.size()) {
     return true;
   }
 
   auto el_it = elem_to_remove.begin();
   auto el_last = elem_to_remove.end();
   std::vector<Element>::iterator el_del;
 
   UInt deletions = 0;
 
   for (; el_it != el_last; ++el_it) {
     el_del = std::find(elem_list.begin(), elem_list.end(), *el_it);
 
     if (el_del != elem_list.end()) {
       elem_list.erase(el_del);
       ++deletions;
     }
   }
 
   AKANTU_DEBUG_ASSERT(deletions == 0 || deletions == elem_to_remove.size(),
                       "Not all elements have been erased");
 
   return deletions == 0;
 }
 
 /* -------------------------------------------------------------------------- */
 void CohesiveElementInserterHelper::updateElementalConnectivity(
     Mesh & mesh, UInt old_node, UInt new_node,
     const std::vector<Element> & element_list,
     const std::vector<Element> * facet_list) {
   AKANTU_DEBUG_IN();
 
-  for (auto & element : element_list) {
+  for (const auto & element : element_list) {
     if (element.type == _not_defined) {
       continue;
     }
 
     Vector<UInt> connectivity = mesh.getConnectivity(element);
 
     if (element.kind() == _ek_cohesive) {
       AKANTU_DEBUG_ASSERT(
           facet_list != nullptr,
           "Provide a facet list in order to update cohesive elements");
 
       const Vector<Element> facets =
           mesh_facets.getSubelementToElement(element);
 
       auto facet_nb_nodes = connectivity.size() / 2;
 
       /// loop over cohesive element's facets
       for (const auto & facet : enumerate(facets)) {
         /// skip facets if not present in the list
         if (std::find(facet_list->begin(), facet_list->end(),
                       std::get<1>(facet)) == facet_list->end()) {
           continue;
         }
 
         auto n = std::get<0>(facet);
 
-        auto begin = connectivity.begin() + static_cast<UInt>(n * facet_nb_nodes);
+        auto begin =
+            connectivity.begin() + static_cast<UInt>(n * facet_nb_nodes);
         auto end = begin + facet_nb_nodes;
 
         auto it = std::find(begin, end, old_node);
         AKANTU_DEBUG_ASSERT(it != end, "Node not found in current element");
 
         *it = new_node;
       }
     } else {
       auto it = std::find(connectivity.begin(), connectivity.end(), old_node);
       AKANTU_DEBUG_ASSERT(it != connectivity.end(),
                           "Node not found in current element");
 
       /// update connectivity
       *it = new_node;
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void CohesiveElementInserterHelper::updateSubelementToElement(UInt dim,
                                                               bool facet_mode) {
   auto & facets_to_double = *facets_to_double_by_dim[dim];
-  auto & facets_to_subfacets =
-      elementsOfDimToElementsOfDim(dim + static_cast<decltype(dim)>(facet_mode), dim);
+  auto & facets_to_subfacets = elementsOfDimToElementsOfDim(
+      dim + static_cast<decltype(dim)>(facet_mode), dim);
 
   for (auto && data :
        zip(make_view(facets_to_double, 2), facets_to_subfacets)) {
     const auto & old_subfacet = std::get<0>(data)[0];
     const auto & new_subfacet = std::get<0>(data)[1];
     auto & facet_to_subfacets = std::get<1>(data);
 
     MeshAccessor mesh_accessor(mesh_facets);
     for (auto & facet : facet_to_subfacets) {
       Vector<Element> subfacets = mesh_accessor.getSubelementToElement(facet);
 
       auto && subfacet_to_update_it =
           std::find(subfacets.begin(), subfacets.end(), old_subfacet);
 
       AKANTU_DEBUG_ASSERT(subfacet_to_update_it != subfacets.end(),
                           "Subfacet not found");
 
       *subfacet_to_update_it = new_subfacet;
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void CohesiveElementInserterHelper::updateElementToSubelement(UInt dim,
                                                               bool facet_mode) {
   auto & facets_to_double = *facets_to_double_by_dim[dim];
   auto & facets_to_subfacets = elementsOfDimToElementsOfDim(
       dim + static_cast<decltype(dim)>(facet_mode), dim);
 
   MeshAccessor mesh_accessor(mesh_facets);
   // resize the arrays
   mesh_accessor.getElementToSubelement().initialize(
       mesh_facets, _spatial_dimension = dim, _with_nb_element = true);
 
   for (auto && data :
        zip(make_view(facets_to_double, 2), facets_to_subfacets)) {
     const auto & new_facet = std::get<0>(data)[1];
     mesh_accessor.getElementToSubelement(new_facet) = std::get<1>(data);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void CohesiveElementInserterHelper::updateQuadraticSegments(UInt dim) {
   AKANTU_DEBUG_IN();
   auto spatial_dimension = mesh.getSpatialDimension();
   auto & facets_to_double = *facets_to_double_by_dim[dim];
 
   MeshAccessor mesh_accessor(mesh_facets);
   auto & connectivities = mesh_accessor.getConnectivities();
 
   /// this ones matter only for segments in 3D
   Array<std::vector<Element>> * element_to_subfacet_double = nullptr;
   Array<std::vector<Element>> * facet_to_subfacet_double = nullptr;
 
   if (dim == spatial_dimension - 2) {
     element_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 2, dim);
     facet_to_subfacet_double = &elementsOfDimToElementsOfDim(dim + 1, dim);
   }
 
-  auto & element_to_subelement = mesh_facets.getElementToSubelement();
+  const auto & element_to_subelement = mesh_facets.getElementToSubelement();
   std::vector<UInt> middle_nodes;
 
   for (auto && facet_to_double : make_view(facets_to_double, 2)) {
     auto & old_facet = facet_to_double[0];
     if (old_facet.type != _segment_3) {
       continue;
     }
 
     auto node = connectivities(old_facet, 2);
     if (not mesh.isPureGhostNode(node)) {
       middle_nodes.push_back(node);
     }
   }
 
   auto n = doubled_nodes.size();
   doubleNodes(middle_nodes);
 
   for (auto && data : enumerate(make_view(facets_to_double, 2))) {
     auto facet = std::get<0>(data);
     auto & old_facet = std::get<1>(data)[0];
     if (old_facet.type != _segment_3) {
       continue;
     }
 
     auto old_node = connectivities(old_facet, 2);
 
     if (mesh.isPureGhostNode(old_node)) {
       continue;
     }
 
     auto new_node = doubled_nodes(n, 1);
 
     auto & new_facet = std::get<1>(data)[1];
     connectivities(new_facet, 2) = new_node;
 
     if (dim == spatial_dimension - 2) {
       updateElementalConnectivity(mesh_facets, old_node, new_node,
                                   element_to_subelement(new_facet, 0));
 
       updateElementalConnectivity(mesh, old_node, new_node,
                                   (*element_to_subfacet_double)(facet),
                                   &(*facet_to_subfacet_double)(facet));
     } else {
       updateElementalConnectivity(mesh, old_node, new_node,
                                   element_to_subelement(new_facet, 0));
     }
 
     ++n;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 UInt CohesiveElementInserterHelper::insertCohesiveElement() {
   auto nb_new_facets = insertFacetsOnly();
   if (nb_new_facets == 0) {
     return 0;
   }
 
   updateCohesiveData();
   return nb_new_facets;
 }
 
 /* -------------------------------------------------------------------------- */
 UInt CohesiveElementInserterHelper::insertFacetsOnly() {
   UInt spatial_dimension = mesh.getSpatialDimension();
 
   if (facets_to_double_by_dim[spatial_dimension - 1]->empty()) {
     return 0;
   }
 
   if (spatial_dimension == 1) {
     doublePointFacet();
   } else if (spatial_dimension == 2) {
     doubleFacets<1>();
     findSubfacetToDouble<1>();
 
     doubleSubfacet<2>();
 
   } else if (spatial_dimension == 3) {
     doubleFacets<2>();
     findSubfacetToDouble<2>();
 
     doubleFacets<1>();
     findSubfacetToDouble<1>();
 
     doubleSubfacet<3>();
   }
 
   return facets_to_double_by_dim[spatial_dimension - 1]->size();
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim> void CohesiveElementInserterHelper::doubleFacets() {
   AKANTU_DEBUG_IN();
   NewElementsEvent new_facets;
   auto spatial_dimension = mesh_facets.getSpatialDimension();
 
   auto & facets_to_double = *facets_to_double_by_dim[dim];
   MeshAccessor mesh_accessor(mesh_facets);
   for (auto && facet_to_double : make_view(facets_to_double, 2)) {
     auto && old_facet = facet_to_double[0];
     auto && new_facet = facet_to_double[1];
 
     auto & facets_connectivities =
         mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type);
     facets_connectivities.resize(
         nb_new_facets(old_facet.type, old_facet.ghost_type));
 
     auto facets_connectivities_begin =
         make_view(facets_connectivities, facets_connectivities.getNbComponent())
             .begin();
 
     // copy the connectivities
     Vector<UInt> new_conn(facets_connectivities_begin[new_facet.element]);
     Vector<UInt> old_conn(facets_connectivities_begin[old_facet.element]);
     new_conn = old_conn;
 
     // this will fail if multiple facet types exists for a given element type
     // \TODO handle multiple sub-facet types
     auto nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(old_facet.type);
 
     auto & subfacets_to_facets = mesh_accessor.getSubelementToElementNC(
         old_facet.type, old_facet.ghost_type);
 
     subfacets_to_facets.resize(
         nb_new_facets(old_facet.type, old_facet.ghost_type), ElementNull);
 
     auto subfacets_to_facets_begin =
         make_view(subfacets_to_facets, nb_subfacet_per_facet).begin();
 
     // copy the subfacet to facets information
     Vector<Element> old_subfacets_to_facet(
         subfacets_to_facets_begin[old_facet.element]);
     Vector<Element> new_subfacet_to_facet(
         subfacets_to_facets_begin[new_facet.element]);
     new_subfacet_to_facet = old_subfacets_to_facet;
 
     for (auto & subfacet : old_subfacets_to_facet) {
       if (subfacet == ElementNull) {
         continue;
       }
 
       /// update facet_to_subfacet array
       mesh_accessor.getElementToSubelement(subfacet).push_back(new_facet);
     }
 
     new_facets.getList().push_back(new_facet);
   }
 
   /// update facet_to_subfacet and _segment_3 facets if any
   if (dim != spatial_dimension - 1) {
     updateSubelementToElement(dim, true);
     updateElementToSubelement(dim, true);
     updateQuadraticSegments(dim);
   } else {
     updateQuadraticSegments(dim);
   }
 
   mesh_facets.sendEvent(new_facets);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim> void CohesiveElementInserterHelper::findSubfacetToDouble() {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh_facets.getSpatialDimension();
 
   MeshAccessor mesh_accessor(mesh_facets);
 
   auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
   auto & subfacets_to_facets = mesh_facets.getSubelementToElement();
   auto & elements_to_facets = mesh_accessor.getElementToSubelement();
 
   /// loop on every facet
   for (auto f : arange(2)) {
     for (auto && facet_to_double : make_view(facets_to_double, 2)) {
       auto old_facet = facet_to_double[f];
       auto new_facet = facet_to_double[1 - f];
 
       const auto & starting_element = elements_to_facets(new_facet, 0)[0];
       auto current_facet = old_facet;
 
       Vector<Element> subfacets_to_facet = subfacets_to_facets.get(old_facet);
       /// loop on every subfacet
       for (auto & subfacet : subfacets_to_facet) {
         if (subfacet == ElementNull) {
           continue;
         }
 
         if (dim == spatial_dimension - 2) {
           Vector<Element> subsubfacets_to_subfacet =
               subfacets_to_facets.get(subfacet);
           /// loop on every subsubfacet
           for (auto & subsubfacet : subsubfacets_to_subfacet) {
             if (subsubfacet == ElementNull) {
               continue;
             }
 
             Vector<UInt> subsubfacet_connectivity =
                 mesh_facets.getConnectivity(subsubfacet);
 
             std::vector<Element> element_list;
             std::vector<Element> facet_list;
             std::vector<Element> subfacet_list;
 
             /// check if subsubfacet is to be doubled
             if (findElementsAroundSubfacet(
                     starting_element, current_facet, subsubfacet_connectivity,
                     element_list, facet_list, &subfacet_list) == false and
                 removeElementsInVector(
                     subfacet_list, elements_to_facets(subsubfacet)) == false) {
               Element new_subsubfacet{
                   subsubfacet.type,
                   nb_new_facets(subsubfacet.type, subsubfacet.ghost_type)++,
                   subsubfacet.ghost_type};
               facets_to_double_by_dim[dim - 1]->push_back(
                   Vector<Element>{subsubfacet, new_subsubfacet});
               elementsOfDimToElementsOfDim(dim - 1, dim - 1)
                   .push_back(subfacet_list);
-              elementsOfDimToElementsOfDim(dim, dim - 1)
-                  .push_back(facet_list);
+              elementsOfDimToElementsOfDim(dim, dim - 1).push_back(facet_list);
               elementsOfDimToElementsOfDim(dim + 1, dim - 1)
                   .push_back(element_list);
             }
           }
         } else {
           std::vector<Element> element_list;
           std::vector<Element> facet_list;
 
           Vector<UInt> subfacet_connectivity =
               mesh_facets.getConnectivity(subfacet);
 
           /// check if subfacet is to be doubled
           if (findElementsAroundSubfacet(starting_element, current_facet,
                                          subfacet_connectivity, element_list,
                                          facet_list) == false and
               removeElementsInVector(facet_list,
                                      elements_to_facets(subfacet)) == false) {
             Element new_subfacet{
                 subfacet.type,
                 nb_new_facets(subfacet.type, subfacet.ghost_type)++,
                 subfacet.ghost_type};
             facets_to_double_by_dim[dim - 1]->push_back(
                 Vector<Element>{subfacet, new_subfacet});
 
             elementsOfDimToElementsOfDim(dim, dim - 1).push_back(facet_list);
             elementsOfDimToElementsOfDim(dim + 1, dim - 1)
                 .push_back(element_list);
           }
         }
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void CohesiveElementInserterHelper::doubleNodes(
     const std::vector<UInt> & old_nodes) {
   auto & position = mesh.getNodes();
   auto spatial_dimension = mesh.getSpatialDimension();
 
   auto old_nb_nodes = position.size();
   position.reserve(old_nb_nodes + old_nodes.size());
   doubled_nodes.reserve(doubled_nodes.size() + old_nodes.size());
 
   auto position_begin = position.begin(spatial_dimension);
   for (auto && data : enumerate(old_nodes)) {
     auto n = std::get<0>(data);
     auto old_node = std::get<1>(data);
     decltype(old_node) new_node = old_nb_nodes + n;
 
     /// store doubled nodes
     doubled_nodes.push_back(Vector<UInt>{old_node, new_node});
 
     /// update position
     Vector<Real> coords = Vector<Real>(position_begin[old_node]);
     position.push_back(coords);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 bool CohesiveElementInserterHelper::findElementsAroundSubfacet(
     const Element & starting_element, const Element & end_facet,
     const Vector<UInt> & subfacet_connectivity,
     std::vector<Element> & element_list, std::vector<Element> & facet_list,
     std::vector<Element> * subfacet_list) {
   bool facet_matched = false;
 
   element_list.push_back(starting_element);
 
   std::queue<Element> elements_to_check;
   elements_to_check.push(starting_element);
 
   /// keep going as long as there are elements to check
   while (not elements_to_check.empty()) {
     /// check current element
     Element & current_element = elements_to_check.front();
 
     const Vector<Element> facets_to_element =
         mesh_facets.getSubelementToElement(current_element);
 
     // for every facet of the element
     for (auto & current_facet : facets_to_element) {
       if (current_facet == ElementNull) {
         continue;
       }
 
       if (current_facet == end_facet) {
         facet_matched = true;
       }
 
       auto find_facet =
           std::find(facet_list.begin(), facet_list.end(), current_facet);
 
       // facet already listed or subfacet_connectivity is not in the
       // connectivity of current_facet;
       if ((find_facet != facet_list.end()) or
           not hasElement(mesh_facets.getConnectivity(current_facet),
                          subfacet_connectivity)) {
         continue;
       }
 
       facet_list.push_back(current_facet);
 
       if (subfacet_list != nullptr) {
         const Vector<Element> subfacets_of_facet =
             mesh_facets.getSubelementToElement(current_facet);
 
         /// check subfacets
         for (const auto & current_subfacet : subfacets_of_facet) {
           if (current_subfacet == ElementNull) {
             continue;
           }
 
           auto find_subfacet = std::find(
               subfacet_list->begin(), subfacet_list->end(), current_subfacet);
           if ((find_subfacet != subfacet_list->end()) or
               not hasElement(mesh_facets.getConnectivity(current_subfacet),
                              subfacet_connectivity)) {
             continue;
           }
 
           subfacet_list->push_back(current_subfacet);
         }
       }
 
       /// consider opposing element
       const auto & elements_to_facet =
           mesh_facets.getElementToSubelement(current_facet);
       UInt opposing = 0;
       if (elements_to_facet[0] == current_element) {
         opposing = 1;
       }
 
       auto & opposing_element = elements_to_facet[opposing];
 
       /// skip null elements since they are on a boundary
       if (opposing_element == ElementNull) {
         continue;
       }
 
       /// skip this element if already added
       if (std::find(element_list.begin(), element_list.end(),
                     opposing_element) != element_list.end()) {
         continue;
       }
 
       /// only regular elements have to be checked
       if (opposing_element.kind() == _ek_regular) {
         elements_to_check.push(opposing_element);
       }
 
       element_list.push_back(opposing_element);
 
       AKANTU_DEBUG_ASSERT(hasElement(mesh.getConnectivity(opposing_element),
                                      subfacet_connectivity),
                           "Subfacet doesn't belong to this element");
     }
 
     /// erased checked element from the list
     elements_to_check.pop();
   }
 
   return facet_matched;
 }
 
 /* -------------------------------------------------------------------------- */
 void CohesiveElementInserterHelper::updateCohesiveData() {
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   bool third_dimension = spatial_dimension == 3;
 
   MeshAccessor mesh_accessor(mesh);
   MeshAccessor mesh_facets_accessor(mesh_facets);
 
   for (auto ghost_type : ghost_types) {
     for (auto cohesive_type : mesh.elementTypes(_element_kind = _ek_cohesive)) {
       auto nb_cohesive_elements = mesh.getNbElement(cohesive_type, ghost_type);
       nb_new_facets(cohesive_type, ghost_type) = nb_cohesive_elements;
 
       mesh_facets_accessor.getSubelementToElement(cohesive_type, ghost_type);
     }
   }
 
   auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
 
   new_elements.reserve(new_elements.size() + facets_to_double.size());
 
   std::array<Element, 2> facets;
 
   auto & element_to_facet = mesh_facets_accessor.getElementToSubelement();
   auto & facets_to_cohesive_elements =
       mesh_facets_accessor.getSubelementToElement();
 
   for (auto && facet_to_double : make_view(facets_to_double, 2)) {
     auto & old_facet = facet_to_double[0];
 
     /// (in 3D cohesive elements connectivity is inverted)
     facets[third_dimension ? 1 : 0] = old_facet;
     facets[third_dimension ? 0 : 1] = facet_to_double[1];
 
     auto type_cohesive = FEEngine::getCohesiveElementType(old_facet.type);
 
     auto & facet_connectivities =
         mesh_facets.getConnectivity(old_facet.type, old_facet.ghost_type);
     auto facet_connectivity_it =
         facet_connectivities.begin(facet_connectivities.getNbComponent());
 
     auto cohesive_element = Element{
         type_cohesive, nb_new_facets(type_cohesive, old_facet.ghost_type)++,
         old_facet.ghost_type};
 
     auto & cohesives_connectivities =
         mesh_accessor.getConnectivity(type_cohesive, old_facet.ghost_type);
     Matrix<UInt> connectivity(facet_connectivities.getNbComponent(), 2);
     Vector<Element> facets_to_cohesive_element(2);
 
     for (auto s : arange(2)) {
       /// store doubled facets
       facets_to_cohesive_element[s] = facets[s];
 
       // update connectivities
       connectivity(s) = Vector<UInt>(facet_connectivity_it[facets[s].element]);
 
       /// update element_to_facet vectors
       element_to_facet(facets[s], 0)[1] = cohesive_element;
     }
 
     cohesives_connectivities.push_back(connectivity);
 
     facets_to_cohesive_elements(type_cohesive, old_facet.ghost_type)
         .push_back(facets_to_cohesive_element);
 
     /// add cohesive element to the element event list
     new_elements.push_back(cohesive_element);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void CohesiveElementInserterHelper::doublePointFacet() {
   UInt spatial_dimension = mesh.getSpatialDimension();
   if (spatial_dimension != 1) {
     return;
   }
 
   NewElementsEvent new_facets_event;
 
   auto & facets_to_double = *facets_to_double_by_dim[spatial_dimension - 1];
-  auto & element_to_facet = mesh_facets.getElementToSubelement();
+  const auto & element_to_facet = mesh_facets.getElementToSubelement();
   auto & position = mesh.getNodes();
   MeshAccessor mesh_accessor(mesh_facets);
 
   for (auto ghost_type : ghost_types) {
     for (auto facet_type : nb_new_facets.elementTypes(
              spatial_dimension - 1, ghost_type, _ek_regular)) {
       auto nb_new_element = nb_new_facets(facet_type, ghost_type);
       auto & connectivities =
           mesh_accessor.getConnectivity(facet_type, ghost_type);
-      connectivities.resize(connectivities.size() + nb_new_element);
+      connectivities.resize(nb_new_element);
     }
   }
 
+  position.reserve(position.size() + facets_to_double.size());
   for (auto facet_to_double : make_view(facets_to_double, 2)) {
     auto & old_facet = facet_to_double[0];
     auto & new_facet = facet_to_double[1];
 
     auto element = element_to_facet(new_facet)[0];
     auto & facet_connectivities =
         mesh_accessor.getConnectivity(old_facet.type, old_facet.ghost_type);
     auto old_node = facet_connectivities(old_facet.element);
     auto new_node = position.size();
 
     /// update position
     position.push_back(position(old_node));
     facet_connectivities(new_facet.element) = new_node;
 
     Vector<UInt> segment_conectivity = mesh.getConnectivity(element);
 
     /// update facet connectivity
     auto it = std::find(segment_conectivity.begin(), segment_conectivity.end(),
                         old_node);
     *it = new_node;
 
     doubled_nodes.push_back(Vector<UInt>{old_node, new_node});
     new_facets_event.getList().push_back(new_facet);
   }
 
   mesh_facets.sendEvent(new_facets_event);
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt spatial_dimension>
 void CohesiveElementInserterHelper::doubleSubfacet() {
   if (spatial_dimension == 1) {
     return;
   }
 
   NewElementsEvent new_facets_event;
 
   std::vector<UInt> nodes_to_double;
   MeshAccessor mesh_accessor(mesh_facets);
   auto & connectivities = mesh_accessor.getConnectivities();
 
   auto & facets_to_double = *facets_to_double_by_dim[0];
   ElementTypeMap<Int> nb_new_elements;
   for (auto && facet_to_double : make_view(facets_to_double, 2)) {
     auto && old_element = facet_to_double[0];
     nb_new_elements(old_element.type, old_element.ghost_type) = 0;
   }
 
   for (auto && facet_to_double : make_view(facets_to_double, 2)) {
     auto && old_element = facet_to_double[0];
     ++nb_new_elements(old_element.type, old_element.ghost_type);
   }
 
   for (auto ghost_type : ghost_types) {
     for (auto facet_type :
          nb_new_elements.elementTypes(0, ghost_type, _ek_regular)) {
       auto & connectivities =
           mesh_accessor.getConnectivity(facet_type, ghost_type);
       connectivities.resize(connectivities.size() +
                             nb_new_elements(facet_type, ghost_type));
     }
   }
 
   for (auto && facet_to_double : make_view(facets_to_double, 2)) {
     auto & old_facet = facet_to_double(0);
     // auto & new_facet = facet_to_double(1);
 
     AKANTU_DEBUG_ASSERT(
         old_facet.type == _point_1,
         "Only _point_1 subfacet doubling is supported at the moment");
 
     /// list nodes double
     nodes_to_double.push_back(connectivities(old_facet));
   }
 
   auto old_nb_doubled_nodes = doubled_nodes.size();
   doubleNodes(nodes_to_double);
 
   auto double_nodes_view = make_view(doubled_nodes, 2);
 
   for (auto && data :
        zip(make_view(facets_to_double, 2),
            range(double_nodes_view.begin() + old_nb_doubled_nodes,
                  double_nodes_view.end()),
            arange(facets_to_double.size()))) {
     // auto & old_facet = std::get<0>(data)[0];
     auto & new_facet = std::get<0>(data)[1];
 
     new_facets_event.getList().push_back(new_facet);
 
     auto & nodes = std::get<1>(data);
     auto old_node = nodes(0);
     auto new_node = nodes(1);
 
     auto f = std::get<2>(data);
 
     /// add new nodes in connectivity
     connectivities(new_facet) = new_node;
 
     updateElementalConnectivity(mesh, old_node, new_node,
                                 elementsOfDimToElementsOfDim(2, 0)(f),
                                 &elementsOfDimToElementsOfDim(1, 0)(f));
 
     updateElementalConnectivity(mesh_facets, old_node, new_node,
                                 elementsOfDimToElementsOfDim(1, 0)(f));
 
     if (spatial_dimension == 3) {
       updateElementalConnectivity(mesh_facets, old_node, new_node,
                                   elementsOfDimToElementsOfDim(0, 0)(f));
     }
   }
 
   updateSubelementToElement(0, spatial_dimension == 2);
   updateElementToSubelement(0, spatial_dimension == 2);
 
   mesh_facets.sendEvent(new_facets_event);
 }
 
 } // namespace akantu
diff --git a/src/model/common/model_solver.cc b/src/model/common/model_solver.cc
index dcbcd6f6f..12531d0f8 100644
--- a/src/model/common/model_solver.cc
+++ b/src/model/common/model_solver.cc
@@ -1,395 +1,406 @@
 /**
  * @file   model_solver.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Aug 18 2015
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  Implementation of ModelSolver
  *
  *
  * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "model_solver.hh"
 #include "dof_manager.hh"
 #include "dof_manager_default.hh"
 #include "mesh.hh"
 #include "non_linear_solver.hh"
 #include "time_step_solver.hh"
 
 #if defined(AKANTU_USE_PETSC)
 #include "dof_manager_petsc.hh"
 #endif
 
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <typename T> static T getOptionToType(const std::string & opt_str) {
   std::stringstream sstr(opt_str);
   T opt;
   sstr >> opt;
 
   return opt;
 }
 
 /* -------------------------------------------------------------------------- */
 ModelSolver::ModelSolver(Mesh & mesh, const ModelType & type, const ID & id,
-                         UInt memory_id,
-                         std::shared_ptr<DOFManager> dof_manager)
+                         UInt memory_id)
     : Parsable(ParserType::_model, id), SolverCallback(), model_type(type),
       parent_id(id), parent_memory_id(memory_id), mesh(mesh) {
+}
+
+/* -------------------------------------------------------------------------- */
+ModelSolver::ModelSolver(Mesh & mesh, const ModelType & type, const ID & id,
+                         UInt memory_id,
+                         std::shared_ptr<DOFManager> dof_manager)
+    : ModelSolver(mesh, type, id, memory_id) {
   if (not dof_manager) {
-    initDOFManager();
+    this->initDOFManager();
   } else {
     this->dof_manager = dof_manager;
     this->setDOFManager(*this->dof_manager);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 ModelSolver::~ModelSolver() = default;
 
 /* -------------------------------------------------------------------------- */
 std::tuple<ParserSection, bool> ModelSolver::getParserSection() {
   auto sub_sections = getStaticParser().getSubSections(ParserType::_model);
 
   auto it = std::find_if(
       sub_sections.begin(), sub_sections.end(), [&](auto && section) {
         auto type = getOptionToType<ModelType>(section.getName());
         // default id should be the model type if not defined
         std::string name = section.getParameter("name", this->parent_id);
         return type == model_type and name == this->parent_id;
       });
 
   if (it == sub_sections.end()) {
     return std::make_tuple(ParserSection(), true);
   }
 
   return std::make_tuple(*it, false);
 }
 
 /* -------------------------------------------------------------------------- */
-void ModelSolver::initDOFManager() {
+std::shared_ptr<DOFManager> ModelSolver::initDOFManager() {
   // default without external solver activated at compilation same as mumps that
   // is the historical solver but with only the lumped solver
   ID solver_type = "default";
 
 #if defined(AKANTU_USE_MUMPS)
   solver_type = "default";
 #elif defined(AKANTU_USE_PETSC)
   solver_type = "petsc";
 #endif
 
   ParserSection section;
   bool is_empty;
   std::tie(section, is_empty) = this->getParserSection();
 
   if (not is_empty) {
     solver_type = section.getOption(solver_type);
-    this->initDOFManager(section, solver_type);
+    return this->initDOFManager(section, solver_type);
   } else {
-    this->initDOFManager(solver_type);
+    return this->initDOFManager(solver_type);
   }
 }
 
 /* -------------------------------------------------------------------------- */
-void ModelSolver::initDOFManager(const ID & solver_type) {
+std::shared_ptr<DOFManager>
+ModelSolver::initDOFManager(const ID & solver_type) {
   if (dof_manager) {
     AKANTU_EXCEPTION("The DOF manager for this model is already initialized !");
   }
 
   try {
     this->dof_manager = DOFManagerFactory::getInstance().allocate(
         solver_type, mesh, this->parent_id + ":dof_manager_" + solver_type,
         this->parent_memory_id);
   } catch (...) {
     AKANTU_EXCEPTION(
         "To use the solver "
         << solver_type
         << " you will have to code it. This is an unknown solver type.");
   }
 
   this->setDOFManager(*this->dof_manager);
+  return this->dof_manager;
 }
 
 /* -------------------------------------------------------------------------- */
-void ModelSolver::initDOFManager(const ParserSection & section,
-                                 const ID & solver_type) {
+std::shared_ptr<DOFManager>
+ModelSolver::initDOFManager(const ParserSection & section,
+                            const ID & solver_type) {
   this->initDOFManager(solver_type);
   auto sub_sections = section.getSubSections(ParserType::_time_step_solver);
 
   // parsing the time step solvers
   for (auto && section : sub_sections) {
     ID type = section.getName();
     ID solver_id = section.getParameter("name", type);
 
     auto tss_type = getOptionToType<TimeStepSolverType>(type);
     auto tss_options = this->getDefaultSolverOptions(tss_type);
 
     auto sub_solvers_sect =
         section.getSubSections(ParserType::_non_linear_solver);
     auto nb_non_linear_solver_section =
         section.getNbSubSections(ParserType::_non_linear_solver);
 
     auto nls_type = tss_options.non_linear_solver_type;
 
     if (nb_non_linear_solver_section == 1) {
       auto && nls_section = *(sub_solvers_sect.first);
       nls_type = getOptionToType<NonLinearSolverType>(nls_section.getName());
     } else if (nb_non_linear_solver_section > 0) {
       AKANTU_EXCEPTION("More than one non linear solver are provided for the "
                        "time step solver "
                        << solver_id);
     }
 
     this->getNewSolver(solver_id, tss_type, nls_type);
     if (nb_non_linear_solver_section == 1) {
       const auto & nls_section = *(sub_solvers_sect.first);
       this->dof_manager->getNonLinearSolver(solver_id).parseSection(
           nls_section);
     }
 
     auto sub_integrator_sections =
         section.getSubSections(ParserType::_integration_scheme);
 
     for (auto && is_section : sub_integrator_sections) {
       const auto & dof_type_str = is_section.getName();
       ID dof_id;
       try {
         ID tmp = is_section.getParameter("name");
         dof_id = tmp;
       } catch (...) {
         AKANTU_EXCEPTION("No degree of freedom name specified for the "
                          "integration scheme of type "
                          << dof_type_str);
       }
 
       auto it_type = getOptionToType<IntegrationSchemeType>(dof_type_str);
 
       IntegrationScheme::SolutionType s_type = is_section.getParameter(
           "solution_type", tss_options.solution_type[dof_id]);
       this->setIntegrationScheme(solver_id, dof_id, it_type, s_type);
     }
 
     for (auto & is_type : tss_options.integration_scheme_type) {
       if (!this->hasIntegrationScheme(solver_id, is_type.first)) {
         this->setIntegrationScheme(solver_id, is_type.first, is_type.second,
                                    tss_options.solution_type[is_type.first]);
       }
     }
   }
 
   if (section.hasParameter("default_solver")) {
     ID default_solver = section.getParameter("default_solver");
     if (this->hasSolver(default_solver)) {
       this->setDefaultSolver(default_solver);
     } else {
       AKANTU_EXCEPTION(
           "The solver \""
           << default_solver
           << "\" was not created, it cannot be set as default solver");
     }
   }
+
+  return this->dof_manager;
 }
 
 /* -------------------------------------------------------------------------- */
 TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) {
   ID tmp_solver_id = solver_id;
   if (tmp_solver_id.empty()) {
     tmp_solver_id = this->default_solver_id;
   }
 
   TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id);
   return tss;
 }
 
 /* -------------------------------------------------------------------------- */
 const TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) const {
   ID tmp_solver_id = solver_id;
   if (solver_id.empty()) {
     tmp_solver_id = this->default_solver_id;
   }
 
   const TimeStepSolver & tss =
       this->dof_manager->getTimeStepSolver(tmp_solver_id);
   return tss;
 }
 
 /* -------------------------------------------------------------------------- */
 TimeStepSolver & ModelSolver::getTimeStepSolver(const ID & solver_id) {
   return this->getSolver(solver_id);
 }
 
 /* -------------------------------------------------------------------------- */
 const TimeStepSolver &
 ModelSolver::getTimeStepSolver(const ID & solver_id) const {
   return this->getSolver(solver_id);
 }
 
 /* -------------------------------------------------------------------------- */
 NonLinearSolver & ModelSolver::getNonLinearSolver(const ID & solver_id) {
   return this->getSolver(solver_id).getNonLinearSolver();
 }
 /* -------------------------------------------------------------------------- */
 const NonLinearSolver &
 ModelSolver::getNonLinearSolver(const ID & solver_id) const {
   return this->getSolver(solver_id).getNonLinearSolver();
 }
 
 /* -------------------------------------------------------------------------- */
 bool ModelSolver::hasSolver(const ID & solver_id) const {
   ID tmp_solver_id = solver_id;
   if (solver_id.empty()) {
     tmp_solver_id = this->default_solver_id;
   }
 
   if (not this->dof_manager) {
     AKANTU_EXCEPTION("No DOF manager was initialized");
   }
   return this->dof_manager->hasTimeStepSolver(tmp_solver_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void ModelSolver::setDefaultSolver(const ID & solver_id) {
   AKANTU_DEBUG_ASSERT(
       this->hasSolver(solver_id),
       "Cannot set the default solver to a solver that does not exists");
   this->default_solver_id = solver_id;
 }
 
 /* -------------------------------------------------------------------------- */
 void ModelSolver::solveStep(SolverCallback & callback, const ID & solver_id) {
   AKANTU_DEBUG_IN();
 
   TimeStepSolver & tss = this->getSolver(solver_id);
   // make one non linear solve
   tss.solveStep(callback);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void ModelSolver::solveStep(const ID & solver_id) {
   solveStep(*this, solver_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void ModelSolver::getNewSolver(const ID & solver_id,
                                TimeStepSolverType time_step_solver_type,
                                NonLinearSolverType non_linear_solver_type) {
   if (this->default_solver_id.empty()) {
     this->default_solver_id = solver_id;
   }
 
   if (non_linear_solver_type == NonLinearSolverType::_auto) {
     switch (time_step_solver_type) {
     case TimeStepSolverType::_dynamic:
     case TimeStepSolverType::_static:
       non_linear_solver_type = NonLinearSolverType::_newton_raphson;
       break;
     case TimeStepSolverType::_dynamic_lumped:
       non_linear_solver_type = NonLinearSolverType::_lumped;
       break;
     case TimeStepSolverType::_not_defined:
       AKANTU_EXCEPTION(time_step_solver_type
                        << " is not a valid time step solver type");
       break;
     }
   }
 
   this->initSolver(time_step_solver_type, non_linear_solver_type);
 
   NonLinearSolver & nls = this->dof_manager->getNewNonLinearSolver(
       solver_id, non_linear_solver_type);
 
   this->dof_manager->getNewTimeStepSolver(solver_id, time_step_solver_type, nls,
                                           *this);
 }
 
 /* -------------------------------------------------------------------------- */
 Real ModelSolver::getTimeStep(const ID & solver_id) const {
   const TimeStepSolver & tss = this->getSolver(solver_id);
 
   return tss.getTimeStep();
 }
 
 /* -------------------------------------------------------------------------- */
 void ModelSolver::setTimeStep(Real time_step, const ID & solver_id) {
   TimeStepSolver & tss = this->getSolver(solver_id);
 
   return tss.setTimeStep(time_step);
 }
 
 /* -------------------------------------------------------------------------- */
 void ModelSolver::setIntegrationScheme(
     const ID & solver_id, const ID & dof_id,
     const IntegrationSchemeType & integration_scheme_type,
     IntegrationScheme::SolutionType solution_type) {
   TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id);
 
   tss.setIntegrationScheme(dof_id, integration_scheme_type, solution_type);
 }
 
 /* -------------------------------------------------------------------------- */
-  void ModelSolver::setIntegrationScheme(const ID & solver_id, const ID & dof_id,
-			  std::unique_ptr<IntegrationScheme> & integration_scheme,
-			  IntegrationScheme::SolutionType solution_type) {
+void ModelSolver::setIntegrationScheme(
+    const ID & solver_id, const ID & dof_id,
+    std::unique_ptr<IntegrationScheme> & integration_scheme,
+    IntegrationScheme::SolutionType solution_type) {
   TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id);
-  
+
   tss.setIntegrationScheme(dof_id, integration_scheme, solution_type);
-  }
-  
+}
 
 /* -------------------------------------------------------------------------- */
 bool ModelSolver::hasDefaultSolver() const {
   return (not this->default_solver_id.empty());
 }
 
 /* -------------------------------------------------------------------------- */
 bool ModelSolver::hasIntegrationScheme(const ID & solver_id,
                                        const ID & dof_id) const {
   TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id);
   return tss.hasIntegrationScheme(dof_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void ModelSolver::predictor() {}
 
 /* -------------------------------------------------------------------------- */
 void ModelSolver::corrector() {}
 
 /* -------------------------------------------------------------------------- */
 TimeStepSolverType ModelSolver::getDefaultSolverType() const {
   return TimeStepSolverType::_dynamic_lumped;
 }
 
 /* -------------------------------------------------------------------------- */
 ModelSolverOptions
 ModelSolver::getDefaultSolverOptions(__attribute__((unused))
                                      const TimeStepSolverType & type) const {
   ModelSolverOptions options;
   options.non_linear_solver_type = NonLinearSolverType::_auto;
   return options;
 }
 
 } // namespace akantu
diff --git a/src/model/common/model_solver.hh b/src/model/common/model_solver.hh
index 9273e0ec0..4f8e79856 100644
--- a/src/model/common/model_solver.hh
+++ b/src/model/common/model_solver.hh
@@ -1,199 +1,204 @@
 /**
  * @file   model_solver.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Wed Jan 31 2018
  *
  * @brief  Class regrouping the common solve interface to the different models
  *
  *
  * 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 "integration_scheme.hh"
 #include "parsable.hh"
 #include "solver_callback.hh"
 #include "synchronizer_registry.hh"
 /* -------------------------------------------------------------------------- */
 #include <set>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_MODEL_SOLVER_HH_
 #define AKANTU_MODEL_SOLVER_HH_
 
 namespace akantu {
 class Mesh;
 class DOFManager;
 class TimeStepSolver;
 class NonLinearSolver;
 struct ModelSolverOptions;
 } // namespace akantu
 
 namespace akantu {
 
 class ModelSolver : public Parsable,
                     public SolverCallback,
                     public SynchronizerRegistry {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   ModelSolver(Mesh & mesh, const ModelType & type, const ID & id,
-              UInt memory_id, std::shared_ptr<DOFManager> dof_manager=nullptr);
+              UInt memory_id);
+  ModelSolver(Mesh & mesh, const ModelType & type, const ID & id,
+              UInt memory_id,
+              std::shared_ptr<DOFManager> dof_manager);
   ~ModelSolver() override;
 
   /// initialize the dof manager based on solver type passed in the input file
-  void initDOFManager();
+  std::shared_ptr<DOFManager> initDOFManager();
   /// initialize the dof manager based on the used chosen solver type
-  void initDOFManager(const ID & solver_type);
+  std::shared_ptr<DOFManager> initDOFManager(const ID & solver_type);
 
 protected:
   /// initialize the dof manager based on the used chosen solver type
-  void initDOFManager(const ParserSection & section, const ID & solver_type);
+  std::shared_ptr<DOFManager> initDOFManager(const ParserSection & section,
+                                             const ID & solver_type);
 
 public:
   /// Callback for the model to instantiate the matricees when needed
   virtual void initSolver(TimeStepSolverType /*time_step_solver_type*/,
                           NonLinearSolverType /*non_linear_solver_type*/) {}
 
   /// get the section in the input file (if it exsits) corresponding to this
   /// model
   std::tuple<ParserSection, bool> getParserSection();
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// solve a step using a given pre instantiated time step solver and
   /// non linear solver
   virtual void solveStep(const ID & solver_id = "");
 
   /// solve a step using a given pre instantiated time step solver and
   /// non linear solver with a user defined callback instead of the
   /// model itself /!\ This can mess up everything
   virtual void solveStep(SolverCallback & callback, const ID & solver_id = "");
 
   /// Initialize a time solver that can be used afterwards with its id
   void getNewSolver(
       const ID & solver_id, TimeStepSolverType time_step_solver_type,
       NonLinearSolverType non_linear_solver_type = NonLinearSolverType::_auto);
 
   /// set an integration scheme for a given dof and a given solver
   void
   setIntegrationScheme(const ID & solver_id, const ID & dof_id,
                        const IntegrationSchemeType & integration_scheme_type,
                        IntegrationScheme::SolutionType solution_type =
                            IntegrationScheme::_not_defined);
 
   /// set an externally instantiated integration scheme
-  void setIntegrationScheme(const ID & solver_id, const ID & dof_id,
-                            std::unique_ptr<IntegrationScheme> & integration_scheme,
-                            IntegrationScheme::SolutionType solution_type =
-                                IntegrationScheme::_not_defined);
+  void
+  setIntegrationScheme(const ID & solver_id, const ID & dof_id,
+                       std::unique_ptr<IntegrationScheme> & integration_scheme,
+                       IntegrationScheme::SolutionType solution_type =
+                           IntegrationScheme::_not_defined);
 
   /* ------------------------------------------------------------------------ */
   /* SolverCallback interface                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   /// Predictor interface for the callback
   void predictor() override;
 
   /// Corrector interface for the callback
   void corrector() override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// Default time step solver to instantiate for this model
   virtual TimeStepSolverType getDefaultSolverType() const;
 
   /// Default configurations for a given time step solver
   virtual ModelSolverOptions
   getDefaultSolverOptions(const TimeStepSolverType & type) const;
 
   /// get access to the internal dof manager
   DOFManager & getDOFManager() { return *this->dof_manager; }
 
   /// get the time step of a given solver
   Real getTimeStep(const ID & solver_id = "") const;
   /// set the time step of a given solver
   virtual void setTimeStep(Real time_step, const ID & solver_id = "");
 
   /// set the parameter 'param' of the solver 'solver_id'
   // template <typename T>
   // void set(const ID & param, const T & value, const ID & solver_id = "");
 
   /// get the parameter 'param' of the solver 'solver_id'
   // const Parameter & get(const ID & param, const ID & solver_id = "") const;
 
   /// answer to the question "does the solver exists ?"
   bool hasSolver(const ID & solver_id) const;
 
   /// changes the current default solver
   void setDefaultSolver(const ID & solver_id);
 
   /// is a default solver defined
   bool hasDefaultSolver() const;
 
   /// is an integration scheme set for a given solver and a given dof
   bool hasIntegrationScheme(const ID & solver_id, const ID & dof_id) const;
 
   TimeStepSolver & getTimeStepSolver(const ID & solver_id = "");
   NonLinearSolver & getNonLinearSolver(const ID & solver_id = "");
 
   const TimeStepSolver & getTimeStepSolver(const ID & solver_id = "") const;
   const NonLinearSolver & getNonLinearSolver(const ID & solver_id = "") const;
 
 private:
   TimeStepSolver & getSolver(const ID & solver_id);
   const TimeStepSolver & getSolver(const ID & solver_id) const;
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   ModelType model_type;
 
   /// Underlying dof_manager (the brain...)
   std::shared_ptr<DOFManager> dof_manager;
 
 private:
   ID parent_id;
   UInt parent_memory_id;
 
   /// Underlying mesh
   Mesh & mesh;
 
   /// Default time step solver to use
   ID default_solver_id{""};
 };
 
 struct ModelSolverOptions {
   NonLinearSolverType non_linear_solver_type;
   std::map<ID, IntegrationSchemeType> integration_scheme_type;
   std::map<ID, IntegrationScheme::SolutionType> solution_type;
 };
 
 } // namespace akantu
 
 #endif /* AKANTU_MODEL_SOLVER_HH_ */
diff --git a/src/model/model.cc b/src/model/model.cc
index bc948054c..02c111a8e 100644
--- a/src/model/model.cc
+++ b/src/model/model.cc
@@ -1,347 +1,349 @@
 /**
  * @file   model.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Mon Oct 03 2011
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  implementation of model common parts
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "model.hh"
 #include "communicator.hh"
 #include "data_accessor.hh"
 #include "element_group.hh"
 #include "element_synchronizer.hh"
 #include "synchronizer_registry.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 Model::Model(Mesh & mesh, const ModelType & type,
              std::shared_ptr<DOFManager> dof_manager, UInt dim, const ID & id,
              const MemoryID & memory_id)
     : Memory(id, memory_id),
-      ModelSolver(mesh, type, id, memory_id, dof_manager), mesh(mesh),
+      ModelSolver(mesh, type, id, memory_id, std::move(dof_manager)),
+      mesh(mesh),
       spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension()
                                                : dim),
       parser(getStaticParser()) {
-  AKANTU_DEBUG_IN();
-
   this->mesh.registerEventHandler(*this, _ehp_model);
-
-  AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 Model::Model(Mesh & mesh, const ModelType & type, UInt dim, const ID & id,
              const MemoryID & memory_id)
-    : Model(mesh, type, nullptr, dim, id, memory_id) {}
+    : Memory(id, memory_id), ModelSolver(mesh, type, id, memory_id), mesh(mesh),
+      spatial_dimension(dim == _all_dimensions ? mesh.getSpatialDimension()
+                                               : dim),
+      parser(getStaticParser()) {
+  this->mesh.registerEventHandler(*this, _ehp_model);
+}
 
 /* -------------------------------------------------------------------------- */
 Model::~Model() = default;
 
 /* -------------------------------------------------------------------------- */
 void Model::initFullImpl(const ModelOptions & options) {
   AKANTU_DEBUG_IN();
 
   method = options.analysis_method;
   if (!this->hasDefaultSolver()) {
     this->initNewSolver(this->method);
   }
 
   initModel();
 
   initFEEngineBoundary();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::initNewSolver(const AnalysisMethod & method) {
   ID solver_name;
   TimeStepSolverType tss_type;
   std::tie(solver_name, tss_type) = this->getDefaultSolverID(method);
 
   if (not this->hasSolver(solver_name)) {
     ModelSolverOptions options = this->getDefaultSolverOptions(tss_type);
     this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type);
 
     for (auto && is_type : options.integration_scheme_type) {
       if (!this->hasIntegrationScheme(solver_name, is_type.first)) {
         this->setIntegrationScheme(solver_name, is_type.first, is_type.second,
                                    options.solution_type[is_type.first]);
       }
     }
   }
 
   this->method = method;
   this->setDefaultSolver(solver_name);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::initFEEngineBoundary() {
   FEEngine & fem_boundary = getFEEngineBoundary();
   fem_boundary.initShapeFunctions(_not_ghost);
   fem_boundary.initShapeFunctions(_ghost);
 
   fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost);
   fem_boundary.computeNormalsOnIntegrationPoints(_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::dumpGroup(const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.dump();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::dumpGroup(const std::string & group_name,
                       const std::string & dumper_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.dump(dumper_name);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::dumpGroup() {
   for (auto & group : mesh.iterateElementGroups()) {
     group.dump();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::setGroupDirectory(const std::string & directory) {
   for (auto & group : mesh.iterateElementGroups()) {
     group.setDirectory(directory);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::setGroupDirectory(const std::string & directory,
                               const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.setDirectory(directory);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::setGroupBaseName(const std::string & basename,
                              const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.setBaseName(basename);
 }
 
 /* -------------------------------------------------------------------------- */
 DumperIOHelper & Model::getGroupDumper(const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   return group.getDumper();
 }
 
 /* -------------------------------------------------------------------------- */
 // DUMPER stuff
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldToDumper(const std::string & field_id,
                                       std::shared_ptr<dumpers::Field> field,
                                       DumperIOHelper & dumper) {
 #ifdef AKANTU_USE_IOHELPER
   dumper.registerField(field_id, std::move(field));
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpField(const std::string & field_id) {
 
   this->addDumpFieldToDumper(mesh.getDefaultDumperName(), field_id);
 }
 /* -------------------------------------------------------------------------- */
 
 void Model::addDumpFieldVector(const std::string & field_id) {
 
   this->addDumpFieldVectorToDumper(mesh.getDefaultDumperName(), field_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpFieldTensor(const std::string & field_id) {
 
   this->addDumpFieldTensorToDumper(mesh.getDefaultDumperName(), field_id);
 }
 
 /* -------------------------------------------------------------------------- */
 
 void Model::setBaseName(const std::string & field_id) {
 
   mesh.setBaseName(field_id);
 }
 /* -------------------------------------------------------------------------- */
 
 void Model::setBaseNameToDumper(const std::string & dumper_name,
                                 const std::string & basename) {
   mesh.setBaseNameToDumper(dumper_name, basename);
 }
 /* -------------------------------------------------------------------------- */
 
 void Model::addDumpFieldToDumper(const std::string & dumper_name,
                                  const std::string & field_id) {
 
   this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
                                   false);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupField(const std::string & field_id,
                               const std::string & group_name) {
 
   ElementGroup & group = mesh.getElementGroup(group_name);
   this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id,
                                   group_name, _ek_regular, false);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::removeDumpGroupField(const std::string & field_id,
                                  const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id,
                                        group_name);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::removeDumpGroupFieldFromDumper(const std::string & dumper_name,
                                            const std::string & field_id,
                                            const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   group.removeDumpFieldFromDumper(dumper_name, field_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpFieldVectorToDumper(const std::string & dumper_name,
                                        const std::string & field_id) {
   this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
                                   true);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldVector(const std::string & field_id,
                                     const std::string & group_name) {
   ElementGroup & group = mesh.getElementGroup(group_name);
   this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id,
                                         group_name);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldVectorToDumper(const std::string & dumper_name,
                                             const std::string & field_id,
                                             const std::string & group_name) {
   this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name,
                                   _ek_regular, true);
 }
 /* -------------------------------------------------------------------------- */
 
 void Model::addDumpFieldTensorToDumper(const std::string & dumper_name,
                                        const std::string & field_id) {
   this->addDumpGroupFieldToDumper(dumper_name, field_id, "all", _ek_regular,
                                   true);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldToDumper(const std::string & dumper_name,
                                       const std::string & field_id,
                                       const std::string & group_name,
                                       ElementKind element_kind,
                                       bool padding_flag) {
   this->addDumpGroupFieldToDumper(dumper_name, field_id, group_name,
                                   this->spatial_dimension, element_kind,
                                   padding_flag);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumpGroupFieldToDumper(const std::string & dumper_name,
                                       const std::string & field_id,
                                       const std::string & group_name,
                                       UInt spatial_dimension,
                                       ElementKind element_kind,
                                       bool padding_flag) {
 
 #ifdef AKANTU_USE_IOHELPER
   std::shared_ptr<dumpers::Field> field;
 
   if (!field) {
     field = this->createNodalFieldReal(field_id, group_name, padding_flag);
   }
   if (!field) {
     field = this->createNodalFieldUInt(field_id, group_name, padding_flag);
   }
   if (!field) {
     field = this->createNodalFieldBool(field_id, group_name, padding_flag);
   }
   if (!field) {
     field = this->createElementalField(field_id, group_name, padding_flag,
                                        spatial_dimension, element_kind);
   }
   if (!field) {
     field = this->mesh.createFieldFromAttachedData<UInt>(field_id, group_name,
                                                          element_kind);
   }
   if (!field) {
     field = this->mesh.createFieldFromAttachedData<Real>(field_id, group_name,
                                                          element_kind);
   }
 
 #ifndef AKANTU_NDEBUG
   if (!field) {
     AKANTU_DEBUG_WARNING("No field could be found based on name: " << field_id);
   }
 #endif
   if (field) {
     DumperIOHelper & dumper = mesh.getGroupDumper(dumper_name, group_name);
     this->addDumpGroupFieldToDumper(field_id, field, dumper);
   }
 
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 
 void Model::dump() { mesh.dump(); }
 
 /* -------------------------------------------------------------------------- */
 
 void Model::setDirectory(const std::string & directory) {
   mesh.setDirectory(directory);
 }
 
 /* -------------------------------------------------------------------------- */
 
 void Model::setDirectoryToDumper(const std::string & dumper_name,
                                  const std::string & directory) {
   mesh.setDirectoryToDumper(dumper_name, directory);
 }
 
 /* -------------------------------------------------------------------------- */
 
 void Model::setTextModeToDumper() { mesh.setTextModeToDumper(); }
 
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu
diff --git a/src/model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh b/src/model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh
index e9bdee79e..8a9e28ed3 100644
--- a/src/model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh
+++ b/src/model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh
@@ -1,121 +1,125 @@
 /* -------------------------------------------------------------------------- */
 #include "material_von_mises_mazars.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 /* -------------------------------------------------------------------------- */
 template <UInt spatial_dimension, template <UInt> class Parent>
 inline void
 MaterialVonMisesMazars<spatial_dimension, Parent>::computeStressOnQuad(
     const Matrix<Real> & grad_u, Matrix<Real> & sigma, Real & dam,
     Real & Ehat) {
   Matrix<Real> epsilon(3, 3);
   epsilon.zero();
 
-  for (UInt i = 0; i < spatial_dimension; ++i)
-    for (UInt j = 0; j < spatial_dimension; ++j)
+  for (UInt i = 0; i < spatial_dimension; ++i) {
+    for (UInt j = 0; j < spatial_dimension; ++j) {
       epsilon(i, j) = .5 * (grad_u(i, j) + grad_u(j, i));
+    }
+  }
 
   Vector<Real> Fdiag(3);
   Math::matrixEig(3, epsilon.storage(), Fdiag.storage());
 
   Ehat = 0.;
   for (UInt i = 0; i < 3; ++i) {
     Real epsilon_p = std::max(Real(0.), Fdiag(i));
     Ehat += epsilon_p * epsilon_p;
   }
   Ehat = sqrt(Ehat);
 
   // MaterialElastic<spatial_dimension>::computeStressOnQuad(grad_u, sigma);
 
   if (damage_in_compute_stress) {
     computeDamageOnQuad(Ehat, sigma, Fdiag, dam);
   }
 
   if (not this->is_non_local) {
     computeDamageAndStressOnQuad(grad_u, sigma, dam, Ehat);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt spatial_dimension, template <UInt> class Parent>
 inline void
 MaterialVonMisesMazars<spatial_dimension, Parent>::computeDamageAndStressOnQuad(
     const Matrix<Real> & grad_u, Matrix<Real> & sigma, Real & dam,
     Real & Ehat) {
   if (!damage_in_compute_stress) {
     Vector<Real> Fdiag(3);
     Fdiag.zero();
 
     Matrix<Real> epsilon(3, 3);
     epsilon.zero();
-    for (UInt i = 0; i < spatial_dimension; ++i)
-      for (UInt j = 0; j < spatial_dimension; ++j)
+    for (UInt i = 0; i < spatial_dimension; ++i) {
+      for (UInt j = 0; j < spatial_dimension; ++j) {
         epsilon(i, j) = .5 * (grad_u(i, j) + grad_u(j, i));
-
+      }
+    }
     Math::matrixEig(3, epsilon.storage(), Fdiag.storage());
 
     computeDamageOnQuad(Ehat, sigma, Fdiag, dam);
   }
 
   sigma *= 1 - dam;
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt spatial_dimension, template <UInt> class Parent>
 inline void
 MaterialVonMisesMazars<spatial_dimension, Parent>::computeDamageOnQuad(
     const Real & epsilon_equ,
     __attribute__((unused)) const Matrix<Real> & sigma,
     const Vector<Real> & epsilon_princ, Real & dam) {
   Real Fs = epsilon_equ - K0;
   if (Fs > 0.) {
     Real dam_t;
     Real dam_c;
     dam_t =
         1 - K0 * (1 - At) / epsilon_equ - At * (exp(-Bt * (epsilon_equ - K0)));
     dam_c =
         1 - K0 * (1 - Ac) / epsilon_equ - Ac * (exp(-Bc * (epsilon_equ - K0)));
 
     Real Cdiag;
     Cdiag = this->E * (1 - this->nu) / ((1 + this->nu) * (1 - 2 * this->nu));
 
     Vector<Real> sigma_princ(3);
     sigma_princ(0) = Cdiag * epsilon_princ(0) +
                      this->lambda * (epsilon_princ(1) + epsilon_princ(2));
     sigma_princ(1) = Cdiag * epsilon_princ(1) +
                      this->lambda * (epsilon_princ(0) + epsilon_princ(2));
     sigma_princ(2) = Cdiag * epsilon_princ(2) +
                      this->lambda * (epsilon_princ(1) + epsilon_princ(0));
 
     Vector<Real> sigma_p(3);
-    for (UInt i = 0; i < 3; i++)
+    for (UInt i = 0; i < 3; i++) {
       sigma_p(i) = std::max(Real(0.), sigma_princ(i));
+    }
     // sigma_p *= 1. - dam;
 
     Real trace_p = this->nu / this->E * (sigma_p(0) + sigma_p(1) + sigma_p(2));
 
     Real alpha_t = 0;
     for (UInt i = 0; i < 3; ++i) {
       Real epsilon_t = (1 + this->nu) / this->E * sigma_p(i) - trace_p;
       Real epsilon_p = std::max(Real(0.), epsilon_princ(i));
       alpha_t += epsilon_t * epsilon_p;
     }
 
     alpha_t /= epsilon_equ * epsilon_equ;
     alpha_t = std::min(alpha_t, Real(1.));
 
     Real alpha_c = 1. - alpha_t;
 
     alpha_t = std::pow(alpha_t, beta);
     alpha_c = std::pow(alpha_c, beta);
 
     Real damtemp;
     damtemp = alpha_t * dam_t + alpha_c * dam_c;
 
     dam = std::max(damtemp, dam);
     dam = std::min(dam, Real(1.));
   }
 }
 
 } // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index 73c9c27c6..f35e1755b 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1242 +1,1242 @@
 /**
  * @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 memory_id the id of the memory
  * @param model_type this is an internal parameter for inheritance purposes
  */
 
 SolidMechanicsModel::SolidMechanicsModel(
     Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id,
     std::shared_ptr<DOFManager> dof_manager, const ModelType model_type)
-    : Model(mesh, model_type, dof_manager, dim, id, memory_id),
+    : Model(mesh, model_type, std::move(dof_manager), dim, id, memory_id),
       BoundaryCondition<SolidMechanicsModel>(),
       material_index("material index", id, memory_id),
       material_local_numbering("material local numbering", id, memory_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
   /*bool need_to_reassemble = true;*/
 
   //bool need_to_reassemble = false;
 
   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 std::string & energy_id) {
   AKANTU_DEBUG_IN();
 
   if (energy_id == "kinetic") {
     return getKineticEnergy();
   }
 
   if (energy_id == "external work") {
     return getExternalWork();
   }
 
   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;
 }
 
 /* -------------------------------------------------------------------------- */
 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(),
                                    this->getMemoryID());
 
   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, this->memory_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());
     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
diff --git a/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh
index 938ae3283..a8c75a40e 100644
--- a/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh
+++ b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh
@@ -1,450 +1,451 @@
 /**
  * @file   test_model_solver_my_model.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Apr 13 2016
  * @date last modification: Tue Feb 20 2018
  *
  * @brief  Test default dof manager
  *
  *
  * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_iterators.hh"
 #include "boundary_condition.hh"
 #include "communicator.hh"
 #include "data_accessor.hh"
 #include "dof_manager_default.hh"
 #include "element_synchronizer.hh"
 #include "mesh.hh"
 #include "model_solver.hh"
 #include "periodic_node_synchronizer.hh"
 #include "solver_vector_default.hh"
 #include "sparse_matrix.hh"
 
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 #ifndef AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH_
 #define AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH_
 
 /**
  *   =\o-----o-----o-> F
  *     |           |
  *     |---- L ----|
  */
 class MyModel : public ModelSolver,
                 public BoundaryCondition<MyModel>,
                 public DataAccessor<Element> {
 public:
   MyModel(Real F, Mesh & mesh, bool lumped,
           const ID & dof_manager_type = "default")
       : ModelSolver(mesh, ModelType::_model, "model_solver", 0),
         nb_dofs(mesh.getNbNodes()), nb_elements(mesh.getNbElement(_segment_2)),
         lumped(lumped), E(1.), A(1.), rho(1.), mesh(mesh),
         displacement(nb_dofs, 1, "disp"), velocity(nb_dofs, 1, "velo"),
         acceleration(nb_dofs, 1, "accel"), blocked(nb_dofs, 1, "blocked"),
         forces(nb_dofs, 1, "force_ext"),
         internal_forces(nb_dofs, 1, "force_int"),
         stresses(nb_elements, 1, "stress"), strains(nb_elements, 1, "strain"),
         initial_lengths(nb_elements, 1, "L0") {
-    
-    this->initBC(*this, displacement, forces);   
+
     this->initDOFManager(dof_manager_type);
 
+    this->initBC(*this, displacement, forces);
+
     this->getDOFManager().registerDOFs("disp", displacement, _dst_nodal);
     this->getDOFManager().registerDOFsDerivative("disp", 1, velocity);
     this->getDOFManager().registerDOFsDerivative("disp", 2, acceleration);
     this->getDOFManager().registerBlockedDOFs("disp", blocked);
 
     displacement.set(0.);
     velocity.set(0.);
     acceleration.set(0.);
 
     forces.set(0.);
     blocked.set(false);
 
     UInt global_nb_nodes = mesh.getNbGlobalNodes();
     for (auto && n : arange(nb_dofs)) {
       auto global_id = mesh.getNodeGlobalId(n);
       if (global_id == (global_nb_nodes - 1))
         forces(n, _x) = F;
 
       if (global_id == 0)
         blocked(n, _x) = true;
     }
 
     auto cit = this->mesh.getConnectivity(_segment_2).begin(2);
     auto cend = this->mesh.getConnectivity(_segment_2).end(2);
     auto L_it = this->initial_lengths.begin();
 
     for (; cit != cend; ++cit, ++L_it) {
       const Vector<UInt> & conn = *cit;
       UInt n1 = conn(0);
       UInt n2 = conn(1);
       Real p1 = this->mesh.getNodes()(n1, _x);
       Real p2 = this->mesh.getNodes()(n2, _x);
 
       *L_it = std::abs(p2 - p1);
     }
 
     this->registerDataAccessor(*this);
     this->registerSynchronizer(
         const_cast<ElementSynchronizer &>(this->mesh.getElementSynchronizer()),
         SynchronizationTag::_user_1);
   }
 
   void assembleLumpedMass() {
     auto & M = this->getDOFManager().getLumpedMatrix("M");
     M.zero();
 
     this->assembleLumpedMass(_not_ghost);
     if (this->mesh.getNbElement(_segment_2, _ghost) > 0)
       this->assembleLumpedMass(_ghost);
 
     is_lumped_mass_assembled = true;
   }
 
   void assembleLumpedMass(GhostType ghost_type) {
     Array<Real> M(nb_dofs, 1, 0.);
 
     Array<Real> m_all_el(this->mesh.getNbElement(_segment_2, ghost_type), 2);
 
     for (auto && data :
          zip(make_view(this->mesh.getConnectivity(_segment_2), 2),
              make_view(m_all_el, 2))) {
       const auto & conn = std::get<0>(data);
       auto & m_el = std::get<1>(data);
 
       UInt n1 = conn(0);
       UInt n2 = conn(1);
 
       Real p1 = this->mesh.getNodes()(n1, _x);
       Real p2 = this->mesh.getNodes()(n2, _x);
 
       Real L = std::abs(p2 - p1);
 
       Real M_n = rho * A * L / 2;
       m_el(0) = m_el(1) = M_n;
     }
 
     this->getDOFManager().assembleElementalArrayLocalArray(
         m_all_el, M, _segment_2, ghost_type);
 
     this->getDOFManager().assembleToLumpedMatrix("disp", M, "M");
   }
 
   void assembleMass() {
     SparseMatrix & M = this->getDOFManager().getMatrix("M");
     M.zero();
 
     Array<Real> m_all_el(this->nb_elements, 4);
 
     Matrix<Real> m(2, 2);
     m(0, 0) = m(1, 1) = 2;
     m(0, 1) = m(1, 0) = 1;
 
     // under integrated
     // m(0, 0) = m(1, 1) = 3./2.;
     // m(0, 1) = m(1, 0) = 3./2.;
 
     // lumping the mass matrix
     // m(0, 0) += m(0, 1);
     // m(1, 1) += m(1, 0);
     // m(0, 1) = m(1, 0) = 0;
 
     for (auto && data :
          zip(make_view(this->mesh.getConnectivity(_segment_2), 2),
              make_view(m_all_el, 2, 2))) {
       const auto & conn = std::get<0>(data);
       auto & m_el = std::get<1>(data);
       UInt n1 = conn(0);
       UInt n2 = conn(1);
 
       Real p1 = this->mesh.getNodes()(n1, _x);
       Real p2 = this->mesh.getNodes()(n2, _x);
 
       Real L = std::abs(p2 - p1);
 
       m_el = m;
       m_el *= rho * A * L / 6.;
     }
 
     this->getDOFManager().assembleElementalMatricesToMatrix(
         "M", "disp", m_all_el, _segment_2);
 
     is_mass_assembled = true;
   }
 
   MatrixType getMatrixType(const ID &) override { return _symmetric; }
 
   void assembleMatrix(const ID & matrix_id) override {
     if (matrix_id == "K") {
       if (not is_stiffness_assembled)
         this->assembleStiffness();
     } else if (matrix_id == "M") {
       if (not is_mass_assembled)
         this->assembleMass();
     } else if (matrix_id == "C") {
       // pass, no damping matrix
     } else {
       AKANTU_EXCEPTION("This solver does not know what to do with a matrix "
                        << matrix_id);
     }
   }
 
   void assembleLumpedMatrix(const ID & matrix_id) override {
     if (matrix_id == "M") {
       if (not is_lumped_mass_assembled)
         this->assembleLumpedMass();
     } else {
       AKANTU_EXCEPTION("This solver does not know what to do with a matrix "
                        << matrix_id);
     }
   }
 
   void assembleStiffness() {
     SparseMatrix & K = this->getDOFManager().getMatrix("K");
     K.zero();
 
     Matrix<Real> k(2, 2);
     k(0, 0) = k(1, 1) = 1;
     k(0, 1) = k(1, 0) = -1;
 
     Array<Real> k_all_el(this->nb_elements, 4);
 
     auto k_it = k_all_el.begin(2, 2);
     auto cit = this->mesh.getConnectivity(_segment_2).begin(2);
     auto cend = this->mesh.getConnectivity(_segment_2).end(2);
 
     for (; cit != cend; ++cit, ++k_it) {
       const auto & conn = *cit;
       UInt n1 = conn(0);
       UInt n2 = conn(1);
 
       Real p1 = this->mesh.getNodes()(n1, _x);
       Real p2 = this->mesh.getNodes()(n2, _x);
 
       Real L = std::abs(p2 - p1);
 
       auto & k_el = *k_it;
       k_el = k;
       k_el *= E * A / L;
     }
 
     this->getDOFManager().assembleElementalMatricesToMatrix(
         "K", "disp", k_all_el, _segment_2);
 
     is_stiffness_assembled = true;
   }
 
   void assembleResidual() override {
     this->getDOFManager().assembleToResidual("disp", forces);
     internal_forces.zero();
 
     this->assembleResidualInternal(_not_ghost);
 
     this->synchronize(SynchronizationTag::_user_1);
 
     this->getDOFManager().assembleToResidual("disp", internal_forces, -1.);
   }
 
   void assembleResidualInternal(GhostType ghost_type) {
     Array<Real> forces_internal_el(
         this->mesh.getNbElement(_segment_2, ghost_type), 2);
 
     auto cit = this->mesh.getConnectivity(_segment_2, ghost_type).begin(2);
     auto cend = this->mesh.getConnectivity(_segment_2, ghost_type).end(2);
 
     auto f_it = forces_internal_el.begin(2);
 
     auto strain_it = this->strains.begin();
     auto stress_it = this->stresses.begin();
     auto L_it = this->initial_lengths.begin();
 
     for (; cit != cend; ++cit, ++f_it, ++strain_it, ++stress_it, ++L_it) {
       const auto & conn = *cit;
       UInt n1 = conn(0);
       UInt n2 = conn(1);
 
       Real u1 = this->displacement(n1, _x);
       Real u2 = this->displacement(n2, _x);
 
       *strain_it = (u2 - u1) / *L_it;
       *stress_it = E * *strain_it;
       Real f_n = A * *stress_it;
       Vector<Real> & f = *f_it;
 
       f(0) = -f_n;
       f(1) = f_n;
     }
 
     this->getDOFManager().assembleElementalArrayLocalArray(
         forces_internal_el, internal_forces, _segment_2, ghost_type);
   }
 
   Real getPotentialEnergy() {
     Real res = 0;
 
     if (not lumped) {
       res = this->mulVectMatVect(this->displacement, "K", this->displacement);
     } else {
       auto strain_it = this->strains.begin();
       auto stress_it = this->stresses.begin();
       auto strain_end = this->strains.end();
       auto L_it = this->initial_lengths.begin();
 
       for (; strain_it != strain_end; ++strain_it, ++stress_it, ++L_it) {
         res += *strain_it * *stress_it * A * *L_it;
       }
 
       mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum);
     }
 
     return res / 2.;
   }
 
   Real getKineticEnergy() {
     Real res = 0;
     if (not lumped) {
       res = this->mulVectMatVect(this->velocity, "M", this->velocity);
     } else {
       Array<Real> & m = dynamic_cast<SolverVectorDefault &>(
           this->getDOFManager().getLumpedMatrix("M"));
       auto it = velocity.begin();
       auto end = velocity.end();
       auto m_it = m.begin();
 
       for (UInt node = 0; it != end; ++it, ++m_it, ++node) {
         if (mesh.isLocalOrMasterNode(node))
           res += *m_it * *it * *it;
       }
 
       mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum);
     }
 
     return res / 2.;
   }
 
   Real getExternalWorkIncrement() {
     Real res = 0;
 
     auto it = velocity.begin();
     auto end = velocity.end();
     auto if_it = internal_forces.begin();
     auto ef_it = forces.begin();
     auto b_it = blocked.begin();
 
     for (UInt node = 0; it != end; ++it, ++if_it, ++ef_it, ++b_it, ++node) {
       if (mesh.isLocalOrMasterNode(node))
         res += (*b_it ? -*if_it : *ef_it) * *it;
     }
 
     mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum);
 
     return res * this->getTimeStep();
   }
 
   Real mulVectMatVect(const Array<Real> & x, const ID & A_id,
                       const Array<Real> & y) {
     Array<Real> Ay(nb_dofs);
     this->getDOFManager().assembleMatMulVectToArray("disp", A_id, y, Ay);
 
     Real res = 0.;
     for (auto && data : zip(arange(nb_dofs), make_view(Ay), make_view(x))) {
       res += std::get<2>(data) * std::get<1>(data) *
              mesh.isLocalOrMasterNode(std::get<0>(data));
     }
 
     mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum);
     return res;
   }
 
   /* ------------------------------------------------------------------------ */
   UInt getNbData(const Array<Element> & elements,
                  const SynchronizationTag &) const override {
     return elements.size() * sizeof(Real);
   }
 
   void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
                 const SynchronizationTag & tag) const override {
     if (tag == SynchronizationTag::_user_1) {
       for (const auto & el : elements) {
         buffer << this->stresses(el.element);
       }
     }
   }
 
   void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
                   const SynchronizationTag & tag) override {
     if (tag == SynchronizationTag::_user_1) {
       auto cit = this->mesh.getConnectivity(_segment_2, _ghost).begin(2);
 
       for (const auto & el : elements) {
         Real stress;
         buffer >> stress;
 
         Real f = A * stress;
 
         Vector<UInt> conn = cit[el.element];
         this->internal_forces(conn(0), _x) += -f;
         this->internal_forces(conn(1), _x) += f;
       }
     }
   }
 
   const Mesh & getMesh() const { return mesh; }
 
   UInt getSpatialDimension() const { return 1; }
 
   auto & getBlockedDOFs() { return blocked; }
 
 private:
   UInt nb_dofs;
   UInt nb_elements;
 
   bool lumped;
 
   bool is_stiffness_assembled{false};
   bool is_mass_assembled{false};
   bool is_lumped_mass_assembled{false};
 
 public:
   Real E, A, rho;
 
   Mesh & mesh;
   Array<Real> displacement;
   Array<Real> velocity;
   Array<Real> acceleration;
   Array<bool> blocked;
   Array<Real> forces;
   Array<Real> internal_forces;
 
   Array<Real> stresses;
   Array<Real> strains;
 
   Array<Real> initial_lengths;
 };
 
 #endif /* AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH_ */
 
 } // namespace akantu