diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc
index 31bcdcf75..1a073da25 100644
--- a/python/py_dof_manager.cc
+++ b/python/py_dof_manager.cc
@@ -1,215 +1,216 @@
 /*
  * Copyright (©) 2018-2021 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 "py_dof_manager.hh"
 #include "py_aka_array.hh"
 #include "py_akantu_pybind11_compatibility.hh"
 /* -------------------------------------------------------------------------- */
 #include <dof_manager.hh>
 #include <non_linear_solver.hh>
 #include <solver_callback.hh>
 #include <time_step_solver.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 namespace {
   class PySolverCallback : public SolverCallback {
   public:
     using SolverCallback::SolverCallback;
 
     /// get the type of matrix needed
     MatrixType getMatrixType(const ID & matrix_id) const override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE_PURE(MatrixType, SolverCallback, getMatrixType,
                              matrix_id);
     }
 
     /// callback to assemble a Matrix
     void assembleMatrix(const ID & matrix_id) override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleMatrix, matrix_id);
     }
 
     /// callback to assemble a lumped Matrix
     void assembleLumpedMatrix(const ID & matrix_id) override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleLumpedMatrix,
                              matrix_id);
     }
 
     /// callback to assemble the residual (rhs)
     void assembleResidual() override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleResidual);
     }
 
     /// callback for the predictor (in case of dynamic simulation)
     void predictor() override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, SolverCallback, predictor);
     }
 
     /// callback for the corrector (in case of dynamic simulation)
     void corrector() override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, SolverCallback, corrector);
     }
 
     void beforeSolveStep() override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, SolverCallback, beforeSolveStep);
     }
 
     void afterSolveStep(bool converged) override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, SolverCallback, afterSolveStep, converged);
     }
   };
 
   class PyInterceptSolverCallback : public InterceptSolverCallback {
   public:
     using InterceptSolverCallback::InterceptSolverCallback;
 
     MatrixType getMatrixType(const ID & matrix_id) const override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(MatrixType, InterceptSolverCallback, getMatrixType,
                         matrix_id);
     }
 
     void assembleMatrix(const ID & matrix_id) override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleMatrix,
                         matrix_id);
     }
 
     /// callback to assemble a lumped Matrix
     void assembleLumpedMatrix(const ID & matrix_id) override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleLumpedMatrix,
                         matrix_id);
     }
 
     void assembleResidual() override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleResidual);
     }
 
     void predictor() override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, InterceptSolverCallback, predictor);
     }
 
     void corrector() override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, InterceptSolverCallback, corrector);
     }
 
     void beforeSolveStep() override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, InterceptSolverCallback, beforeSolveStep);
     }
 
     void afterSolveStep(bool converged) override {
       // NOLINTNEXTLINE
       PYBIND11_OVERRIDE(void, InterceptSolverCallback, afterSolveStep,
                         converged);
     }
   };
 
 } // namespace
 
 /* -------------------------------------------------------------------------- */
 void register_dof_manager(py::module & mod) {
   py::class_<DOFManager, std::shared_ptr<DOFManager>>(mod, "DOFManager")
       .def("getMatrix", &DOFManager::getMatrix,
            py::return_value_policy::reference)
       .def(
           "getNewMatrix",
           [](DOFManager & self, const std::string & name,
              const std::string & matrix_to_copy_id) -> decltype(auto) {
             return self.getNewMatrix(name, matrix_to_copy_id);
           },
           py::return_value_policy::reference)
       .def(
           "getResidual",
           [](DOFManager & self) -> decltype(auto) {
             return self.getResidual();
           },
           py::return_value_policy::reference)
       .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs)
       .def(
           "hasMatrix",
           [](DOFManager & self, const ID & name) -> bool {
             return self.hasMatrix(name);
           },
           py::arg("name"))
       .def("assembleToResidual", &DOFManager::assembleToResidual,
            py::arg("dof_id"), py::arg("array_to_assemble"),
            py::arg("scale_factor") = 1.)
       .def("assembleToLumpedMatrix", &DOFManager::assembleToLumpedMatrix,
            py::arg("dof_id"), py::arg("array_to_assemble"),
            py::arg("lumped_mtx"), py::arg("scale_factor") = 1.)
       .def("assemblePreassembledMatrix",
            &DOFManager::assemblePreassembledMatrix, py::arg("matrix_id"),
-           py::arg("terms"));
+           py::arg("terms"))
+      .def("zeroResidual", &DOFManager::zeroResidual);
 
   py::class_<NonLinearSolver>(mod, "NonLinearSolver")
       .def(
           "set",
           [](NonLinearSolver & self, const std::string & id, const Real & val) {
             if (id == "max_iterations") {
               self.set(id, int(val));
             } else {
               self.set(id, val);
             }
           })
       .def("set",
            [](NonLinearSolver & self, const std::string & id,
               const SolveConvergenceCriteria & val) { self.set(id, val); });
 
   py::class_<TimeStepSolver>(mod, "TimeStepSolver")
       .def("getIntegrationScheme", &TimeStepSolver::getIntegrationScheme);
 
   py::class_<SolverCallback, PySolverCallback>(mod, "SolverCallback")
       .def(py::init_alias<DOFManager &>())
       .def("getMatrixType", &SolverCallback::getMatrixType)
       .def("assembleMatrix", &SolverCallback::assembleMatrix)
       .def("assembleLumpedMatrix", &SolverCallback::assembleLumpedMatrix)
       .def("assembleResidual",
            [](SolverCallback & self) { self.assembleResidual(); })
       .def("predictor", &SolverCallback::predictor)
       .def("corrector", &SolverCallback::corrector)
       .def("beforeSolveStep", &SolverCallback::beforeSolveStep)
       .def("afterSolveStep", &SolverCallback::afterSolveStep)
       .def_property_readonly("dof_manager", &SolverCallback::getSCDOFManager,
                              py::return_value_policy::reference);
 
   py::class_<InterceptSolverCallback, SolverCallback,
              PyInterceptSolverCallback>(mod, "InterceptSolverCallback")
       .def(py::init_alias<SolverCallback &>());
 }
 
 } // namespace akantu
diff --git a/python/py_solver.cc b/python/py_solver.cc
index bbcd931c3..d5e20670d 100644
--- a/python/py_solver.cc
+++ b/python/py_solver.cc
@@ -1,113 +1,119 @@
 /**
  * @file   py_solver.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Sep 29 2020
  * @date last modification: Sat Mar 06 2021
  *
  * @brief  pybind11 interface to Solver and SparseMatrix
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-2021 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 "py_solver.hh"
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <model.hh>
 #include <non_linear_solver.hh>
 #include <sparse_matrix_aij.hh>
 #include <terms_to_assemble.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 void register_solvers(py::module & mod) {
   py::class_<SparseMatrix>(mod, "SparseMatrix")
       .def("getMatrixType", &SparseMatrix::getMatrixType)
       .def("size", &SparseMatrix::size)
       .def("zero", &SparseMatrix::zero)
       .def("saveProfile", &SparseMatrix::saveProfile)
       .def("saveMatrix", &SparseMatrix::saveMatrix)
       .def(
           "add", [](SparseMatrix & self, UInt i, UInt j) { self.add(i, j); },
           "Add entry in the profile")
       .def(
           "add",
           [](SparseMatrix & self, UInt i, UInt j, Real value) {
             self.add(i, j, value);
           },
           "Add the value to the matrix")
       .def(
           "add",
           [](SparseMatrix & self, SparseMatrix & A, Real alpha) {
             self.add(A, alpha);
           },
           "Add a matrix to the matrix", py::arg("A"), py::arg("alpha") = 1.)
+
+      .def("isFinite", &SparseMatrix::isFinite)
+
+      .def("getRelease",
+           [](const SparseMatrix & self) -> UInt { return self.getRelease(); })
       .def("__call__",
            [](const SparseMatrix & self, UInt i, UInt j) { return self(i, j); })
       .def("getRelease", &SparseMatrix::getRelease);
 
   py::class_<SparseMatrixAIJ, SparseMatrix>(mod, "SparseMatrixAIJ")
       .def("getIRN", &SparseMatrixAIJ::getIRN)
       .def("getJCN", &SparseMatrixAIJ::getJCN)
       .def("getA", &SparseMatrixAIJ::getA);
 
   py::class_<SolverVector>(mod, "SolverVector")
-        .def("getValues",
-            [](SolverVector& self) -> decltype(auto) {
-        	return static_cast<const Array<Real>& >(self);
-            },
-	    py::return_value_policy::reference_internal,
-	    "Transform this into a vector, Is not copied.")
-	.def("isDistributed", [](const SolverVector& self) { return self.isDistributed(); })
-  ;
+      .def(
+          "getValues",
+          [](SolverVector & self) -> decltype(auto) {
+            return static_cast<const Array<Real> &>(self);
+          },
+          py::return_value_policy::reference_internal,
+          "Transform this into a vector, Is not copied.")
+      .def("isDistributed",
+           [](const SolverVector & self) { return self.isDistributed(); });
 
   py::class_<TermsToAssemble::TermToAssemble>(mod, "TermToAssemble")
       .def(py::init<Int, Int>())
       .def(py::self += Real())
       .def_property_readonly("i", &TermsToAssemble::TermToAssemble::i)
       .def_property_readonly("j", &TermsToAssemble::TermToAssemble::j);
 
   py::class_<TermsToAssemble>(mod, "TermsToAssemble")
       .def(py::init<const ID &, const ID &>())
       .def("getDOFIdM", &TermsToAssemble::getDOFIdM)
       .def("getDOFIdN", &TermsToAssemble::getDOFIdN)
       .def(
           "__call__",
           [](TermsToAssemble & self, UInt i, UInt j, Real val) {
             auto & term = self(i, j);
             term = val;
             return term;
           },
           py::arg("i"), py::arg("j"), py::arg("val") = 0.,
           py::return_value_policy::reference);
 }
 
 } // namespace akantu
diff --git a/python/py_structural_mechanics_model.cc b/python/py_structural_mechanics_model.cc
index a80b6a529..366953881 100644
--- a/python/py_structural_mechanics_model.cc
+++ b/python/py_structural_mechanics_model.cc
@@ -1,162 +1,179 @@
 /**
  * @file   py_structural_mechanics_model.cc
  *
  * @author Philip Mueller <philip.paul.mueller@bluemail.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Feb 03 2021
  * @date last modification: Thu Apr 01 2021
  *
  * @brief  pybind11 interface to StructuralMechanicsModel
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-2021 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 "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <structural_mechanics_model.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 #define def_deprecated(func_name, mesg)                                        \
   def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); })
 
 #define def_function_nocopy(func_name)                                         \
   def(                                                                         \
       #func_name,                                                              \
       [](StructuralMechanicsModel & self) -> decltype(auto) {                  \
         return self.func_name();                                               \
       },                                                                       \
       py::return_value_policy::reference)
 
 #define def_function_(func_name)                                               \
   def(#func_name, [](StructuralMechanicsModel & self) -> decltype(auto) {      \
     return self.func_name();                                                   \
   })
 
 #define def_plainmember(M) def_readwrite(#M, &StructuralMaterial::M)
 /* -------------------------------------------------------------------------- */
 
 void register_structural_mechanics_model(pybind11::module & mod) {
   /* First we have to register the material class
    * The wrapper aims to mimic the behaviour of the real material.
    */
   py::class_<StructuralMaterial>(mod, "StructuralMaterial")
       .def(py::init<>())
       .def(py::init<const StructuralMaterial &>())
       .def_plainmember(E)
       .def_plainmember(A)
       .def_plainmember(I)
       .def_plainmember(Iz)
       .def_plainmember(Iy)
       .def_plainmember(GJ)
       .def_plainmember(rho)
       .def_plainmember(t)
       .def_plainmember(nu);
 
   /* Now we create the structural model wrapper
    * Note that this is basically a port from the solid mechanic part.
    */
   py::class_<StructuralMechanicsModel, Model>(mod, "StructuralMechanicsModel")
       .def(py::init<Mesh &, UInt, const ID &>(), py::arg("mesh"),
            py::arg("spatial_dimension") = _all_dimensions,
            py::arg("id") = "structural_mechanics_model")
       .def(
           "initFull",
           [](StructuralMechanicsModel & self,
              const AnalysisMethod & analysis_method) -> void {
             self.initFull(_analysis_method = analysis_method);
           },
           py::arg("_analysis_method"))
       .def("initFull",
            [](StructuralMechanicsModel & self) -> void { self.initFull(); })
       .def_function_nocopy(getExternalForce)
       .def_function_nocopy(getDisplacement)
       .def_function_nocopy(getInternalForce)
       .def_function_nocopy(getVelocity)
       .def_function_nocopy(getAcceleration)
       .def_function_nocopy(getInternalForce)
       .def_function_nocopy(getBlockedDOFs)
       .def_function_nocopy(getMesh)
 
       .def("setTimeStep", &StructuralMechanicsModel::setTimeStep,
            py::arg("time_step"), py::arg("solver_id") = "")
       .def(
           "getElementMaterial",
           [](StructuralMechanicsModel & self, const ElementType & type,
              GhostType ghost_type) -> decltype(auto) {
             return self.getElementMaterial(type, ghost_type);
           },
           "This function returns the map that maps elements to materials.",
           py::arg("type"), py::arg("ghost_type") = _not_ghost,
           py::return_value_policy::reference)
       .def(
           "getMaterialByElement",
           [](StructuralMechanicsModel & self, Element element)
               -> decltype(auto) { return self.getMaterialByElement(element); },
           "This function returns the `StructuralMaterial` instance that is "
           "associated with element `element`.",
           py::arg("element"), py::return_value_policy::reference)
       .def(
           "addMaterial",
           [](StructuralMechanicsModel & self, StructuralMaterial & mat,
              const ID & name) -> UInt { return self.addMaterial(mat, name); },
           "This function adds the `StructuralMaterial` `mat` to `self`."
           " The function returns the ID of the new material.",
           py::arg("mat"), py::arg("name") = "")
       .def(
           "getMaterial",
-          [](StructuralMechanicsModel & self, UInt material_index)
-              -> decltype(auto) { return self.getMaterial(material_index); },
+          [](const StructuralMechanicsModel & self,
+             UInt material_index) -> StructuralMaterial {
+            return self.getMaterial(material_index);
+          },
           "This function returns the `i`th material of `self`."
-          " Note a reference is returned which allows to modify the material "
-          "inside `self`.",
-          py::arg("i"), py::return_value_policy::reference)
+          " The function returns a copy of the material.",
+          py::arg("material_index"), py::return_value_policy::copy)
       .def(
           "getMaterial",
-          [](StructuralMechanicsModel & self, const ID & name)
-              -> decltype(auto) { return self.getMaterial(name); },
-          "This function returns the material with name `i` of `self`."
-          " Note a reference is returned which allows to modify the material "
-          "inside `self`.",
-          py::arg("i"), py::return_value_policy::reference)
+          [](const StructuralMechanicsModel & self, const ID & name)
+              -> StructuralMaterial { return self.getMaterial(name); },
+          "This function returns the `i`th material of `self`."
+          " The function returns a copy of the material.",
+          py::arg("material_index"), py::return_value_policy::copy)
       .def(
           "getNbMaterials",
           [](StructuralMechanicsModel & self) { return self.getNbMaterials(); },
           "Returns the number of different materials inside `self`.")
+
       .def("getKineticEnergy", &StructuralMechanicsModel::getKineticEnergy,
            "Compute kinetic energy")
       .def("getPotentialEnergy", &StructuralMechanicsModel::getPotentialEnergy,
            "Compute potential energy")
       .def("getEnergy", &StructuralMechanicsModel::getEnergy,
-           "Compute the specified energy");
+           "Compute the specified energy")
 
-} // End: register structural mechanical model
+      .def(
+          "getLumpedMass",
+          [](const StructuralMechanicsModel & self) -> decltype(auto) {
+            return self.getLumpedMass();
+          },
+          py::return_value_policy::reference_internal)
+      .def(
+          "getMass",
+          [](const StructuralMechanicsModel & self) -> decltype(auto) {
+            return self.getMass();
+          },
+          py::return_value_policy::reference_internal)
+      .def("assembleLumpedMassMatrix",
+           &StructuralMechanicsModel::assembleLumpedMassMatrix,
+           "Assembles the lumped mass matrix")
+      .def("hasLumpedMass", &StructuralMechanicsModel::hasLumpedMass);
+}
 
 } // namespace akantu
diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh
index 3a346fc5d..21d980519 100644
--- a/src/common/aka_array.hh
+++ b/src/common/aka_array.hh
@@ -1,445 +1,450 @@
 /**
  * @file   aka_array.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Sun Nov 22 2020
  *
  * @brief  Array container for Akantu This container differs from the
  * std::vector from the fact it as 2 dimensions a main dimension and the size
  * stored per entries
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 #include "aka_types.hh"
 /* -------------------------------------------------------------------------- */
 #include <typeinfo>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #ifndef AKANTU_ARRAY_HH_
 #define AKANTU_ARRAY_HH_
 
 namespace akantu {
 
 /// class that afford to store vectors in static memory
 // NOLINTNEXTLINE(cppcoreguidelines-special-member-functions)
 class ArrayBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   explicit ArrayBase(const ID & id = "") : id(id) {}
   ArrayBase(const ArrayBase & other, const ID & id = "") {
     this->id = (id.empty()) ? other.id : id;
   }
 
   ArrayBase(ArrayBase && other) = default;
   ArrayBase & operator=(const ArrayBase & other) = default;
   ArrayBase & operator=(ArrayBase && other) noexcept = default;
 
   virtual ~ArrayBase() = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the amount of space allocated in bytes
   virtual UInt getMemorySize() const = 0;
 
   // changed empty to match std::vector empty
   inline bool empty() const __attribute__((warn_unused_result)) {
     return size_ == 0;
   }
 
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors */
   /* ------------------------------------------------------------------------ */
 public:
   /// Get the Size of the Array
   UInt size() const { return size_; }
   /// Get the number of components
   AKANTU_GET_MACRO(NbComponent, nb_component, UInt);
   /// Get the name of th array
   AKANTU_GET_MACRO(ID, id, const ID &);
   /// Set the name of th array
   AKANTU_SET_MACRO(ID, id, const ID &);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// id of the vector
   ID id;
 
   /// the size used
   UInt size_{0};
 
   /// number of components
   UInt nb_component{1};
 };
 
 /* -------------------------------------------------------------------------- */
 namespace {
   template <std::size_t dim, typename T> struct IteratorHelper {};
 
   template <typename T> struct IteratorHelper<0, T> { using type = T; };
   template <typename T> struct IteratorHelper<1, T> { using type = Vector<T>; };
   template <typename T> struct IteratorHelper<2, T> { using type = Matrix<T>; };
   template <typename T> struct IteratorHelper<3, T> {
     using type = Tensor3<T>;
   };
 
   template <std::size_t dim, typename T>
   using IteratorHelper_t = typename IteratorHelper<dim, T>::type;
 } // namespace
 
 /* -------------------------------------------------------------------------- */
 /* Memory handling layer                                                      */
 /* -------------------------------------------------------------------------- */
 enum class ArrayAllocationType {
   _default,
   _pod,
 };
 
 template <typename T>
 struct ArrayAllocationTrait
     : public std::conditional_t<
           std::is_scalar<T>::value,
           std::integral_constant<ArrayAllocationType,
                                  ArrayAllocationType::_pod>,
           std::integral_constant<ArrayAllocationType,
                                  ArrayAllocationType::_default>> {};
 
 /* -------------------------------------------------------------------------- */
 template <typename T,
           ArrayAllocationType allocation_trait = ArrayAllocationTrait<T>::value>
 class ArrayDataLayer : 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 = default;
 
   /// Allocation of a new vector
   explicit ArrayDataLayer(UInt size = 0, UInt nb_component = 1,
                           const ID & id = "");
 
   /// Allocation of a new vector with a default value
   ArrayDataLayer(UInt size, UInt nb_component, const_reference value,
                  const ID & id = "");
 
   /// Copy constructor (deep copy)
   ArrayDataLayer(const ArrayDataLayer & vect, const ID & id = "");
 
   /// Copy constructor (deep copy)
   explicit ArrayDataLayer(const std::vector<value_type> & vect);
 
   // copy operator
   ArrayDataLayer & operator=(const ArrayDataLayer & other);
 
   // move constructor
   ArrayDataLayer(ArrayDataLayer && other) noexcept = default;
 
   // move assign
   ArrayDataLayer & operator=(ArrayDataLayer && other) noexcept = default;
 
 protected:
   // deallocate the memory
   virtual void deallocate() {}
 
   // allocate the memory
   virtual void allocate(UInt size, UInt nb_component);
 
   // allocate and initialize the memory
   virtual void allocate(UInt size, UInt nb_component, const T & value);
 
 public:
   /// append a tuple of size nb_component containing value
   inline void push_back(const_reference value);
   /// append a vector
   // inline void push_back(const value_type new_elem[]);
 
   /// 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);
 
   /// changes the allocated size but not the size, if new_size = 0, the size is
   /// set to min(current_size and reserve size)
   virtual void reserve(UInt size, UInt new_size = UInt(-1));
 
   /// change the size of the Array
   virtual void resize(UInt size);
 
   /// change the size of the Array and initialize the values
   virtual void resize(UInt size, const T & val);
 
   /// get the amount of space allocated in bytes
   inline UInt getMemorySize() const override;
 
   /// Get the real size allocated in memory
   inline UInt getAllocatedSize() const;
 
   /// give the address of the memory allocated for this vector
   T * storage() const { return values; };
 
 protected:
   /// allocation type agnostic  data access
   T * values{nullptr};
 
   /// data storage
   std::vector<T> data_storage;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Actual Array                                                               */
 /* -------------------------------------------------------------------------- */
 template <typename T, bool is_scal> class Array : public ArrayDataLayer<T> {
 private:
   using parent = ArrayDataLayer<T>;
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   using value_type = typename parent::value_type;
   using reference = typename parent::reference;
   using pointer_type = typename parent::pointer_type;
   using const_reference = typename parent::const_reference;
   using array_type = Array<T>;
 
   ~Array() override;
 
   Array() : Array(0){};
 
   /// Allocation of a new vector
   explicit Array(UInt size, UInt nb_component = 1, const ID & id = "");
 
   /// Allocation of a new vector with a default value
   explicit Array(UInt size, UInt nb_component, const_reference value,
                  const ID & id = "");
 
   /// Copy constructor
   Array(const Array & vect, const ID & id = "");
 
   /// Copy constructor (deep copy)
   explicit Array(const std::vector<T> & vect);
 
   // copy operator
   Array & operator=(const Array & other);
 
   // move constructor
   Array(Array && other) noexcept = default;
 
   // move assign
   Array & operator=(Array && other) noexcept = default;
 
   /* ------------------------------------------------------------------------ */
   /* Iterator                                                                 */
   /* ------------------------------------------------------------------------ */
   /// \todo protected: does not compile with intel  check why
 public:
   template <class R, class it, class IR = R,
             bool is_tensor_ = aka::is_tensor<std::decay_t<R>>::value>
   class iterator_internal;
 
 public:
   /* ------------------------------------------------------------------------ */
 
   /* ------------------------------------------------------------------------ */
   template <typename R = T> class const_iterator;
   template <typename R = T> class iterator;
 
   /* ------------------------------------------------------------------------ */
 
   /// iterator for Array of nb_component = 1
   using scalar_iterator = iterator<T>;
   /// const_iterator for Array of nb_component = 1
   using const_scalar_iterator = const_iterator<T>;
 
   /// iterator returning Vectors of size n  on entries of Array with
   /// nb_component = n
   using vector_iterator = iterator<Vector<T>>;
   /// const_iterator returning Vectors of n size on entries of Array with
   /// nb_component = n
   using const_vector_iterator = const_iterator<Vector<T>>;
 
   /// iterator returning Matrices of size (m, n) on entries of Array with
   /// nb_component = m*n
   using matrix_iterator = iterator<Matrix<T>>;
   /// const iterator returning Matrices of size (m, n) on entries of Array with
   /// nb_component = m*n
   using const_matrix_iterator = const_iterator<Matrix<T>>;
 
   /// iterator returning Tensor3 of size (m, n, k) on entries of Array with
   /// nb_component = m*n*k
   using tensor3_iterator = iterator<Tensor3<T>>;
   /// const iterator returning Tensor3 of size (m, n, k) on entries of Array
   /// with nb_component = m*n*k
   using const_tensor3_iterator = const_iterator<Tensor3<T>>;
 
   /* ------------------------------------------------------------------------ */
   template <typename... Ns> inline decltype(auto) begin(Ns &&... n);
   template <typename... Ns> inline decltype(auto) end(Ns &&... n);
 
   template <typename... Ns> inline decltype(auto) begin(Ns &&... n) const;
   template <typename... Ns> inline decltype(auto) end(Ns &&... n) const;
 
   template <typename... Ns> inline decltype(auto) begin_reinterpret(Ns &&... n);
   template <typename... Ns> inline decltype(auto) end_reinterpret(Ns &&... n);
 
   template <typename... Ns>
   inline decltype(auto) begin_reinterpret(Ns &&... n) const;
   template <typename... Ns>
   inline decltype(auto) end_reinterpret(Ns &&... n) const;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// search elem in the vector, return  the position of the first occurrence or
   /// -1 if not found
   UInt find(const_reference elem) const;
 
   /// @see Array::find(const_reference elem) const
   //  UInt find(T elem[]) const;
 
   /// append a value to the end of the Array
   inline void push_back(const_reference value) { parent::push_back(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) {
     parent::push_back(new_elem);
   }
 
   /// append the content of the iterator at the end of the Array
   template <typename Ret> inline void push_back(const iterator<Ret> & it) {
     push_back(*it);
   }
 
   /// erase the value at position i
   inline void erase(UInt i);
   /// ask Nico, clarify
   template <typename R> inline iterator<R> erase(const iterator<R> & it);
 
   /// @see Array::find(const_reference elem) const
   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 UInt find(const C<T> & elem);
 
   /// set all entries of the array to the value t
   /// @param t value to fill the array with
   inline void set(T t) {
     std::fill_n(this->values, this->size_ * this->nb_component, t);
   }
 
   /// set the array to T{}
   inline void zero() { this->set({}); }
 
   /// resize the array to 0
   inline void clear() { this->resize(0); }
 
   /// set all tuples of the array to a given vector or matrix
   /// @param vm Matrix or Vector to fill the array with
   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 set(const C<T> & vm);
 
   /// Append the content of the other array to the current one
   void append(const Array<T> & other);
 
   /// copy another Array in the current Array, the no_sanity_check allows you to
   /// force the copy in cases where you know what you do with two non matching
   /// Arrays in terms of n
   void copy(const Array<T, is_scal> & other, bool no_sanity_check = false);
 
   /// function to print the containt of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
+  /// Tests if all elements are finite.
+  template <typename OT = T,
+            std::enable_if_t<std::is_arithmetic<OT>::value> * = nullptr>
+  bool isFinite() const noexcept;
+
   /* ------------------------------------------------------------------------ */
   /* Operators                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// substraction entry-wise
   Array<T, is_scal> & operator-=(const Array<T, is_scal> & other);
   /// addition entry-wise
   Array<T, is_scal> & operator+=(const Array<T, is_scal> & other);
   /// multiply evry entry by alpha
   Array<T, is_scal> & operator*=(const T & alpha);
 
   /// check if the array are identical entry-wise
   bool operator==(const Array<T, is_scal> & other) const;
   /// @see Array::operator==(const Array<T, is_scal> & other) const
   bool operator!=(const Array<T, is_scal> & other) const;
 
   /// return a reference to the j-th entry of the i-th tuple
   inline reference operator()(UInt i, UInt j = 0);
   /// return a const reference to the j-th entry of the i-th tuple
   inline const_reference operator()(UInt i, UInt j = 0) const;
 
   /// return a reference to the ith component of the 1D array
   inline reference operator[](UInt i);
   /// return a const reference to the ith component of the 1D array
   inline const_reference operator[](UInt i) const;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline Functions Array<T, is_scal>                                         */
 /* -------------------------------------------------------------------------- */
 template <typename T, bool is_scal>
 inline std::ostream & operator<<(std::ostream & stream,
                                  const Array<T, is_scal> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Inline Functions ArrayBase                                                 */
 /* -------------------------------------------------------------------------- */
 inline std::ostream & operator<<(std::ostream & stream,
                                  const ArrayBase & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "aka_array_tmpl.hh"
 
 #endif /* AKANTU_ARRAY_HH_ */
diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh
index 6813a95d7..664bbba09 100644
--- a/src/common/aka_array_tmpl.hh
+++ b/src/common/aka_array_tmpl.hh
@@ -1,1363 +1,1372 @@
 /**
  * @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: Fri Feb 26 2021
  *
  * @brief  Inline functions of the classes Array<T> and ArrayBase
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 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 <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) {
         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;
 }
 
-/* --------------------------------------------------------------------------
- */
+/* -------------------------------------------------------------------------- */
+template <typename T, bool is_scal>
+template <typename OT, std::enable_if_t<std::is_arithmetic<OT>::value> *>
+bool Array<T, is_scal>::isFinite() const noexcept {
+  return std::all_of(this->values,
+                     this->values + this->size_ * this->nb_component,
+                     [](auto && a) { return std::isfinite(a); });
+}
+
+/* -------------------------------------------------------------------------- */
 /* 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 daughter operator-(difference_type n) const {
+    return daughter(ret - n);
+  }
+  inline daughter operator+(difference_type n) const {
+    return daughter(ret + n);
+  }
 
-  inline difference_type operator-(const iterator_internal & b) {
+  inline difference_type operator-(const iterator_internal & b) const {
     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;
 
 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));
+  inline daughter operator+(difference_type n) const {
+    daughter tmp(static_cast<const daughter &>(*this));
     tmp += n;
     return tmp;
   }
-  inline daughter operator-(difference_type n) {
-    daughter tmp(static_cast<daughter &>(*this));
+  inline daughter operator-(difference_type n) const {
+    daughter tmp(static_cast<const daughter &>(*this));
     tmp -= n;
     return tmp;
   }
 
-  inline difference_type operator-(const iterator_internal & b) {
+  inline difference_type operator-(const iterator_internal & b) const {
     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/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh
index 545629934..296ca2117 100644
--- a/src/fe_engine/fe_engine_template.hh
+++ b/src/fe_engine/fe_engine_template.hh
@@ -1,431 +1,431 @@
 /**
  * @file   fe_engine_template.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Fri May 14 2021
  *
  * @brief  templated class that calls integration and shape objects
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "fe_engine.hh"
 #include "integrator.hh"
 #include "shape_functions.hh"
 /* -------------------------------------------------------------------------- */
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #ifndef AKANTU_FE_ENGINE_TEMPLATE_HH_
 #define AKANTU_FE_ENGINE_TEMPLATE_HH_
 
 namespace akantu {
 class DOFManager;
 namespace fe_engine {
   namespace details {
     template <ElementKind> struct AssembleLumpedTemplateHelper;
     template <ElementKind> struct AssembleFieldMatrixHelper;
   } // namespace details
 } // namespace fe_engine
 
 template <ElementKind, typename> struct AssembleFieldMatrixStructHelper;
 
 struct DefaultIntegrationOrderFunctor {
   template <ElementType type> static inline constexpr int getOrder() {
     return ElementClassProperty<type>::polynomial_degree;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <template <ElementKind, class> class I, template <ElementKind> class S,
           ElementKind kind = _ek_regular,
           class IntegrationOrderFunctor = DefaultIntegrationOrderFunctor>
 class FEEngineTemplate : public FEEngine {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   using Integ = I<kind, IntegrationOrderFunctor>;
   using Shape = S<kind>;
 
   FEEngineTemplate(Mesh & mesh, UInt spatial_dimension = _all_dimensions,
                    const ID & id = "fem");
 
   ~FEEngineTemplate() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// pre-compute all the shape functions, their derivatives and the jacobians
   void initShapeFunctions(GhostType ghost_type = _not_ghost) override;
   void initShapeFunctions(const Array<Real> & nodes,
                           GhostType ghost_type = _not_ghost);
 
   /* ------------------------------------------------------------------------ */
   /* Integration method bridges                                               */
   /* ------------------------------------------------------------------------ */
   /// integrate f for all elements of type "type"
   void
   integrate(const Array<Real> & f, Array<Real> & intf,
             UInt nb_degree_of_freedom, ElementType type,
             GhostType ghost_type = _not_ghost,
             const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// integrate a scalar value on all elements of type "type"
   Real
   integrate(const Array<Real> & f, ElementType type,
             GhostType ghost_type = _not_ghost,
             const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// integrate one element scalar value on all elements of type "type"
   Real integrate(const Vector<Real> & f, ElementType type, UInt index,
                  GhostType ghost_type = _not_ghost) const override;
 
   /// integrate partially around an integration point (@f$ intf_q = f_q * J_q *
   /// w_q @f$)
   void integrateOnIntegrationPoints(
       const Array<Real> & f, Array<Real> & intf, UInt nb_degree_of_freedom,
       ElementType type, GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// interpolate on a phyiscal point inside an element
   void interpolate(const Vector<Real> & real_coords,
                    const Matrix<Real> & nodal_values,
                    Vector<Real> & interpolated,
                    const Element & element) const override;
 
   /// get the number of integration points
   UInt getNbIntegrationPoints(ElementType type,
                               GhostType ghost_type = _not_ghost) const override;
 
   /// get shapes precomputed
   const Array<Real> & getShapes(ElementType type,
                                 GhostType ghost_type = _not_ghost,
                                 UInt id = 0) const override;
 
   /// get the derivatives of shapes
   const Array<Real> & getShapesDerivatives(ElementType type,
                                            GhostType ghost_type = _not_ghost,
                                            UInt id = 0) const override;
 
   /// get integration points
   const inline Matrix<Real> &
   getIntegrationPoints(ElementType type,
                        GhostType ghost_type = _not_ghost) const override;
 
   /* ------------------------------------------------------------------------ */
   /* Shape method bridges                                                     */
   /* ------------------------------------------------------------------------ */
 
   /// compute the gradient of a nodal field on the integration points
   void gradientOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom,
       ElementType type, GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// interpolate a nodal field on the integration points
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom,
       ElementType type, GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// interpolate a nodal field on the integration points given a
   /// by_element_type
   void interpolateOnIntegrationPoints(
       const Array<Real> & u, ElementTypeMapArray<Real> & uq,
       const ElementTypeMapArray<UInt> * filter_elements =
           nullptr) const override;
 
   /// pre multiplies a tensor by the shapes derivaties
   void
-  computeBtD(const Array<Real> & Ds, Array<Real> & BtDs,
-             ElementType type, GhostType ghost_type,
+  computeBtD(const Array<Real> & Ds, Array<Real> & BtDs, ElementType type,
+             GhostType ghost_type,
              const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// left and right  multiplies a tensor by the shapes derivaties
   void computeBtDB(
       const Array<Real> & Ds, Array<Real> & BtDBs, UInt order_d,
       ElementType type, GhostType ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// left multiples a vector by the shape functions
   void
   computeNtb(const Array<Real> & bs, Array<Real> & Ntbs, ElementType type,
              GhostType ghost_type,
              const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// left and right  multiplies a tensor by the shapes
   void computeNtbN(
       const Array<Real> & bs, Array<Real> & NtbNs, ElementType type,
       GhostType ghost_type,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// compute the position of integration points given by an element_type_map
   /// from nodes position
   inline void computeIntegrationPointsCoordinates(
       ElementTypeMapArray<Real> & quadrature_points_coordinates,
       const ElementTypeMapArray<UInt> * filter_elements =
           nullptr) const override;
 
   /// compute the position of integration points from nodes position
   inline void computeIntegrationPointsCoordinates(
       Array<Real> & quadrature_points_coordinates, ElementType type,
       GhostType ghost_type = _not_ghost,
       const Array<UInt> & filter_elements = empty_filter) const override;
 
   /// interpolate field at given position (interpolation_points) from given
   /// values of this field at integration points (field)
   inline void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & result, GhostType ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const override;
 
   /// Interpolate field at given position from given values of this field at
   /// integration points (field)
   /// using matrices precomputed with
   /// initElementalFieldInterplationFromIntegrationPoints
   inline void interpolateElementalFieldFromIntegrationPoints(
       const ElementTypeMapArray<Real> & field,
       const ElementTypeMapArray<Real> &
           interpolation_points_coordinates_matrices,
       const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       ElementTypeMapArray<Real> & result, GhostType ghost_type,
       const ElementTypeMapArray<UInt> * element_filter) const override;
 
   /// Build pre-computed matrices for interpolation of field form integration
   /// points at other given positions (interpolation_points)
   inline void initElementalFieldInterpolationFromIntegrationPoints(
       const ElementTypeMapArray<Real> & interpolation_points_coordinates,
       ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
       ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
       const ElementTypeMapArray<UInt> * element_filter =
           nullptr) const override;
 
   /// find natural coords from real coords provided an element
   void inverseMap(const Vector<Real> & real_coords, UInt element,
                   ElementType type, Vector<Real> & natural_coords,
                   GhostType ghost_type = _not_ghost) const;
 
   /// return true if the coordinates provided are inside the element, false
   /// otherwise
   inline bool contains(const Vector<Real> & real_coords, UInt element,
                        ElementType type,
                        GhostType ghost_type = _not_ghost) const;
 
   /// compute the shape on a provided point
   inline void computeShapes(const Vector<Real> & real_coords, UInt element,
                             ElementType type, Vector<Real> & shapes,
                             GhostType ghost_type = _not_ghost) const override;
 
   /// compute the shape derivatives on a provided point
   inline void
   computeShapeDerivatives(const Vector<Real> & real_coords, UInt element,
                           ElementType type, Matrix<Real> & shape_derivatives,
                           GhostType ghost_type = _not_ghost) const override;
 
   /* ------------------------------------------------------------------------ */
   /* Other methods                                                            */
   /* ------------------------------------------------------------------------ */
   /// pre-compute normals on integration points
   void
   computeNormalsOnIntegrationPoints(GhostType ghost_type = _not_ghost) override;
   void
   computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                     GhostType ghost_type = _not_ghost) override;
   void computeNormalsOnIntegrationPoints(
       const Array<Real> & field, Array<Real> & normal, ElementType type,
       GhostType ghost_type = _not_ghost) const override;
   template <ElementType type>
   void computeNormalsOnIntegrationPoints(const Array<Real> & field,
                                          Array<Real> & normal,
                                          GhostType ghost_type) const;
 
 private:
   // To avoid a weird full specialization of a method in a non specalized class
   void computeNormalsOnIntegrationPointsPoint1(const Array<Real> & /*unused*/,
                                                Array<Real> & normal,
                                                GhostType ghost_type) const;
 
 public:
   /// function to print the contain of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   void assembleFieldLumped(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       ElementType type, GhostType ghost_type) const override;
 
   /// assemble a field as a matrix (ex. rho to mass matrix)
   void assembleFieldMatrix(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       ElementType type, GhostType ghost_type) const override;
 
   /// assemble a field as a lumped matrix (ex. rho in lumped mass)
   // template <class Functor>
   // void assembleFieldLumped(const Functor & field_funct, const ID & matrix_id,
   //                          const ID & dof_id, DOFManager & dof_manager,
   //                          ElementType type,
   //                          GhostType ghost_type) const;
 
   // /// assemble a field as a matrix (ex. rho to mass matrix)
   // template <class Functor>
   // void assembleFieldMatrix(const Functor & field_funct, const ID & matrix_id,
   //                          const ID & dof_id, DOFManager & dof_manager,
   //                          ElementType type,
   //                          GhostType ghost_type) const;
 
   // #ifdef AKANTU_STRUCTURAL_MECHANICS
   //   /// assemble a field as a matrix (ex. rho to mass matrix)
   //   void assembleFieldMatrix(const Array<Real> & field_1,
   //                            UInt nb_degree_of_freedom, SparseMatrix & M,
   //                            Array<Real> * n,
   //                            ElementTypeMapArray<Real> & rotation_mat,
   //                            ElementType type,
   //                            GhostType ghost_type = _not_ghost)
   //                            const;
 
   //   /// compute shapes function in a matrix for structural elements
   //   void
   //   computeShapesMatrix(ElementType type, UInt nb_degree_of_freedom,
   //                       UInt nb_nodes_per_element, Array<Real> * n, UInt id,
   //                       UInt degree_to_interpolate, UInt degree_interpolated,
   //                       const bool sign,
   //                       GhostType ghost_type = _not_ghost) const
   //                       override;
   // #endif
 
 private:
   friend struct fe_engine::details::AssembleLumpedTemplateHelper<kind>;
   friend struct fe_engine::details::AssembleFieldMatrixHelper<kind>;
   friend struct AssembleFieldMatrixStructHelper<kind, void>;
 
   /// templated function to compute the scaling to assemble a lumped matrix
   template <ElementType type>
   void assembleFieldLumped(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       GhostType ghost_type) const;
 
   /// @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j
   /// dV = \int \rho \varphi_i dV @f$
   template <ElementType type>
   void assembleLumpedRowSum(const Array<Real> & field, const ID & matrix_id,
                             const ID & dof_id, DOFManager & dof_manager,
                             GhostType ghost_type) const;
 
   /// @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
   template <ElementType type>
   void assembleLumpedDiagonalScaling(const Array<Real> & field,
                                      const ID & matrix_id, const ID & dof_id,
                                      DOFManager & dof_manager,
                                      GhostType ghost_type) const;
 
   /// assemble a field as a matrix (ex. rho to mass matrix)
   template <ElementType type>
   void assembleFieldMatrix(
       const std::function<void(Matrix<Real> &, const Element &)> & field_funct,
       const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
       GhostType ghost_type) const;
 
 #ifdef AKANTU_STRUCTURAL_MECHANICS
 
   /// assemble a field as a matrix for structural elements (ex. rho to mass
   /// matrix)
   template <ElementType type>
   void assembleFieldMatrix(const Array<Real> & field_1,
                            UInt nb_degree_of_freedom, SparseMatrix & M,
                            Array<Real> * n,
                            ElementTypeMapArray<Real> & rotation_mat,
                            __attribute__((unused)) GhostType ghost_type) const;
 
 #endif
 
   /* ------------------------------------------------------------------------ */
   /* Mesh Event Handler interface                                             */
   /* ------------------------------------------------------------------------ */
 public:
   void onElementsAdded(const Array<Element> & /*new_elements*/,
                        const NewElementsEvent & /*unused*/) override;
   void onElementsRemoved(const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const RemovedElementsEvent & /*unused*/) override;
   void onElementsChanged(const Array<Element> & /*unused*/,
                          const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const ChangedElementsEvent & /*unused*/) override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the shape class (probably useless: see getShapeFunction)
   const ShapeFunctions & getShapeFunctionsInterface() const override {
     return shape_functions;
   };
   /// get the shape class
   const Shape & getShapeFunctions() const { return shape_functions; };
 
   /// get the integrator class (probably useless: see getIntegrator)
   const Integrator & getIntegratorInterface() const override {
     return integrator;
   };
   /// get the integrator class
   const Integ & getIntegrator() const { return integrator; };
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   Integ integrator;
   Shape shape_functions;
 };
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "fe_engine_template_tmpl.hh"
 #include "fe_engine_template_tmpl_field.hh"
 /* -------------------------------------------------------------------------- */
 /* Shape Linked specialization                                                */
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 #include "fe_engine_template_tmpl_struct.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 /* Shape IGFEM specialization                                                 */
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_IGFEM)
 #include "fe_engine_template_tmpl_igfem.hh"
 #endif
 
 #endif /* AKANTU_FE_ENGINE_TEMPLATE_HH_ */
diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc
index 0abb6da9e..04e05dc5e 100644
--- a/src/model/structural_mechanics/structural_mechanics_model.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model.cc
@@ -1,627 +1,775 @@
 /**
  * @file   structural_mechanics_model.cc
  *
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Damien Spielmann <damien.spielmann@epfl.ch>
  *
  * @date creation: Fri Jul 15 2011
  * @date last modification: Mon Mar 15 2021
  *
  * @brief  Model implementation for Structural Mechanics elements
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "structural_mechanics_model.hh"
 #include "dof_manager.hh"
 #include "integrator_gauss.hh"
 #include "mesh.hh"
 #include "shape_structural.hh"
 #include "sparse_matrix.hh"
 #include "time_step_solver.hh"
 /* -------------------------------------------------------------------------- */
 #include "dumpable_inline_impl.hh"
 #include "dumper_elemental_field.hh"
 #include "dumper_internal_material_field.hh"
 #include "dumper_iohelper_paraview.hh"
 #include "group_manager_inline_impl.hh"
 /* -------------------------------------------------------------------------- */
 #include "structural_element_bernoulli_beam_2.hh"
 #include "structural_element_bernoulli_beam_3.hh"
 #include "structural_element_kirchhoff_shell.hh"
 /* -------------------------------------------------------------------------- */
 //#include "structural_mechanics_model_inline_impl.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 inline UInt StructuralMechanicsModel::getNbDegreeOfFreedom(ElementType type) {
   UInt ndof = 0;
 #define GET_(type) ndof = ElementClass<type>::getNbDegreeOfFreedom()
   AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_, _ek_structural);
 #undef GET_
 
   return ndof;
 }
 
 /* -------------------------------------------------------------------------- */
 StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim,
                                                    const ID & id)
     : Model(mesh, ModelType::_structural_mechanics_model, dim, id), f_m2a(1.0),
       stress("stress", id), element_material("element_material", id),
       set_ID("beam sets", id) {
   AKANTU_DEBUG_IN();
 
   registerFEEngineObject<MyFEEngineType>("StructuralMechanicsFEEngine", mesh,
                                          spatial_dimension);
 
   if (spatial_dimension == 2) {
     nb_degree_of_freedom = 3;
   } else if (spatial_dimension == 3) {
     nb_degree_of_freedom = 6;
   } else {
     AKANTU_TO_IMPLEMENT();
   }
 
   this->mesh.registerDumper<DumperParaview>("structural_mechanics_model", id,
                                             true);
   this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural);
 
   this->initDOFManager();
 
   this->dumper_default_element_kind = _ek_structural;
 
   mesh.getElementalData<Real>("extra_normal")
       .initialize(mesh, _element_kind = _ek_structural,
                   _nb_component = spatial_dimension, _with_nb_element = true,
                   _default_value = 0.);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 StructuralMechanicsModel::~StructuralMechanicsModel() = default;
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) {
   Model::initFullImpl(options);
 
   // Initializing stresses
   ElementTypeMap<UInt> stress_components;
 
   /// TODO this is ugly af, maybe add a function to FEEngine
   for (auto && type : mesh.elementTypes(_spatial_dimension = _all_dimensions,
                       _element_kind = _ek_structural)) {
     UInt nb_components = 0;
 
 // Getting number of components for each element type
 #define GET_(type) nb_components = ElementClass<type>::getNbStressComponents()
     AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_);
 #undef GET_
 
     stress_components(nb_components, type);
   }
 
   stress.initialize(
       getFEEngine(), _spatial_dimension = _all_dimensions,
       _element_kind = _ek_structural,
       _nb_component = [&stress_components](ElementType type,
                                            GhostType /*unused*/) -> UInt {
         return stress_components(type);
       });
 }
 
 /* -------------------------------------------------------------------------- */
 
 void StructuralMechanicsModel::initFEEngineBoundary() {
   /// TODO: this function should not be reimplemented
   /// we're just avoiding a call to Model::initFEEngineBoundary()
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::setTimeStep(Real time_step,
                                            const ID & solver_id) {
   Model::setTimeStep(time_step, solver_id);
   this->mesh.getDumper().setTimeStep(time_step);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Initialisation                                                             */
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::initSolver(
     TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) {
   AKANTU_DEBUG_IN();
 
   this->allocNodalField(displacement_rotation, nb_degree_of_freedom,
                         "displacement");
   this->allocNodalField(external_force, nb_degree_of_freedom, "external_force");
   this->allocNodalField(internal_force, nb_degree_of_freedom, "internal_force");
   this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs");
 
   auto & dof_manager = this->getDOFManager();
 
   if (not dof_manager.hasDOFs("displacement")) {
     dof_manager.registerDOFs("displacement", *displacement_rotation,
                              _dst_nodal);
     dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
   }
 
   if (time_step_solver_type == TimeStepSolverType::_dynamic ||
       time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
     this->allocNodalField(velocity, nb_degree_of_freedom, "velocity");
     this->allocNodalField(acceleration, nb_degree_of_freedom, "acceleration");
 
     if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
       dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
       dof_manager.registerDOFsDerivative("displacement", 2,
                                          *this->acceleration);
     }
+
+    /* Only allocate the mass if the "lumped" mode is ennabled.
+     *  Also it is not a 1D array, but has an element for every
+     *  DOF, which are most of the time equal, but makes handling
+     *  some operations a bit simpler. */
+    if (time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
+      this->allocateLumpedMassArray();
+    }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::initModel() {
   element_material.initialize(mesh, _element_kind = _ek_structural,
                               _default_value = 0, _with_nb_element = true);
 
   getFEEngine().initShapeFunctions(_not_ghost);
   getFEEngine().initShapeFunctions(_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleStiffnessMatrix() {
   AKANTU_DEBUG_IN();
 
   if (not need_to_reassemble_stiffness) {
     return;
   }
 
   if (not getDOFManager().hasMatrix("K")) {
     getDOFManager().getNewMatrix("K", getMatrixType("K"));
   }
 
   this->getDOFManager().zeroMatrix("K");
 
   for (const auto & type :
        mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) {
 #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix<type>();
 
     AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX);
 #undef ASSEMBLE_STIFFNESS_MATRIX
   }
 
   need_to_reassemble_stiffness = false;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::computeStresses() {
   AKANTU_DEBUG_IN();
 
   for (const auto & type :
        mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) {
 #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad<type>();
 
     AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD);
 #undef COMPUTE_STRESS_ON_QUAD
   }
 
   AKANTU_DEBUG_OUT();
 }
 
+/* -------------------------------------------------------------------------- */
+bool StructuralMechanicsModel::allocateLumpedMassArray() {
+  if (this->mass != nullptr) // Already allocated, so nothing to do.
+  {
+    return true;
+  };
+
+  // now allocate it
+  this->allocNodalField(this->mass, this->nb_degree_of_freedom, "lumped_mass");
+
+  return true;
+};
+
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field> StructuralMechanicsModel::createNodalFieldBool(
     const std::string & field_name, const std::string & group_name,
     __attribute__((unused)) bool padding_flag) {
 
   std::map<std::string, Array<bool> *> uint_nodal_fields;
   uint_nodal_fields["blocked_dofs"] = blocked_dofs.get();
 
   return mesh.createNodalField(uint_nodal_fields[field_name], group_name);
 }
 
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field>
 StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name,
                                                const std::string & group_name,
                                                bool padding_flag) {
 
   UInt n;
   if (spatial_dimension == 2) {
     n = 2;
   } else {
     n = 3;
   }
 
   UInt padding_size = 0;
   if (padding_flag) {
     padding_size = 3;
   }
 
   if (field_name == "displacement") {
     return mesh.createStridedNodalField(displacement_rotation.get(), group_name,
                                         n, 0, padding_size);
   }
 
   if (field_name == "velocity") {
     return mesh.createStridedNodalField(velocity.get(), group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "acceleration") {
     return mesh.createStridedNodalField(acceleration.get(), group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "rotation") {
     return mesh.createStridedNodalField(displacement_rotation.get(), group_name,
                                         nb_degree_of_freedom - n, n,
                                         padding_size);
   }
 
   if (field_name == "force") {
     return mesh.createStridedNodalField(external_force.get(), group_name, n, 0,
                                         padding_size);
   }
   if (field_name == "external_force") {
     return mesh.createStridedNodalField(external_force.get(), group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "momentum") {
     return mesh.createStridedNodalField(external_force.get(), group_name,
                                         nb_degree_of_freedom - n, n,
                                         padding_size);
   }
 
   if (field_name == "internal_force") {
     return mesh.createStridedNodalField(internal_force.get(), group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "internal_momentum") {
     return mesh.createStridedNodalField(internal_force.get(), group_name,
                                         nb_degree_of_freedom - n, n,
                                         padding_size);
   }
 
+  if (field_name == "mass") {
+    AKANTU_DEBUG_ASSERT(this->mass.get() != nullptr,
+                        "The lumped mass matrix was not allocated.");
+    return mesh.createStridedNodalField(this->mass.get(), group_name, n, 0,
+                                        padding_size);
+  }
+
   return nullptr;
 }
 
 /* -------------------------------------------------------------------------- */
 std::shared_ptr<dumpers::Field> StructuralMechanicsModel::createElementalField(
     const std::string & field_name, const std::string & group_name,
     bool /*unused*/, UInt spatial_dimension, ElementKind kind) {
 
   std::shared_ptr<dumpers::Field> field;
 
   if (field_name == "element_index_by_material") {
     field = mesh.createElementalField<UInt, Vector, dumpers::ElementalField>(
         field_name, group_name, spatial_dimension, kind);
   }
   if (field_name == "stress") {
     ElementTypeMap<UInt> nb_data_per_elem = this->mesh.getNbDataPerElem(stress);
 
     field = mesh.createElementalField<Real, dumpers::InternalMaterialField>(
         stress, group_name, this->spatial_dimension, kind, nb_data_per_elem);
   }
 
   return field;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Virtual methods from SolverCallback */
 /* -------------------------------------------------------------------------- */
 /// get the type of matrix needed
 MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) const {
   return _symmetric;
 }
 
 /// callback to assemble a Matrix
 void StructuralMechanicsModel::assembleMatrix(const ID & id) {
   if (id == "K") {
     assembleStiffnessMatrix();
   } else if (id == "M") {
     assembleMassMatrix();
   }
 }
 
 /// callback to assemble a lumped Matrix
-void StructuralMechanicsModel::assembleLumpedMatrix(const ID & /*id*/) {}
+void StructuralMechanicsModel::assembleLumpedMatrix(const ID & id) {
+  if ("M" == id) {
+    this->assembleLumpedMassMatrix();
+  }
+  return;
+}
 
 /// callback to assemble the residual StructuralMechanicsModel::(rhs)
 void StructuralMechanicsModel::assembleResidual() {
   AKANTU_DEBUG_IN();
 
   auto & dof_manager = getDOFManager();
 
   assembleInternalForce();
 
+  // Ensures that the matrix are assembled.
+  if (dof_manager.hasMatrix("K")) {
+    this->assembleMatrix("K");
+  }
+  if (dof_manager.hasMatrix("M")) {
+    this->assembleMatrix("M");
+  }
+  if (dof_manager.hasLumpedMatrix("M")) {
+    this->assembleLumpedMassMatrix();
+  }
+
+  /* This is essentially a summing up of forces
+   * first the external forces are counted for and then stored inside the
+   * residual.
+   */
   dof_manager.assembleToResidual("displacement", *external_force, 1);
   dof_manager.assembleToResidual("displacement", *internal_force, 1);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleResidual(const ID & residual_part) {
   AKANTU_DEBUG_IN();
 
   auto & dof_manager = this->getDOFManager();
 
   if ("external" == residual_part) {
-    dof_manager.assembleToResidual("displacement",
-                                             *this->external_force, 1);
+    dof_manager.assembleToResidual("displacement", *this->external_force, 1);
     AKANTU_DEBUG_OUT();
     return;
   }
 
   if ("internal" == residual_part) {
     this->assembleInternalForce();
-    dof_manager.assembleToResidual("displacement",
-                                             *this->internal_force, 1);
+    dof_manager.assembleToResidual("displacement", *this->internal_force, 1);
     AKANTU_DEBUG_OUT();
     return;
   }
 
   AKANTU_CUSTOM_EXCEPTION(
       debug::SolverCallbackResidualPartUnknown(residual_part));
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /* Virtual methods from Model                                                 */
 /* -------------------------------------------------------------------------- */
 /// get some default values for derived classes
 std::tuple<ID, TimeStepSolverType>
 StructuralMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
   switch (method) {
   case _static: {
     return std::make_tuple("static", TimeStepSolverType::_static);
   }
   case _implicit_dynamic: {
     return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
   }
+  case _explicit_lumped_mass: { // Taken from the solid mechanics part
+    return std::make_tuple("explicit_lumped",
+                           TimeStepSolverType::_dynamic_lumped);
+  }
+  case _explicit_consistent_mass: { // Taken from the solid mechanics part
+    return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
+  }
   default:
+    std::cout << "UNKOWN." << std::endl;
     return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
   }
 }
 
 /* ------------------------------------------------------------------------ */
 ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions(
     const TimeStepSolverType & type) const {
   ModelSolverOptions options;
 
   switch (type) {
+  case TimeStepSolverType::_dynamic_lumped: { // Taken from the solid mechanic
+                                              // part
+    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.non_linear_solver_type =
+        NonLinearSolverType::_newton_raphson; // _linear;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_pseudo_time;
     options.solution_type["displacement"] = IntegrationScheme::_not_defined;
     break;
   }
-  case TimeStepSolverType::_dynamic: {
+#if 1
+  case TimeStepSolverType::_dynamic: { // Copied from solid
+    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;
+  }
+#else
+  case TimeStepSolverType::_dynamic: { // Original
     options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_trapezoidal_rule_2;
     options.solution_type["displacement"] = IntegrationScheme::_displacement;
     break;
   }
+#endif
   default:
     AKANTU_EXCEPTION(type << " is not a valid time step solver type");
   }
 
   return options;
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleInternalForce() {
+
   internal_force->zero();
   computeStresses();
 
   for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions,
                    _element_kind = _ek_structural)) {
     assembleInternalForce(type, _not_ghost);
     // assembleInternalForce(type, _ghost);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleInternalForce(ElementType type,
                                                      GhostType gt) {
   auto & fem = getFEEngine();
   auto & sigma = stress(type, gt);
   auto ndof = getNbDegreeOfFreedom(type);
   auto nb_nodes = mesh.getNbNodesPerElement(type);
   auto ndof_per_elem = ndof * nb_nodes;
 
   Array<Real> BtSigma(fem.getNbIntegrationPoints(type) *
                           mesh.getNbElement(type),
                       ndof_per_elem, "BtSigma");
   fem.computeBtD(sigma, BtSigma, type, gt);
 
   Array<Real> intBtSigma(0, ndof_per_elem, "intBtSigma");
   fem.integrate(BtSigma, intBtSigma, ndof_per_elem, type, gt);
 
   getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force,
                                                    type, gt, -1.);
 }
 
 /* -------------------------------------------------------------------------- */
 Real StructuralMechanicsModel::getKineticEnergy() {
-  if (not this->getDOFManager().hasMatrix("M")) {
-    return 0.;
-  }
 
-  //if mass matrix was not assembled, assemble it now
+  const UInt nb_nodes = mesh.getNbNodes();
+  const UInt nb_degree_of_freedom = this->nb_degree_of_freedom;
+  Real ekin = 0.; // used to sum up energy (is divided by two at the very end)
+
+  // if mass matrix was not assembled, assemble it now
   this->assembleMassMatrix();
-  Real ekin = 0.;
-  UInt nb_nodes = mesh.getNbNodes();
 
-  Array<Real> Mv(nb_nodes, nb_degree_of_freedom);
-  this->getDOFManager().assembleMatMulVectToArray("displacement", "M",
-                                                  *this->velocity, Mv);
+  if (this->getDOFManager().hasLumpedMatrix("M")) {
+    /* This code computes the kinetic energy for the case when the mass is
+     * lumped. It is based on the solid mechanic equivalent.
+     */
+    AKANTU_DEBUG_ASSERT(this->mass != nullptr,
+                        "The lumped mass is not allocated.");
 
-  for (auto && data : zip(arange(nb_nodes), make_view(Mv, nb_degree_of_freedom),
-                          make_view(*this->velocity, nb_degree_of_freedom))) {
-    ekin += std::get<2>(data).dot(std::get<1>(data)) *
-            static_cast<Real>(mesh.isLocalOrMasterNode(std::get<0>(data)));
+    if (this->need_to_reassemble_lumpedMass) {
+      this->assembleLumpedMatrix("M");
+    }
+
+    /* Iterating over all nodes.
+     *   Important the velocity and mass also contains the rotational parts.
+     *   However, they can be handled in an uniform way. */
+    for (auto && data :
+         zip(arange(nb_nodes), make_view(*this->velocity, nb_degree_of_freedom),
+             make_view(*this->mass, nb_degree_of_freedom))) {
+      const UInt n = std::get<0>(data); // This is the ID of the current node
+
+      if (not mesh.isLocalOrMasterNode(
+              n)) // Only handle the node if it belongs to us.
+      {
+        continue;
+      }
+
+      const auto & v =
+          std::get<1>(data); // Get the velocity and mass of that node.
+      const auto & m = std::get<2>(data);
+      Real mv2 = 0.; // Contribution of this node.
+
+      for (UInt i = 0; i < nb_degree_of_freedom; ++i) {
+        /* In the solid mechanics part, only masses that are above a certain
+         * value are considered.
+         * However, the structural part, does not do this. */
+        const Real v_ = v(i);
+        const Real m_ = m(i);
+        mv2 += v_ * v_ * m_;
+      } // end for(i): going through the components
+
+      ekin += mv2; // add continution
+    }              // end for(n): iterating through all nodes
+  } else if (this->getDOFManager().hasMatrix("M")) {
+    /* Handle the case where no lumped mass is there.
+     * This is basically the original code.
+     */
+    if (this->need_to_reassemble_mass) {
+      this->assembleMassMatrix();
+    }
+
+    Array<Real> Mv(nb_nodes, nb_degree_of_freedom);
+    this->getDOFManager().assembleMatMulVectToArray("displacement", "M",
+                                                    *this->velocity, Mv);
+
+    for (auto && data :
+         zip(arange(nb_nodes), make_view(Mv, nb_degree_of_freedom),
+             make_view(*this->velocity, nb_degree_of_freedom))) {
+      if (mesh.isLocalOrMasterNode(std::get<0>(
+              data))) // only consider the node if we are blonging to it
+      {
+        ekin += std::get<2>(data).dot(std::get<1>(data));
+      }
+    }
+  } else {
+    /* This is the case where no mass is present, for whatever reason, such as
+     * the static case. We handle it specially be returning directly zero.
+     * However, by doing that there will not be a syncronizing event as in the
+     * other cases. Which is faster, but could be a problem in case the user
+     * expects this.
+     *
+     * Another not is, that the solid mechanics part, would generate an error in
+     * this clause. But, since the original implementation of the structural
+     * part, did not do that, I, Philip, decided to refrain from that. However,
+     * it is an option that should be considered.
+     */
+    return 0.;
   }
 
+  // Sum up across the comunicator
   mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum);
 
-  return ekin / 2.;
+  return ekin / 2.; // finally divide the energy by two
 }
 
 /* -------------------------------------------------------------------------- */
 Real StructuralMechanicsModel::getPotentialEnergy() {
   Real epot = 0.;
   UInt nb_nodes = mesh.getNbNodes();
 
-  //if stiffness matrix is not assembled, do it
+  // if stiffness matrix is not assembled, do it
   this->assembleStiffnessMatrix();
 
   Array<Real> Ku(nb_nodes, nb_degree_of_freedom);
   this->getDOFManager().assembleMatMulVectToArray(
       "displacement", "K", *this->displacement_rotation, Ku);
 
   for (auto && data :
        zip(arange(nb_nodes), make_view(Ku, nb_degree_of_freedom),
            make_view(*this->displacement_rotation, nb_degree_of_freedom))) {
     epot += std::get<2>(data).dot(std::get<1>(data)) *
             static_cast<Real>(mesh.isLocalOrMasterNode(std::get<0>(data)));
   }
 
   mesh.getCommunicator().allReduce(epot, SynchronizerOperation::_sum);
 
   return epot / 2.;
 }
 
 /* -------------------------------------------------------------------------- */
 Real StructuralMechanicsModel::getEnergy(const ID & energy) {
   if (energy == "kinetic") {
     return getKineticEnergy();
   }
 
   if (energy == "potential") {
     return getPotentialEnergy();
   }
 
   return 0;
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::computeForcesByLocalTractionArray(
     const Array<Real> & tractions, ElementType type) {
   AKANTU_DEBUG_IN();
 
   auto nb_element = getFEEngine().getMesh().getNbElement(type);
   auto nb_nodes_per_element =
       getFEEngine().getMesh().getNbNodesPerElement(type);
   auto nb_quad = getFEEngine().getNbIntegrationPoints(type);
 
   // check dimension match
   AKANTU_DEBUG_ASSERT(
       Mesh::getSpatialDimension(type) == getFEEngine().getElementDimension(),
       "element type dimension does not match the dimension of boundaries : "
           << getFEEngine().getElementDimension()
           << " != " << Mesh::getSpatialDimension(type));
 
   // check size of the vector
   AKANTU_DEBUG_ASSERT(
       tractions.size() == nb_quad * nb_element,
       "the size of the vector should be the total number of quadrature points");
 
   // check number of components
   AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == nb_degree_of_freedom,
                       "the number of components should be the spatial "
                       "dimension of the problem");
 
   Array<Real> Ntbs(nb_element * nb_quad,
                    nb_degree_of_freedom * nb_nodes_per_element);
 
   auto & fem = getFEEngine();
   fem.computeNtb(tractions, Ntbs, type);
 
   // allocate the vector that will contain the integrated values
   auto name = id + std::to_string(type) + ":integral_boundary";
   Array<Real> int_funct(nb_element, nb_degree_of_freedom * nb_nodes_per_element,
                         name);
 
   // do the integration
   getFEEngine().integrate(Ntbs, int_funct,
                           nb_degree_of_freedom * nb_nodes_per_element, type);
 
   // assemble the result into force vector
   getDOFManager().assembleElementalArrayLocalArray(int_funct, *external_force,
                                                    type, _not_ghost, 1);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::computeForcesByGlobalTractionArray(
     const Array<Real> & traction_global, ElementType type) {
   AKANTU_DEBUG_IN();
 
   UInt nb_element = mesh.getNbElement(type);
   UInt nb_quad = getFEEngine().getNbIntegrationPoints(type);
 
   Array<Real> traction_local(nb_element * nb_quad, nb_degree_of_freedom,
                              id + ":structuralmechanics:imposed_linear_load");
 
   auto R_it = getFEEngineClass<MyFEEngineType>()
                   .getShapeFunctions()
                   .getRotations(type)
                   .begin(nb_degree_of_freedom, nb_degree_of_freedom);
 
   auto Te_it = traction_global.begin(nb_degree_of_freedom);
   auto te_it = traction_local.begin(nb_degree_of_freedom);
 
   for (UInt e = 0; e < nb_element; ++e, ++R_it) {
     for (UInt q = 0; q < nb_quad; ++q, ++Te_it, ++te_it) {
       // turn the traction in the local referential
       te_it->template mul<false>(*R_it, *Te_it);
     }
   }
 
   computeForcesByLocalTractionArray(traction_local, type);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::afterSolveStep(bool converged) {
   if (converged) {
     assembleInternalForce();
   }
 }
 
 } // namespace akantu
diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh
index b3ac31d3f..3c80f2bd2 100644
--- a/src/model/structural_mechanics/structural_mechanics_model.hh
+++ b/src/model/structural_mechanics/structural_mechanics_model.hh
@@ -1,335 +1,398 @@
 /**
  * @file   structural_mechanics_model.hh
  *
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Philip Mueller <philip.paul.mueller@bluemail.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Damien Spielmann <damien.spielmann@epfl.ch>
  *
  * @date creation: Fri Jul 15 2011
  * @date last modification: Thu Apr 01 2021
  *
  * @brief  Particular implementation of the structural elements in the
  * StructuralMechanicsModel
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_named_argument.hh"
 #include "boundary_condition.hh"
 #include "model.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_
 #define AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 class Material;
 class MaterialSelector;
 class DumperIOHelper;
 class NonLocalManager;
 template <ElementKind kind, class IntegrationOrderFunctor>
 class IntegratorGauss;
 template <ElementKind kind> class ShapeStructural;
 } // namespace akantu
 
 namespace akantu {
 
 struct StructuralMaterial {
   Real E{0};
   Real A{1};
   Real I{0};
   Real Iz{0};
   Real Iy{0};
   Real GJ{0};
   Real rho{0};
   Real t{0};
   Real nu{0};
 };
 
 class StructuralMechanicsModel : public Model {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   using MyFEEngineType =
       FEEngineTemplate<IntegratorGauss, ShapeStructural, _ek_structural>;
 
   StructuralMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions,
                            const ID & id = "structural_mechanics_model");
 
   ~StructuralMechanicsModel() override;
 
   /// Init full model
   void initFullImpl(const ModelOptions & options) override;
 
   /// Init boundary FEEngine
   void initFEEngineBoundary() override;
 
   /* ------------------------------------------------------------------------ */
   /* Virtual methods from SolverCallback                                      */
   /* ------------------------------------------------------------------------ */
   /// get the type of matrix needed
   MatrixType getMatrixType(const ID & matrix_id) const override;
 
   /// callback to assemble a Matrix
   void assembleMatrix(const ID & matrix_id) override;
 
   /// callback to assemble a lumped Matrix
   void assembleLumpedMatrix(const ID & matrix_id) override;
 
   /// callback to assemble the residual (rhs)
   void assembleResidual() override;
 
   void assembleResidual(const ID & residual_part) override;
 
   bool canSplitResidual() const override { return true; }
 
   void afterSolveStep(bool converged) override;
 
   /// compute kinetic energy
   Real getKineticEnergy();
 
   /// compute potential energy
   Real getPotentialEnergy();
 
   /// compute the specified energy
   Real getEnergy(const ID & energy);
 
+  /**
+   * \brief This function computes the an approximation of the lumped mass.
+   *
+   * The mass is computed by looping over all beams and computing their mass.
+   * The mass of a single beam is computed by the (initial) length of the beam,
+   * its cross sectional area and its density.
+   * The beam mass is then equaly distributed among the two nodes.
+   *
+   * For computing the rotational inertia, the function assumes that the mass of
+   * a node is uniformaly distributed inside a disc (2D) or a sphere (3D). The
+   * size of that disc, depends on the volume of the beam.
+   *
+   * Note that the computation of the mass is not unambigius.
+   * The reason for this is, that the units of `StructralMaterial::rho` are not
+   * clear. By default the function assumes that its unit are 'Mass per Volume'.
+   * However, this makes the computed mass different than the consistent mass,
+   * which seams to assume that its units are 'mass per unit length'.
+   * The main difference between thge two are not the values, but that the
+   * first version depends on `StructuralMaterial::A` while the later does not.
+   * By defining the macro `AKANTU_STRUCTURAL_MECHANICS_CONSISTENT_LUMPED_MASS`
+   * the function will compute the mass in a way that is consistent with the
+   * consistent mass matrix.
+   *
+   * \note	The lumped mass is not stored inside the DOFManager.
+   *
+   * \param  ghost_type 	Should ghost types be computed.
+   */
+  void assembleLumpedMassMatrix();
+
   /* ------------------------------------------------------------------------ */
   /* Virtual methods from Model                                               */
   /* ------------------------------------------------------------------------ */
 protected:
   /// get some default values for derived classes
   std::tuple<ID, TimeStepSolverType>
   getDefaultSolverID(const AnalysisMethod & method) override;
 
   ModelSolverOptions
   getDefaultSolverOptions(const TimeStepSolverType & type) const override;
 
   static UInt getNbDegreeOfFreedom(ElementType type);
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
   void initSolver(TimeStepSolverType time_step_solver_type,
                   NonLinearSolverType non_linear_solver_type) override;
 
   /// initialize the model
   void initModel() override;
 
   /// compute the stresses per elements
   void computeStresses();
 
   /// compute the nodal forces
   void assembleInternalForce();
 
   /// compute the nodal forces for an element type
   void assembleInternalForce(ElementType type, GhostType gt);
 
   /// assemble the stiffness matrix
   void assembleStiffnessMatrix();
 
   /// assemble the mass matrix for consistent mass resolutions
   void assembleMassMatrix();
 
 protected:
   /// assemble the mass matrix for either _ghost or _not_ghost elements
   void assembleMassMatrix(GhostType ghost_type);
 
   /// computes rho
   void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type);
 
   /// finish the computation of residual to solve in increment
   void updateResidualInternal();
 
   /* ------------------------------------------------------------------------ */
 private:
   template <ElementType type> void assembleStiffnessMatrix();
   template <ElementType type> void computeStressOnQuad();
   template <ElementType type>
   void computeTangentModuli(Array<Real> & tangent_moduli);
 
   /* ------------------------------------------------------------------------ */
   /* Dumpable interface                                                       */
   /* ------------------------------------------------------------------------ */
 public:
   std::shared_ptr<dumpers::Field>
   createNodalFieldReal(const std::string & field_name,
                        const std::string & group_name,
                        bool padding_flag) override;
 
   std::shared_ptr<dumpers::Field>
   createNodalFieldBool(const std::string & field_name,
                        const std::string & group_name,
                        bool padding_flag) override;
 
   std::shared_ptr<dumpers::Field>
   createElementalField(const std::string & field_name,
                        const std::string & group_name, bool padding_flag,
                        UInt spatial_dimension, ElementKind kind) override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// set the value of the time step
   void setTimeStep(Real time_step, const ID & solver_id = "") override;
 
   /// get the StructuralMechanicsModel::displacement vector
   AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array<Real> &);
 
   /// get the StructuralMechanicsModel::velocity vector
   AKANTU_GET_MACRO(Velocity, *velocity, Array<Real> &);
 
   /// get    the    StructuralMechanicsModel::acceleration    vector,   updated
   /// by
   /// StructuralMechanicsModel::updateAcceleration
   AKANTU_GET_MACRO(Acceleration, *acceleration, Array<Real> &);
 
   /// get the StructuralMechanicsModel::external_force vector
   AKANTU_GET_MACRO(ExternalForce, *external_force, Array<Real> &);
 
   /// get the StructuralMechanicsModel::internal_force vector (boundary forces)
   AKANTU_GET_MACRO(InternalForce, *internal_force, Array<Real> &);
 
   /// get the StructuralMechanicsModel::boundary vector
   AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<bool> &);
 
+  /**
+   * Returns a const reference to the array that stores the lumped mass.
+   *
+   * The returned array has dimension `N x d` where `N` is the number of nodes
+   * and `d`, is the number of degrees of freedom per node.
+   */
+  inline const Array<Real> & getLumpedMass() const {
+    if (this->mass == nullptr) {
+      AKANTU_EXCEPTION("The pointer to the mass was not allocated.");
+    };
+    return *(this->mass);
+  };
+
+  // These function is an alias, for compability with the solid mechanics
+  inline const Array<Real> & getMass() const { return this->getLumpedMass(); }
+
+  // Creates the array for storing the mass
+  bool allocateLumpedMassArray();
+
+  /**
+   * Tests if *this has a lumped mass pointer.
+   */
+  inline bool hasLumpedMass() const { return (this->mass != nullptr); };
+
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(RotationMatrix, rotation_matrix, Real);
 
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
 
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt);
 
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Set_ID, set_ID, UInt);
 
   /**
    * \brief This function adds the `StructuralMaterial` material to the list of
    * materials managed by *this.
    *
    * It is important that this function might invalidate all references to
    * structural materials, that were previously optained by `getMaterial()`.
    *
    * \param  material The new material.
    *
    * \return The ID of the material that was added.
    *
    * \note The return type is is new.
    */
   UInt addMaterial(StructuralMaterial & material, const ID & name = "");
 
   const StructuralMaterial &
   getMaterialByElement(const Element & element) const;
 
   /**
    * \brief Returns the ith material of *this.
    * \param i The ith material
    */
   const StructuralMaterial & getMaterial(UInt material_index) const;
 
   const StructuralMaterial & getMaterial(const ID & name) const;
 
   /**
    * \brief Returns the number of the different materials inside *this.
    */
   UInt getNbMaterials() const { return materials.size(); }
 
   /* ------------------------------------------------------------------------ */
   /* Boundaries (structural_mechanics_model_boundary.cc)                      */
   /* ------------------------------------------------------------------------ */
 public:
   /// Compute Linear load function set in global axis
   void computeForcesByGlobalTractionArray(const Array<Real> & traction_global,
                                           ElementType type);
 
   /// Compute Linear load function set in local axis
   void computeForcesByLocalTractionArray(const Array<Real> & tractions,
                                          ElementType type);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   /// time step
   Real time_step;
 
   /// conversion coefficient form force/mass to acceleration
   Real f_m2a;
 
   /// displacements array
   std::unique_ptr<Array<Real>> displacement_rotation;
 
   /// velocities array
   std::unique_ptr<Array<Real>> velocity;
 
   /// accelerations array
   std::unique_ptr<Array<Real>> acceleration;
 
   /// forces array
   std::unique_ptr<Array<Real>> internal_force;
 
   /// forces array
   std::unique_ptr<Array<Real>> external_force;
 
-  /// lumped mass array
+  /**
+   * \brief	This is the "lumped" mass array.
+   *
+   * It is a bit special, since it is not a one dimensional array, bit it is
+   * actually a matrix. The number of rows equals the number of nodes. The
+   * number of colums equals the number of degrees of freedoms per node. This
+   * layout makes the thing a bit more simple.
+   *
+   * Note that it is only allocated in case, the "Lumped" mode is enabled.
+   */
   std::unique_ptr<Array<Real>> mass;
 
   /// boundaries array
   std::unique_ptr<Array<bool>> blocked_dofs;
 
   /// stress array
   ElementTypeMapArray<Real> stress;
 
   ElementTypeMapArray<UInt> element_material;
 
   // Define sets of beams
   ElementTypeMapArray<UInt> set_ID;
 
   /// number of degre of freedom
   UInt nb_degree_of_freedom;
 
   // Rotation matrix
   ElementTypeMapArray<Real> rotation_matrix;
 
   // /// analysis method check the list in akantu::AnalysisMethod
   // AnalysisMethod method;
 
   /// flag defining if the increment must be computed or not
   bool increment_flag;
 
   bool need_to_reassemble_mass{true};
   bool need_to_reassemble_stiffness{true};
+  bool need_to_reassemble_lumpedMass{true};
 
   /* ------------------------------------------------------------------------ */
   std::vector<StructuralMaterial> materials;
   std::map<std::string, UInt> materials_names_to_id;
 };
 
 } // namespace akantu
 
 #include "structural_mechanics_model_inline_impl.hh"
 
 #endif /* AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ */
diff --git a/src/model/structural_mechanics/structural_mechanics_model_mass.cc b/src/model/structural_mechanics/structural_mechanics_model_mass.cc
index d444e27e5..9aad97666 100644
--- a/src/model/structural_mechanics/structural_mechanics_model_mass.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model_mass.cc
@@ -1,91 +1,219 @@
 /**
  * @file   structural_mechanics_model_mass.cc
  *
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Mon Jul 07 2014
  * @date last modification: Thu Mar 04 2021
  *
  * @brief  function handling mass computation
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2014-2021 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_math.hh"
 #include "integrator_gauss.hh"
 #include "material.hh"
 #include "shape_structural.hh"
 #include "structural_mechanics_model.hh"
-
 /* -------------------------------------------------------------------------- */
+
 namespace akantu {
 
 class ComputeRhoFunctorStruct {
 public:
   explicit ComputeRhoFunctorStruct(const StructuralMechanicsModel & model)
       : model(model){};
 
   void operator()(Matrix<Real> & rho, const Element & element) const {
-    Real mat_rho = model.getMaterialByElement(element).rho;
-    rho.set(mat_rho);
+    const auto & mat = model.getMaterialByElement(element);
+    const Real mat_rho = mat.rho; // This is the density
+    const Real mat_A = mat.A;     // This is the cross section
+    const Real mat_mass_per_length =
+        mat_rho * mat_A; // Mass of the beam per unit length
+
+    rho.set(
+        mat_mass_per_length); // The integrator _recquiers_ mass per unit length
   }
 
 private:
   const StructuralMechanicsModel & model;
 };
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleMassMatrix() {
   AKANTU_DEBUG_IN();
 
   if (not need_to_reassemble_mass) {
     return;
   }
 
   if (not getDOFManager().hasMatrix("M")) {
     getDOFManager().getNewMatrix("M", getMatrixType("M"));
   }
 
   this->getDOFManager().zeroMatrix("M");
   assembleMassMatrix(_not_ghost);
 
   need_to_reassemble_mass = false;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void StructuralMechanicsModel::assembleMassMatrix(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
   auto & fem = getFEEngineClass<MyFEEngineType>();
   ComputeRhoFunctorStruct compute_rho(*this);
 
   for (auto type :
        mesh.elementTypes(spatial_dimension, ghost_type, _ek_structural)) {
     fem.assembleFieldMatrix(compute_rho, "M", "displacement",
                             this->getDOFManager(), type, ghost_type);
   }
   AKANTU_DEBUG_OUT();
 }
 
+/* -------------------------------------------------------------------------- */
+void StructuralMechanicsModel::assembleLumpedMassMatrix() {
+
+  if (not this->need_to_reassemble_lumpedMass)
+    return;
+
+  allocNodalField(this->mass, nb_degree_of_freedom, "lumped_mass");
+
+  if (!this->getDOFManager().hasLumpedMatrix("M")) {
+    this->getDOFManager().getNewLumpedMatrix("M");
+  }
+
+  this->getDOFManager().zeroLumpedMatrix("M");
+
+  // Preparing
+  const UInt nb_nodes = mesh.getNbNodes();
+  const auto & nodes = mesh.getNodes();
+  const auto & node_it = make_view(nodes, spatial_dimension).begin();
+  auto & lumped_mass = *(this->mass);
+
+  lumped_mass.zero(); // Set the lumped mass to zero
+
+  Array<Real> volumes(nb_nodes, 1, 0., "volumes");
+
+  /// We now compute the mass and the volume, but not the inertia
+  for (auto type :
+       mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) {
+    const auto & element_material_id = this->element_material(type);
+    const auto & connectivity = mesh.getConnectivity(type);
+    const UInt nb_element = connectivity.size();
+
+    if (type != _bernoulli_beam_3 or type != _bernoulli_beam_2) {
+      AKANTU_EXCEPTION(
+          "The lumped mass was not implemented for non Bernoulli Beams");
+    }
+
+    // Now iterate through all the connectivity
+    for (auto && data : zip(make_view(connectivity, 2), element_material_id)) {
+      const auto & conn = std::get<0>(data);
+      auto node1 = conn(0);
+      auto node2 = conn(0);
+
+      const auto & material = this->materials.at(std::get<1>(data));
+      const Real rho = material.rho;
+      const Real cross_section = material.A;
+
+      const Vector<Real> coord_node_1 = node_it[node1];
+      const Vector<Real> coord_node_2 = node_it[node2];
+
+      const Real length = coord_node_1.distance(coord_node_2);
+      const Real volume = length * cross_section;
+      const Real mass = volume * rho;
+
+      // Now distribute the mass at the right places
+      for (UInt d = 0; d != spatial_dimension; ++d) {
+        lumped_mass(node1, d) += mass / 2.;
+        lumped_mass(node2, d) += mass / 2.;
+      }; // end for(d):
+
+      volumes(node1) += volume / 2.; // this entry never point to mass
+      volumes(node2) += volume / 2.;
+    }
+  }
+
+  const Real pi = std::atan(1.) * 4;
+
+  // We now compute the inertia.
+  if (spatial_dimension == 2) {
+    /* This is the 2D case, so we are assuming that the mass is inside a disc.
+     * Which is given as
+     * \begin{align}
+     *   I := m \cdot r^2
+     * \end{align}
+     *
+     * The radius is obtained by assuming that the volume, that is associated to
+     * a beam, forms a uniformly disc. From this volume we can compute the
+     * radius.
+     */
+    for (auto && data :
+         zip(make_view(lumped_mass, nb_degree_of_freedom), volumes)) {
+      const Real volume = std::get<1>(data);
+      auto & masses = std::get<0>(data);
+
+      const Real r2 = volume / pi; // The square of the volume
+
+      const Real inertia = masses(0) * r2;
+      masses(spatial_dimension) = inertia;
+    }
+  } else if (spatial_dimension == 3) {
+    /* This is essentially the same as for 2D only that we assume here,
+     * that the mass is uniformly distributed inside a sphere.
+     * And thus we have
+     * \begin{align}
+     *   I := \frac{2}{5} m \cdot r^2
+     * \end{align}
+     *
+     * We also have to set three values, for the three axis.
+     * But since we assume a sphere, they are all the same.
+     */
+    for (auto && data :
+         zip(make_view(lumped_mass, nb_degree_of_freedom), volumes)) {
+      const Real volume = std::get<1>(data);
+      auto & masses = std::get<0>(data);
+
+      const Real r2 = std::pow((volume * 3.) / (4 * pi), 2. / 3.);
+      const Real inertia = (2. / 5.) * masses(0) * r2;
+
+      for (UInt d = spatial_dimension; d < masses.size(); ++d) {
+        masses(d) = inertia;
+      }
+    }
+  } else {
+    AKANTU_EXCEPTION("The dimension " << spatial_dimension << " is not known.");
+  }
+
+  this->getDOFManager().assembleToLumpedMatrix("displacement", lumped_mass,
+                                               "M");
+  this->need_to_reassemble_lumpedMass = false;
+  return;
+}
+
 } // namespace akantu
diff --git a/src/solver/solver_petsc.cc b/src/solver/solver_petsc.cc
index 696c75c24..baa032bde 100644
--- a/src/solver/solver_petsc.cc
+++ b/src/solver/solver_petsc.cc
@@ -1,93 +1,100 @@
 /**
  * @file   solver_petsc.cc
  *
  * @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue May 13 2014
  * @date last modification: Tue May 26 2020
  *
  * @brief  Solver class implementation for the petsc solver
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2014-2021 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 "solver_petsc.hh"
 #include "dof_manager_petsc.hh"
 #include "mpi_communicator_data.hh"
 #include "solver_vector_petsc.hh"
 #include "sparse_matrix_petsc.hh"
 /* -------------------------------------------------------------------------- */
 #include <petscksp.h>
 //#include <petscsys.h>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 SolverPETSc::SolverPETSc(DOFManagerPETSc & dof_manager, const ID & matrix_id,
                          const ID & id)
     : SparseSolver(dof_manager, matrix_id, id), dof_manager(dof_manager),
       matrix(dof_manager.getMatrix(matrix_id)) {
   auto && mpi_comm = dof_manager.getMPIComm();
 
   /// create a solver context
   PETSc_call(KSPCreate, mpi_comm, &this->ksp);
 }
 
 /* -------------------------------------------------------------------------- */
 SolverPETSc::~SolverPETSc() {
   if (ksp != nullptr) {
     PETSc_call(KSPDestroy, &ksp);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolverPETSc::setOperators() {
   // set the matrix that defines the linear system and the matrix for
 // preconditioning (here they are the same)
 #if PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 5
   PETSc_call(KSPSetOperators, ksp, this->matrix.getMat(),
              this->matrix.getMat());
 #else
   PETSc_call(KSPSetOperators, ksp, this->matrix.getMat(), this->matrix.getMat(),
              SAME_NONZERO_PATTERN);
 #endif
 
   // If this is not called the solution vector is zeroed in the call to
   // KSPSolve().
   PETSc_call(KSPSetInitialGuessNonzero, ksp, PETSC_TRUE);
   PETSc_call(KSPSetFromOptions, ksp);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolverPETSc::solve() {
   Vec & rhs(this->dof_manager.getResidual());
   Vec & solution(this->dof_manager.getSolution());
 
   PETSc_call(KSPSolve, ksp, rhs, solution);
 }
 
+/* -------------------------------------------------------------------------- */
+bool SolverPETSc::isFinite() const {
+  PetscReal norm;
+  PETSc_call(VecNorm, NORM_INFINITY, &norm);
+  return std::isfinite(norm);
+}
+
 } // namespace akantu
diff --git a/src/solver/solver_petsc.hh b/src/solver/solver_petsc.hh
index e0b573f35..d14feca41 100644
--- a/src/solver/solver_petsc.hh
+++ b/src/solver/solver_petsc.hh
@@ -1,82 +1,84 @@
 /**
  * @file   solver_petsc.hh
  *
  * @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue May 13 2014
  * @date last modification: Tue May 26 2020
  *
  * @brief  Solver class interface for the petsc solver
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2014-2021 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 "sparse_solver.hh"
 /* -------------------------------------------------------------------------- */
 #include <petscksp.h>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SOLVER_PETSC_HH_
 #define AKANTU_SOLVER_PETSC_HH_
 
 namespace akantu {
 class SparseMatrixPETSc;
 class DOFManagerPETSc;
 } // namespace akantu
 
 namespace akantu {
 
 class SolverPETSc : public SparseSolver {
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   SolverPETSc(DOFManagerPETSc & dof_manager, const ID & matrix_id,
               const ID & id = "solver_petsc");
 
   ~SolverPETSc() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// create the solver context and set the matrices
   virtual void setOperators();
   void solve() override;
 
+  bool isFinite() const override;
+
 private:
   /// DOFManager correctly typed
   DOFManagerPETSc & dof_manager;
 
   /// PETSc linear solver
   KSP ksp;
 
   /// Matrix defining the system of equations
   SparseMatrixPETSc & matrix;
 };
 
 } // namespace akantu
 
 #endif /* AKANTU_SOLVER_PETSC_HH_ */
diff --git a/src/solver/solver_vector.hh b/src/solver/solver_vector.hh
index fd43c5ac5..d05727fa3 100644
--- a/src/solver/solver_vector.hh
+++ b/src/solver/solver_vector.hh
@@ -1,105 +1,96 @@
 /**
  * @file   solver_vector.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Feb 21 2013
  * @date last modification: Tue May 26 2020
  *
  * @brief  Solver vector interface base class
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SOLVER_VECTOR_HH_
 #define AKANTU_SOLVER_VECTOR_HH_
 
 namespace akantu {
 class DOFManager;
 }
 
 namespace akantu {
 
 class SolverVector {
 public:
   SolverVector(DOFManager & dof_manager, const ID & id = "solver_vector")
       : id(id), _dof_manager(dof_manager) {}
 
   SolverVector(const SolverVector & vector, const ID & id = "solver_vector")
       : id(id), _dof_manager(vector._dof_manager) {}
 
   virtual ~SolverVector() = default;
 
   // resize the vector to the size of the problem
   virtual void resize() = 0;
 
   // clear the vector
   virtual void set(Real val) = 0;
   void zero() { this->set({}); }
 
   virtual operator const Array<Real> &() const = 0;
 
   virtual Int size() const = 0;
   virtual Int localSize() const = 0;
 
   virtual SolverVector & operator+(const SolverVector & y) = 0;
   virtual SolverVector & operator=(const SolverVector & y) = 0;
 
   UInt & release() { return release_; }
   UInt release() const { return release_; }
 
   virtual void printself(std::ostream & stream, int indent = 0) const = 0;
 
+  virtual bool isFinite() const = 0;
 
-  /**
-   * \brief	Returns `true` if `*this` is distributed or not.
-   *
-   * The default implementation returns `false`.
-   *
-   */
-  virtual
-  bool
-  isDistributed()
-    const
-      { return false; }
-
+  /// Returns `true` if `*this` is distributed or not.
+  virtual bool isDistributed() const { return false; }
 
 protected:
   ID id;
 
   /// Underlying dof manager
   DOFManager & _dof_manager;
 
   UInt release_{0};
 };
 
 inline std::ostream & operator<<(std::ostream & stream, SolverVector & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #endif /* AKANTU_SOLVER_VECTOR_HH_ */
diff --git a/src/solver/solver_vector_default.hh b/src/solver/solver_vector_default.hh
index 932b87288..49e11dd24 100644
--- a/src/solver/solver_vector_default.hh
+++ b/src/solver/solver_vector_default.hh
@@ -1,154 +1,156 @@
 /**
  * @file   solver_vector_default.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Jan 01 2019
  * @date last modification: Sat May 23 2020
  *
  * @brief  Solver vector interface to Array
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-2021 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 "solver_vector.hh"
 /* -------------------------------------------------------------------------- */
 #include <utility>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SOLVER_VECTOR_DEFAULT_HH_
 #define AKANTU_SOLVER_VECTOR_DEFAULT_HH_
 
 namespace akantu {
 class DOFManagerDefault;
 } // namespace akantu
 
 namespace akantu {
 
 class SolverVectorArray : public SolverVector {
 public:
   SolverVectorArray(DOFManagerDefault & dof_manager, const ID & id);
   SolverVectorArray(const SolverVectorArray & vector, const ID & id);
 
   ~SolverVectorArray() override = default;
 
   virtual Array<Real> & getVector() = 0;
   virtual const Array<Real> & getVector() const = 0;
 
   void printself(std::ostream & stream, int indent = 0) const override {
     std::string space(indent, AKANTU_INDENT);
     stream << space << "SolverVectorArray [" << std::endl;
     stream << space << " + id: " << id << std::endl;
     this->getVector().printself(stream, indent + 1);
     stream << space << "]" << std::endl;
   }
 
   using SolverVector::isDistributed;
 };
 
 /* -------------------------------------------------------------------------- */
 template <class Array_> class SolverVectorArrayTmpl : public SolverVectorArray {
 public:
   SolverVectorArrayTmpl(DOFManagerDefault & dof_manager, Array_ & vector,
                         const ID & id = "solver_vector_default")
       : SolverVectorArray(dof_manager, id), dof_manager(dof_manager),
         vector(vector) {}
 
   template <class A = Array_,
             std::enable_if_t<not std::is_reference<A>::value> * = nullptr>
   SolverVectorArrayTmpl(DOFManagerDefault & dof_manager,
                         const ID & id = "solver_vector_default")
       : SolverVectorArray(dof_manager, id), dof_manager(dof_manager),
         vector(0, 1, id + ":vector") {}
 
   SolverVectorArrayTmpl(const SolverVectorArrayTmpl & vector,
                         const ID & id = "solver_vector_default")
       : SolverVectorArray(vector, id), dof_manager(vector.dof_manager),
         vector(vector.vector) {}
 
   operator const Array<Real> &() const override { return getVector(); };
   virtual operator Array<Real> &() { return getVector(); };
 
   SolverVector & operator+(const SolverVector & y) override;
   SolverVector & operator=(const SolverVector & y) override;
 
   void resize() override {
     static_assert(not std::is_const<std::remove_reference_t<Array_>>::value,
                   "Cannot resize a const Array");
     this->vector.resize(this->localSize(), 0.);
     ++this->release_;
   }
 
   void set(Real val) override {
     static_assert(not std::is_const<std::remove_reference_t<Array_>>::value,
                   "Cannot clear a const Array");
     this->vector.set(val);
     ++this->release_;
   }
 
   virtual
   bool
   isDistributed()
     const
     override
       { return false; }
 
 
 public:
   Array<Real> & getVector() override { return vector; }
   const Array<Real> & getVector() const override { return vector; }
 
   Int size() const override;
   Int localSize() const override;
 
   virtual Array<Real> & getGlobalVector() { return this->vector; }
   virtual void setGlobalVector(const Array<Real> & solution) {
     this->vector.copy(solution);
   }
 
+  bool isFinite() const override { return vector.isFinite(); }
+
 protected:
   DOFManagerDefault & dof_manager;
   Array_ vector;
 
   template <class A> friend class SolverVectorArrayTmpl;
 };
 
 /* -------------------------------------------------------------------------- */
 using SolverVectorDefault = SolverVectorArrayTmpl<Array<Real>>;
 
 /* -------------------------------------------------------------------------- */
 template <class Array>
 using SolverVectorDefaultWrap = SolverVectorArrayTmpl<Array &>;
 
 template <class Array>
 decltype(auto) make_solver_vector_default_wrap(DOFManagerDefault & dof_manager,
                                                Array & vector) {
   return SolverVectorDefaultWrap<Array>(dof_manager, vector);
 }
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 #include "solver_vector_default_tmpl.hh"
 /* -------------------------------------------------------------------------- */
 
 #endif /* AKANTU_SOLVER_VECTOR_DEFAULT_HH_ */
diff --git a/src/solver/solver_vector_distributed.cc b/src/solver/solver_vector_distributed.cc
index 4eacf65a1..7d2865eef 100644
--- a/src/solver/solver_vector_distributed.cc
+++ b/src/solver/solver_vector_distributed.cc
@@ -1,81 +1,90 @@
 /**
  * @file   solver_vector_distributed.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Jan 01 2019
  * @date last modification: Sat May 23 2020
  *
  * @brief  Solver vector interface for distributed arrays
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-2021 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 "solver_vector_distributed.hh"
 #include "dof_manager_default.hh"
 #include "dof_synchronizer.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 SolverVectorDistributed::SolverVectorDistributed(
     DOFManagerDefault & dof_manager, const ID & id)
     : SolverVectorDefault(dof_manager, id) {}
 
 /* -------------------------------------------------------------------------- */
 SolverVectorDistributed::SolverVectorDistributed(
     const SolverVectorDefault & vector, const ID & id)
     : SolverVectorDefault(vector, id) {}
 
 /* -------------------------------------------------------------------------- */
 Array<Real> & SolverVectorDistributed::getGlobalVector() {
   auto & synchronizer = dof_manager.getSynchronizer();
 
   if (not this->global_vector) {
     this->global_vector =
         std::make_unique<Array<Real>>(0, 1, "global_residual");
   }
 
   if (synchronizer.getCommunicator().whoAmI() == 0) {
     this->global_vector->resize(dof_manager.getSystemSize());
     synchronizer.gather(this->vector, *this->global_vector);
   } else {
     synchronizer.gather(this->vector);
   }
 
   return *this->global_vector;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolverVectorDistributed::setGlobalVector(const Array<Real> & solution) {
   auto & synchronizer = dof_manager.getSynchronizer();
   if (synchronizer.getCommunicator().whoAmI() == 0) {
     synchronizer.scatter(this->vector, solution);
   } else {
     synchronizer.scatter(this->vector);
   }
 }
 
+/* -------------------------------------------------------------------------- */
+bool SolverVectorDistributed::isFinite() const {
+  auto & synchronizer = dof_manager.getSynchronizer();
+  bool is_finite = this->vector.isFinite();
+  synchronizer.getCommunicator().allReduce(is_finite,
+                                           SynchronizerOperation::_land);
+  return is_finite;
+}
+
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu
diff --git a/src/solver/solver_vector_distributed.hh b/src/solver/solver_vector_distributed.hh
index e7b988144..a928e0ac8 100644
--- a/src/solver/solver_vector_distributed.hh
+++ b/src/solver/solver_vector_distributed.hh
@@ -1,68 +1,63 @@
 /**
  * @file   solver_vector_distributed.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Feb 21 2013
  * @date last modification: Tue Jan 01 2019
  *
  * @brief  Solver vector interface for distributed arrays
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-2021 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 "solver_vector_default.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_
 #define AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_
 
 namespace akantu {
 
 class SolverVectorDistributed : public SolverVectorDefault {
 public:
   SolverVectorDistributed(DOFManagerDefault & dof_manager,
                           const ID & id = "solver_vector_mumps");
 
   SolverVectorDistributed(const SolverVectorDefault & vector,
                           const ID & id = "solver_vector_mumps");
 
   Array<Real> & getGlobalVector() override;
   void setGlobalVector(const Array<Real> & solution) override;
 
+  bool isFinite() const override;
 
-  virtual
-  bool
-  isDistributed()
-    const
-    override
-      { return true; }
-
+  virtual bool isDistributed() const override { return true; }
 
 protected:
   // full vector in case it needs to be centralized on master
   std::unique_ptr<Array<Real>> global_vector;
 };
 
 } // namespace akantu
 
 #endif /* AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_ */
diff --git a/src/solver/sparse_matrix.hh b/src/solver/sparse_matrix.hh
index 3c27783fc..1be77de04 100644
--- a/src/solver/sparse_matrix.hh
+++ b/src/solver/sparse_matrix.hh
@@ -1,168 +1,171 @@
 /**
  * @file   sparse_matrix.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Mon Dec 13 2010
  * @date last modification: Wed Sep 16 2020
  *
  * @brief  sparse matrix storage class (distributed assembled matrix)
  * This is a COO format (Coordinate List)
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 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"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SPARSE_MATRIX_HH_
 #define AKANTU_SPARSE_MATRIX_HH_
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 class DOFManager;
 class TermsToAssemble;
 class SolverVector;
 } // namespace akantu
 
 namespace akantu {
 
 class SparseMatrix {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   SparseMatrix(DOFManager & dof_manager, const MatrixType & matrix_type,
                const ID & id = "sparse_matrix");
 
   SparseMatrix(const SparseMatrix & matrix, const ID & id = "sparse_matrix");
 
   virtual ~SparseMatrix();
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// remove the existing profile
   virtual void clearProfile();
 
   /// set the matrix to 0
   virtual void set(Real val) = 0;
 
   virtual void zero() { this->set(0); }
 
   /// add a non-zero element to the profile
   virtual UInt add(UInt i, UInt j) = 0;
 
   /// assemble a local matrix in the sparse one
   virtual void add(UInt i, UInt j, Real value) = 0;
 
   /// save the profil in a file using the MatrixMarket file format
   virtual void saveProfile(const std::string & /* filename */) const {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// save the matrix in a file using the MatrixMarket file format
   virtual void saveMatrix(const std::string & /* filename */) const {
     AKANTU_TO_IMPLEMENT();
   };
 
   /// multiply the matrix by a coefficient
   virtual void mul(Real alpha) = 0;
 
   /// add matrices
   virtual void add(const SparseMatrix & B, Real alpha = 1.);
 
   /// Equivalent of *gemv in blas
   virtual void matVecMul(const SolverVector & x, SolverVector & y,
                          Real alpha = 1., Real beta = 0.) const = 0;
 
   /// modify the matrix to "remove" the blocked dof
   virtual void applyBoundary(Real block_val = 1.) = 0;
 
   /// copy the profile of another matrix
   virtual void copyProfile(const SparseMatrix & other) = 0;
 
   /// operator *=
   SparseMatrix & operator*=(Real alpha) {
     this->mul(alpha);
     return *this;
   }
 
+  /// Check if all entries are finite. The default implementation throws.
+  virtual bool isFinite() const { AKANTU_TO_IMPLEMENT(); }
+
 protected:
   /// This is the revert of add \f[B += \alpha * *this\f];
   virtual void addMeTo(SparseMatrix & B, Real alpha) const = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// return the values at potition i, j
   virtual inline Real operator()(UInt /*i*/, UInt /*j*/) const {
     AKANTU_TO_IMPLEMENT();
   }
   /// return the values at potition i, j
   virtual inline Real & operator()(UInt /*i*/, UInt /*j*/) {
     AKANTU_TO_IMPLEMENT();
   }
 
   /// return the minimum value
   virtual inline Real min() { AKANTU_TO_IMPLEMENT(); }
 
   AKANTU_GET_MACRO(NbNonZero, nb_non_zero, UInt);
   UInt size() const { return size_; }
   AKANTU_GET_MACRO(MatrixType, matrix_type, const MatrixType &);
 
   virtual UInt getRelease() const = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   ID id;
 
   /// Underlying dof manager
   DOFManager & _dof_manager;
 
   /// sparce matrix type
   MatrixType matrix_type;
 
   /// Size of the matrix
   UInt size_;
 
   /// number of processors
   UInt nb_proc;
 
   /// number of non zero element
   UInt nb_non_zero;
 };
 
 // Array<Real> & operator*=(Array<Real> & vect, const SparseMatrix & mat);
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "sparse_matrix_inline_impl.hh"
 
 #endif /* AKANTU_SPARSE_MATRIX_HH_ */
diff --git a/src/solver/sparse_matrix_aij.cc b/src/solver/sparse_matrix_aij.cc
index cb4db6eb2..526449aa1 100644
--- a/src/solver/sparse_matrix_aij.cc
+++ b/src/solver/sparse_matrix_aij.cc
@@ -1,301 +1,305 @@
 /**
  * @file   sparse_matrix_aij.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Aug 21 2015
  * @date last modification: Fri Jul 24 2020
  *
  * @brief  Implementation of the AIJ sparse matrix
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2015-2021 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 "sparse_matrix_aij.hh"
 #include "aka_iterators.hh"
+#include "aka_math.hh"
 #include "dof_manager_default.hh"
 #include "dof_synchronizer.hh"
 #include "solver_vector_default.hh"
 #include "terms_to_assemble.hh"
 /* -------------------------------------------------------------------------- */
 #include <fstream>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 SparseMatrixAIJ::SparseMatrixAIJ(DOFManagerDefault & dof_manager,
                                  const MatrixType & matrix_type, const ID & id)
     : SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager),
       irn(0, 1, id + ":irn"), jcn(0, 1, id + ":jcn"), a(0, 1, id + ":a") {}
 
 /* -------------------------------------------------------------------------- */
 SparseMatrixAIJ::SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id)
     : SparseMatrix(matrix, id), dof_manager(matrix.dof_manager),
       irn(matrix.irn, id + ":irn"), jcn(matrix.jcn, id + ":jcn"),
       a(matrix.a, id + ":a") {}
 
 /* -------------------------------------------------------------------------- */
 SparseMatrixAIJ::~SparseMatrixAIJ() = default;
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::applyBoundary(Real block_val) {
   AKANTU_DEBUG_IN();
 
   const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs();
   auto begin = blocked_dofs.begin();
   auto end = blocked_dofs.end();
 
   auto is_blocked = [&](auto && i) -> bool {
     auto il = this->dof_manager.globalToLocalEquationNumber(i);
     return std::binary_search(begin, end, il);
   };
 
   for (auto && ij_a : zip(irn, jcn, a)) {
     UInt ni = std::get<0>(ij_a) - 1;
     UInt nj = std::get<1>(ij_a) - 1;
 
     if (is_blocked(ni) or is_blocked(nj)) {
 
       std::get<2>(ij_a) =
           std::get<0>(ij_a) != std::get<1>(ij_a) ? 0.
           : this->dof_manager.isLocalOrMasterDOF(
                 this->dof_manager.globalToLocalEquationNumber(ni))
               ? block_val
               : 0.;
     }
   }
 
   this->value_release++;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::saveProfile(const std::string & filename) const {
   AKANTU_DEBUG_IN();
 
   std::ofstream outfile;
   outfile.open(filename.c_str());
 
   UInt m = this->size_;
 
   auto & comm = dof_manager.getCommunicator();
 
   // write header
   if (comm.whoAmI() == 0) {
 
     outfile << "%%MatrixMarket matrix coordinate pattern";
     if (this->matrix_type == _symmetric) {
       outfile << " symmetric";
     } else {
       outfile << " general";
     }
     outfile << std::endl;
     outfile << m << " " << m << " " << this->nb_non_zero << std::endl;
   }
 
   for (auto p : arange(comm.getNbProc())) {
     // write content
     if (comm.whoAmI() == p) {
       for (UInt i = 0; i < this->nb_non_zero; ++i) {
         outfile << this->irn.storage()[i] << " " << this->jcn.storage()[i]
                 << " 1" << std::endl;
       }
     }
     comm.barrier();
   }
 
   outfile.close();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::saveMatrix(const std::string & filename) const {
   AKANTU_DEBUG_IN();
   auto & comm = dof_manager.getCommunicator();
 
   // open and set the properties of the stream
   std::ofstream outfile;
 
   if (0 == comm.whoAmI()) {
     outfile.open(filename.c_str());
   } else {
     outfile.open(filename.c_str(), std::ios_base::app);
   }
 
   outfile.precision(std::numeric_limits<Real>::digits10);
   // write header
   decltype(nb_non_zero) nnz = this->nb_non_zero;
   comm.allReduce(nnz);
 
   if (comm.whoAmI() == 0) {
     outfile << "%%MatrixMarket matrix coordinate real";
     if (this->matrix_type == _symmetric) {
       outfile << " symmetric";
     } else {
       outfile << " general";
     }
     outfile << std::endl;
     outfile << this->size_ << " " << this->size_ << " " << nnz << std::endl;
   }
 
   for (auto p : arange(comm.getNbProc())) {
     // write content
     if (comm.whoAmI() == p) {
       for (UInt i = 0; i < this->nb_non_zero; ++i) {
         outfile << this->irn(i) << " " << this->jcn(i) << " " << this->a(i)
                 << std::endl;
       }
     }
     comm.barrier();
   }
   // time to end
   outfile.close();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::matVecMul(const Array<Real> & x, Array<Real> & y,
                                 Real alpha, Real beta) const {
   AKANTU_DEBUG_IN();
 
   y *= beta;
 
   auto i_it = this->irn.begin();
   auto j_it = this->jcn.begin();
   auto a_it = this->a.begin();
   auto a_end = this->a.end();
   auto x_it = x.begin_reinterpret(x.size() * x.getNbComponent());
   auto y_it = y.begin_reinterpret(x.size() * x.getNbComponent());
 
   for (; a_it != a_end; ++i_it, ++j_it, ++a_it) {
     Int i = this->dof_manager.globalToLocalEquationNumber(*i_it - 1);
     Int j = this->dof_manager.globalToLocalEquationNumber(*j_it - 1);
     const Real & A = *a_it;
 
     y_it[i] += alpha * A * x_it[j];
 
     if ((this->matrix_type == _symmetric) && (i != j)) {
       y_it[j] += alpha * A * x_it[i];
     }
   }
 
   if (this->dof_manager.hasSynchronizer()) {
     this->dof_manager.getSynchronizer().reduceSynchronizeArray<AddOperation>(y);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::matVecMul(const SolverVector & _x, SolverVector & _y,
                                 Real alpha, Real beta) const {
   AKANTU_DEBUG_IN();
 
   auto && x = aka::as_type<SolverVectorArray>(_x).getVector();
   auto && y = aka::as_type<SolverVectorArray>(_y).getVector();
   this->matVecMul(x, y, alpha, beta);
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::copyContent(const SparseMatrix & matrix) {
   AKANTU_DEBUG_IN();
   const auto & mat = aka::as_type<SparseMatrixAIJ>(matrix);
   AKANTU_DEBUG_ASSERT(nb_non_zero == mat.getNbNonZero(),
                       "The to matrix don't have the same profiles");
   memcpy(a.storage(), mat.getA().storage(), nb_non_zero * sizeof(Real));
 
   this->value_release++;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::copyProfile(const SparseMatrix & other) {
   const auto & A = aka::as_type<SparseMatrixAIJ>(other);
 
   SparseMatrix::clearProfile();
 
   this->irn.copy(A.irn);
   this->jcn.copy(A.jcn);
 
   this->irn_jcn_k.clear();
 
   UInt i;
   UInt j;
   UInt k;
   for (auto && data : enumerate(irn, jcn)) {
     std::tie(k, i, j) = data;
 
     this->irn_jcn_k[this->key(i - 1, j - 1)] = k;
   }
 
   this->nb_non_zero = this->irn.size();
   this->a.resize(this->nb_non_zero);
 
   this->a.set(0.);
   this->size_ = A.size_;
 
   this->profile_release = A.profile_release;
   this->value_release++;
 }
 
 /* -------------------------------------------------------------------------- */
 template <class MatrixType>
 void SparseMatrixAIJ::addMeToTemplated(MatrixType & B, Real alpha) const {
   UInt i;
   UInt j;
   Real A_ij;
   for (auto && tuple : zip(irn, jcn, a)) {
     std::tie(i, j, A_ij) = tuple;
     B.add(i - 1, j - 1, alpha * A_ij);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::addMeTo(SparseMatrix & B, Real alpha) const {
 
   if (aka::is_of_type<SparseMatrixAIJ>(B)) {
     this->addMeToTemplated<SparseMatrixAIJ>(aka::as_type<SparseMatrixAIJ>(B),
                                             alpha);
   } else {
     //    this->addMeToTemplated<SparseMatrix>(*this, alpha);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::mul(Real alpha) {
   this->a *= alpha;
   this->value_release++;
 }
 
 /* -------------------------------------------------------------------------- */
 void SparseMatrixAIJ::set(Real val) {
   a.set(val);
 
   this->value_release++;
 }
 
+/* -------------------------------------------------------------------------- */
+bool SparseMatrixAIJ::isFinite() const { return this->a.isFinite(); };
+
 } // namespace akantu
diff --git a/src/solver/sparse_matrix_aij.hh b/src/solver/sparse_matrix_aij.hh
index ebcbd1dee..771d48e8e 100644
--- a/src/solver/sparse_matrix_aij.hh
+++ b/src/solver/sparse_matrix_aij.hh
@@ -1,203 +1,206 @@
 /**
  * @file   sparse_matrix_aij.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Mon Dec 13 2010
  * @date last modification: Wed Sep 16 2020
  *
  * @brief  AIJ implementation of the SparseMatrix (this the format used by
  * Mumps)
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 #include "aka_common.hh"
 #include "sparse_matrix.hh"
 /* -------------------------------------------------------------------------- */
 #include <unordered_map>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SPARSE_MATRIX_AIJ_HH_
 #define AKANTU_SPARSE_MATRIX_AIJ_HH_
 
 namespace akantu {
 class DOFManagerDefault;
 class TermsToAssemble;
 } // namespace akantu
 
 namespace akantu {
 
 class SparseMatrixAIJ : public SparseMatrix {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   SparseMatrixAIJ(DOFManagerDefault & dof_manager,
                   const MatrixType & matrix_type,
                   const ID & id = "sparse_matrix_aij");
 
   SparseMatrixAIJ(const SparseMatrixAIJ & matrix,
                   const ID & id = "sparse_matrix_aij");
 
   ~SparseMatrixAIJ() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// remove the existing profile
   inline void clearProfile() override;
 
   /// add a non-zero element
   inline UInt add(UInt i, UInt j) override;
 
   /// set the matrix to 0
   void set(Real val) override;
 
   /// assemble a local matrix in the sparse one
   inline void add(UInt i, UInt j, Real value) override;
 
   /// add a block of values
   inline void addValues(const Vector<Int> & is, const Vector<Int> & js,
                         const Matrix<Real> & values, MatrixType values_type);
 
   /// set the size of the matrix
   void resize(UInt size) { this->size_ = size; }
 
   /// modify the matrix to "remove" the blocked dof
   void applyBoundary(Real block_val = 1.) override;
 
   /// save the profil in a file using the MatrixMarket file format
   void saveProfile(const std::string & filename) const override;
 
   /// save the matrix in a file using the MatrixMarket file format
   void saveMatrix(const std::string & filename) const override;
 
   /// copy assuming the profile are the same
   virtual void copyContent(const SparseMatrix & matrix);
 
   /// multiply the matrix by a scalar
   void mul(Real alpha) override;
 
   /// Equivalent of *gemv in blas
   void matVecMul(const SolverVector & x, SolverVector & y, Real alpha = 1.,
                  Real beta = 0.) const override;
 
   void matVecMul(const Array<Real> & x, Array<Real> & y, Real alpha = 1.,
                  Real beta = 0.) const;
 
   /// copy the profile of another matrix
   void copyProfile(const SparseMatrix & other) override;
 
+  /// Check if all entries are finite
+  virtual bool isFinite() const override;
+
   /* ------------------------------------------------------------------------ */
   /// accessor to A_{ij} - if (i, j) not present it returns 0
   inline Real operator()(UInt i, UInt j) const override;
 
   /// accessor to A_{ij} - if (i, j) not present it fails, (i, j) should be
   /// first added to the profile
   inline Real & operator()(UInt i, UInt j) override;
 
   /// accessor to get the minimum value of A_{ij}
   inline Real min() override;
 
 protected:
   void addMeTo(SparseMatrix & B, Real alpha) const override;
 
   inline void addSymmetricValuesToSymmetric(const Vector<Int> & is,
                                             const Vector<Int> & js,
                                             const Matrix<Real> & values);
   inline void addUnsymmetricValuesToSymmetric(const Vector<Int> & is,
                                               const Vector<Int> & js,
                                               const Matrix<Real> & values);
   inline void addValuesToUnsymmetric(const Vector<Int> & is,
                                      const Vector<Int> & js,
                                      const Matrix<Real> & values);
 
 private:
   /// This is just to inline the addToMatrix function
   template <class MatrixType>
   void addMeToTemplated(MatrixType & B, Real alpha) const;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   AKANTU_GET_MACRO(IRN, irn, const Array<Int> &);
 
   AKANTU_GET_MACRO(JCN, jcn, const Array<Int> &);
 
   AKANTU_GET_MACRO(A, a, const Array<Real> &);
 
   /// The release changes at each call of a function that changes the profile,
   /// it in increasing but could overflow so it should be checked as
   /// (my_release != release) and not as (my_release < release)
   AKANTU_GET_MACRO(ProfileRelease, profile_release, UInt);
   AKANTU_GET_MACRO(ValueRelease, value_release, UInt);
   UInt getRelease() const override { return value_release; }
 
 protected:
   using KeyCOO = std::pair<UInt, UInt>;
   using coordinate_list_map = std::unordered_map<KeyCOO, UInt>;
 
   /// get the pair corresponding to (i, j)
   inline KeyCOO key(UInt i, UInt j) const {
     if (this->matrix_type == _symmetric && (i > j)) {
       return std::make_pair(j, i);
     }
     return std::make_pair(i, j);
   }
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   DOFManagerDefault & dof_manager;
 
   /// row indexes
   Array<Int> irn;
 
   /// column indexes
   Array<Int> jcn;
 
   /// values : A[k] = Matrix[irn[k]][jcn[k]]
   Array<Real> a;
 
   /// Profile release
   UInt profile_release{1};
 
   /// Value release
   UInt value_release{1};
 
   /// map for (i, j) ->  k correspondence
   coordinate_list_map irn_jcn_k;
 };
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "sparse_matrix_aij_inline_impl.hh"
 
 #endif /* AKANTU_SPARSE_MATRIX_AIJ_HH_ */