diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 703ad6391..76f8580d0 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -1,110 +1,111 @@
 
 #===============================================================================
 # @file   CMakeLists.txt
 #
 # @author Nicolas Richart <nicolas.richart@epfl.ch>
 #
 # @date creation: Fri Dec 12 2014
 # @date last modification: Mon Jan 18 2016
 #
 # @brief  CMake file for the python wrapping of akantu
 #
 # @section LICENSE
 #
 # Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory
 # (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 #
 # Akantu is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Akantu. If not, see <http://www.gnu.org/licenses/>.
 #
 #===============================================================================
 
 package_get_all_include_directories(
   AKANTU_LIBRARY_INCLUDE_DIRS
   )
 
 package_get_all_external_informations(
   PRIVATE_INCLUDE AKANTU_PRIVATE_EXTERNAL_INCLUDE_DIR
   INTERFACE_INCLUDE AKANTU_INTERFACE_EXTERNAL_INCLUDE_DIR
   LIBRARIES AKANTU_EXTERNAL_LIBRARIES
   )
 
 
 set(PYAKANTU_SRCS
   py_aka_common.cc
   py_aka_error.cc
   py_akantu.cc
   py_boundary_conditions.cc
   py_fe_engine.cc
   py_group_manager.cc
   py_mesh.cc
   py_model.cc
   py_parser.cc
+  py_solver.cc
   )
 
 package_is_activated(iohelper _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_dumpable.cc
     )
 endif()
 
 
 package_is_activated(solid_mechanics _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_solid_mechanics_model.cc
     py_material.cc
     )
 endif()
 
 package_is_activated(cohesive_element _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_solid_mechanics_model_cohesive.cc
     )
 endif()
 
 package_is_activated(heat_transfer _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_heat_transfer_model.cc
     )
 endif()
 
 package_is_activated(structural_mechanics _is_activated)
 if (_is_activated)
   list(APPEND PYAKANTU_SRCS
     py_structural_mechanics_model.cc
     )
 endif()
 
 
 pybind11_add_module(py11_akantu ${PYAKANTU_SRCS})
 target_link_libraries(py11_akantu PUBLIC akantu)
 set_target_properties(py11_akantu PROPERTIES DEBUG_POSTFIX "")
 
 set(_python_install_dir
   ${CMAKE_INSTALL_LIBDIR}/python${PYTHON_VERSION_MAJOR}.${PYTHON_VERSION_MINOR}/site-packages)
 
 install(TARGETS py11_akantu
   LIBRARY DESTINATION ${_python_install_dir})
 
 install(DIRECTORY akantu
   DESTINATION ${_python_install_dir}
   FILES_MATCHING PATTERN "*.py")
 
 
 configure_file(setup.py.in
   "${CMAKE_CURRENT_BINARY_DIR}/setup.py" @ONLY)
 
diff --git a/python/py_akantu.cc b/python/py_akantu.cc
index 87c9e1fe5..d3d3a93a8 100644
--- a/python/py_akantu.cc
+++ b/python/py_akantu.cc
@@ -1,111 +1,113 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_config.hh"
 /* -------------------------------------------------------------------------- */
 #include "py_aka_common.hh"
 #include "py_aka_error.hh"
 #include "py_boundary_conditions.hh"
 #include "py_fe_engine.hh"
 #include "py_group_manager.hh"
 #include "py_mesh.hh"
 #include "py_model.hh"
 #include "py_parser.hh"
+#include "py_solver.hh"
 
 #if defined(AKANTU_USE_IOHELPER)
 #include "py_dumpable.hh"
 #endif
 
 #if defined(AKANTU_SOLID_MECHANICS)
 #include "py_material.hh"
 #include "py_solid_mechanics_model.hh"
 #endif
 
 #if defined(AKANTU_HEAT_TRANSFER)
 #include "py_heat_transfer_model.hh"
 #endif
 
 #if defined(AKANTU_COHESIVE_ELEMENT)
 #include "py_solid_mechanics_model_cohesive.hh"
 #endif
 
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
 #include "py_structural_mechanics_model.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 #include <aka_error.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 
 namespace py = pybind11;
 
 namespace akantu {
 void register_all(pybind11::module & mod) {
   register_initialize(mod);
   register_enums(mod);
   register_error(mod);
   register_functions(mod);
   register_parser(mod);
+  register_solvers(mod);
 
   register_group_manager(mod);
 #if defined(AKANTU_USE_IOHELPER)
   register_dumpable(mod);
 #endif
   register_mesh(mod);
 
   register_fe_engine(mod);
 
   register_boundary_conditions(mod);
   register_model(mod);
 #if defined(AKANTU_HEAT_TRANSFER)
   register_heat_transfer_model(mod);
 #endif
 
 #if defined(AKANTU_SOLID_MECHANICS)
   register_solid_mechanics_model(mod);
   register_material(mod);
 #endif
 
 #if defined(AKANTU_COHESIVE_ELEMENT)
   register_solid_mechanics_model_cohesive(mod);
 #endif
 
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
   register_structural_mechanics_model(mod);
 #endif
 }
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 PYBIND11_MODULE(py11_akantu, mod) {
   mod.doc() = "Akantu python interface";
 
   static py::exception<akantu::debug::Exception>
       akantu_exception(mod, "Exception");
 
   py::register_exception_translator([](std::exception_ptr ptr) {
     try {
       if (ptr) {
         std::rethrow_exception(ptr);
       }
     } catch (akantu::debug::Exception & e) {
       if (akantu::debug::debugger.printBacktrace()) {
         akantu::debug::printBacktrace();
       }
       akantu_exception(e.info().c_str());
     }
   });
 
   akantu::register_all(mod);
 
   mod.def("has_mpi", []() {
 #if defined(AKANTU_USE_MPI)
     return true;
 #else
     return false;
 #endif
   });
 
 } // Module akantu
diff --git a/python/py_mesh.cc b/python/py_mesh.cc
index be3c85931..5ed06f8df 100644
--- a/python/py_mesh.cc
+++ b/python/py_mesh.cc
@@ -1,126 +1,141 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_config.hh"
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <mesh.hh>
-#include <mesh_utils.hh>
 #include <mesh_accessor.hh>
+#include <mesh_utils.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
-
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void register_element_type_map_array(py::module & mod,
                                      const std::string & name) {
   py::class_<ElementTypeMapArray<T>, std::shared_ptr<ElementTypeMapArray<T>>>(
       mod, ("ElementTypeMapArray" + name).c_str())
       .def(
           "__call__",
           [](ElementTypeMapArray<T> & self, ElementType type,
              GhostType ghost_type) -> decltype(auto) {
             return self(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost,
           py::return_value_policy::reference)
       .def(
           "elementTypes",
           [](ElementTypeMapArray<T> & self, UInt _dim, GhostType _ghost_type,
              ElementKind _kind) -> decltype(auto) {
             auto types = self.elementTypes(_dim, _ghost_type, _kind);
             std::vector<ElementType> _types;
             for (auto && t : types) {
               _types.push_back(t);
             }
             return _types;
           },
           py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost,
           py::arg("kind") = _ek_regular);
 }
 
-
 /* -------------------------------------------------------------------------- */
 void register_mesh(py::module & mod) {
 
   register_element_type_map_array<Real>(mod, "Real");
   register_element_type_map_array<UInt>(mod, "UInt");
 
   py::class_<MeshData>(mod, "MeshData")
       .def(
           "getElementalDataUInt",
-          [](MeshData & _this, const ID & name) -> ElementTypeMapArray<UInt> & {
+          [](MeshData & _this, const ID & name) -> decltype(auto) {
             return _this.getElementalData<UInt>(name);
           },
+          py::return_value_policy::reference)
+      .def(
+          "getElementalDataReal",
+          [](MeshData & _this, const ID & name) -> decltype(auto) {
+            return _this.getElementalData<Real>(name);
+          },
           py::return_value_policy::reference);
 
   py::class_<Mesh, GroupManager, Dumpable, MeshData>(mod, "Mesh",
                                                      py::multiple_inheritance())
       .def(py::init<UInt, const ID &, const MemoryID &>(),
            py::arg("spatial_dimension"), py::arg("id") = "mesh",
            py::arg("memory_id") = 0)
       .def("read", &Mesh::read, py::arg("filename"),
            py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file")
       .def(
           "getNodes",
           [](Mesh & self) -> decltype(auto) { return self.getNodes(); },
           py::return_value_policy::reference)
       .def("getNbNodes", &Mesh::getNbNodes)
       .def(
           "getConnectivity",
           [](Mesh & self, ElementType type) -> decltype(auto) {
             return self.getConnectivity(type);
           },
           py::return_value_policy::reference)
-      .def("addConnectivityType",
+      .def(
+          "addConnectivityType",
           [](Mesh & self, ElementType type, GhostType ghost_type) -> void {
             self.addConnectivityType(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def("distribute", [](Mesh & self) { self.distribute(); })
       .def("makePeriodic",
            [](Mesh & self, const SpatialDirection & direction) {
              self.makePeriodic(direction);
            })
       .def(
           "getNbElement",
-          [](Mesh & self, const UInt spatial_dimension,
-             GhostType ghost_type, ElementKind kind) {
+          [](Mesh & self, const UInt spatial_dimension, GhostType ghost_type,
+             ElementKind kind) {
             return self.getNbElement(spatial_dimension, ghost_type, kind);
           },
           py::arg("spatial_dimension") = _all_dimensions,
           py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined)
       .def(
           "getNbElement",
-          [](Mesh & self, ElementType type,
-             GhostType ghost_type) {
+          [](Mesh & self, ElementType type, GhostType ghost_type) {
             return self.getNbElement(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost)
-      .def_static("getSpatialDimension", [](ElementType & type) {
-        return Mesh::getSpatialDimension(type);
-      });
+      .def_static(
+          "getSpatialDimension",
+          [](ElementType & type) { return Mesh::getSpatialDimension(type); })
+      .def(
+          "getDataReal",
+          [](Mesh & _this, const ID & name, ElementType type,
+             GhostType ghost_type) -> decltype(auto) {
+            return _this.getData<Real>(name, type, ghost_type);
+          },
+          py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost,
+          py::return_value_policy::reference);
 
   /* ------------------------------------------------------------------------ */
   py::class_<MeshUtils>(mod, "MeshUtils")
       .def_static("buildFacets", &MeshUtils::buildFacets);
 
   py::class_<MeshAccessor>(mod, "MeshAccessor")
+      .def(py::init<Mesh &>(), py::arg("mesh"))
+      .def(
+          "resizeConnectivity",
+          [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt)
+              -> void { self.resizeConnectivity(new_size, type, gt); },
+          py::arg("new_size"), py::arg("type"),
+          py::arg("ghost_type") = _not_ghost)
       .def(
-      	  py::init<Mesh &>(), py::arg("mesh"))
-      .def("resizeConnectivity",
-      	   [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt) -> void {
-      	     self.resizeConnectivity(new_size, type, gt);},
-      	   py::arg("new_size"), py::arg("type"), py::arg("ghost_type") = _not_ghost)
-      .def("resizeNodes",
-      	   [](MeshAccessor& self, UInt new_size) -> void {
-      	     self.resizeNodes(new_size);},
-           py::arg("new_size"))
+          "resizeNodes",
+          [](MeshAccessor & self, UInt new_size) -> void {
+            self.resizeNodes(new_size);
+          },
+          py::arg("new_size"))
       .def("makeReady", &MeshAccessor::makeReady);
 }
 } // namespace akantu
diff --git a/python/py_model.cc b/python/py_model.cc
index 8b461ee07..6863419f4 100644
--- a/python/py_model.cc
+++ b/python/py_model.cc
@@ -1,82 +1,83 @@
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <model.hh>
 #include <non_linear_solver.hh>
 #include <sparse_matrix_aij.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 void register_model(py::module & mod) {
   py::class_<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("assembleToResidual", &DOFManager::assembleToResidual);
 
   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_<ModelSolver, Parsable>(mod, "ModelSolver",
                                     py::multiple_inheritance())
       .def("getNonLinearSolver",
            (NonLinearSolver & (ModelSolver::*)(const ID &)) &
                ModelSolver::getNonLinearSolver,
            py::arg("solver_id") = "", py::return_value_policy::reference)
       .def("solveStep", [](ModelSolver & self) { self.solveStep(); })
       .def("solveStep", [](ModelSolver & self, const ID & solver_id) {
         self.solveStep(solver_id);
       });
 
   py::class_<Model, ModelSolver>(mod, "Model", py::multiple_inheritance())
       .def("setBaseName", &Model::setBaseName)
       .def("getFEEngine", &Model::getFEEngine, py::arg("name") = "",
            py::return_value_policy::reference)
       .def("getFEEngineBoundary", &Model::getFEEngine, py::arg("name") = "",
            py::return_value_policy::reference)
       .def("addDumpFieldVector", &Model::addDumpFieldVector)
       .def("addDumpField", &Model::addDumpField)
       .def("setBaseNameToDumper", &Model::setBaseNameToDumper)
       .def("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper)
       .def("addDumpFieldToDumper", &Model::addDumpFieldToDumper)
       .def("dump", py::overload_cast<>(&Model::dump))
       .def("dump", py::overload_cast<UInt>(&Model::dump))
       .def("dump", py::overload_cast<Real, UInt>(&Model::dump))
       .def("initNewSolver", &Model::initNewSolver)
       .def("getDOFManager", &Model::getDOFManager,
-           py::return_value_policy::reference);
+           py::return_value_policy::reference)
+      .def("assembleMatrix", &Model::assembleMatrix);
 }
 
 } // namespace akantu
diff --git a/python/py_solver.cc b/python/py_solver.cc
index 6e292867d..7ac4674e2 100644
--- a/python/py_solver.cc
+++ b/python/py_solver.cc
@@ -1,45 +1,47 @@
 /* -------------------------------------------------------------------------- */
+#include "py_solver.hh"
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
-#include "py_solver.h"
-/* -------------------------------------------------------------------------- */#include <model.hh>
+#include <model.hh>
 #include <non_linear_solver.hh>
 #include <sparse_matrix_aij.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
-void register_solver(py::module & mod) {
+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("__call__", [](const SparseMatrix & self, UInt i, UInt j) {
         return self(i, j);
       });
 
   py::class_<SparseMatrixAIJ, SparseMatrix>(mod, "SparseMatrixAIJ")
       .def("getIRN", &SparseMatrixAIJ::getIRN)
       .def("getJCN", &SparseMatrixAIJ::getJCN)
       .def("getA", &SparseMatrixAIJ::getA);
 
   py::class_<SolverVector>(mod, "SolverVector");
 }
+
+} // namespace akantu
diff --git a/python/py_solver.hh b/python/py_solver.hh
index 892606998..41e1124e8 100644
--- a/python/py_solver.hh
+++ b/python/py_solver.hh
@@ -1,12 +1,12 @@
 #include <pybind11/pybind11.h>
 
 #ifndef AKANTU_PY_AKA_SOLVER_HH_
 #define AKANTU_PY_AKA_SOLVER_HH_
 
 namespace akantu {
 
-void register_solver(pybind11::module & mod);
+void register_solvers(pybind11::module & mod);
 
 }
 
 #endif
diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc
index 2c6a9d885..3098b1b22 100644
--- a/src/model/structural_mechanics/structural_mechanics_model.cc
+++ b/src/model/structural_mechanics/structural_mechanics_model.cc
@@ -1,479 +1,483 @@
 /**
  * @file   structural_mechanics_model.cc
  *
  * @author Fabian Barras <fabian.barras@epfl.ch>
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  * @author Sébastien Hartmann <sebastien.hartmann@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Damien Spielmann <damien.spielmann@epfl.ch>
  *
  * @date creation: Fri Jul 15 2011
  * @date last modification: Wed Feb 21 2018
  *
  * @brief  Model implementation for Structural Mechanics elements
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "structural_mechanics_model.hh"
 #include "dof_manager.hh"
 #include "integrator_gauss.hh"
 #include "mesh.hh"
 #include "shape_structural.hh"
 #include "sparse_matrix.hh"
 #include "time_step_solver.hh"
 /* -------------------------------------------------------------------------- */
 #ifdef AKANTU_USE_IOHELPER
 #include "dumpable_inline_impl.hh"
 #include "dumper_elemental_field.hh"
 #include "dumper_internal_material_field.hh"
 #include "dumper_iohelper_paraview.hh"
 #include "group_manager_inline_impl.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 #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,
                                                    const MemoryID & memory_id)
     : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id),
       f_m2a(1.0), stress("stress", id, memory_id),
       element_material("element_material", id, memory_id),
       set_ID("beam sets", id, memory_id),
       rotation_matrix("rotation_matices", id, memory_id) {
   AKANTU_DEBUG_IN();
 
   registerFEEngineObject<MyFEEngineType>("StructuralMechanicsFEEngine", mesh,
                                          spatial_dimension);
 
   if (spatial_dimension == 2) {
     nb_degree_of_freedom = 3;
   } else if (spatial_dimension == 3) {
     nb_degree_of_freedom = 6;
   } else {
     AKANTU_TO_IMPLEMENT();
   }
 
 #ifdef AKANTU_USE_IOHELPER
   this->mesh.registerDumper<DumperParaview>("structural_mechanics_model", id,
                                             true);
 #endif
   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);
 #if defined(AKANTU_USE_IOHELPER)
   this->mesh.getDumper().setTimeStep(time_step);
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 /* 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);
     }
   }
 
   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();
 }
 
 /* -------------------------------------------------------------------------- */
 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;
 
   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, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "velocity") {
     return mesh.createStridedNodalField(velocity, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "acceleration") {
     return mesh.createStridedNodalField(acceleration, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "rotation") {
     return mesh.createStridedNodalField(displacement_rotation, group_name,
                                         nb_degree_of_freedom - n, n,
                                         padding_size);
   }
 
   if (field_name == "force") {
     return mesh.createStridedNodalField(external_force, group_name, n, 0,
                                         padding_size);
   }
   if (field_name == "external_force") {
     return mesh.createStridedNodalField(external_force, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "momentum") {
     return mesh.createStridedNodalField(
         external_force, group_name, nb_degree_of_freedom - n, n, padding_size);
   }
 
   if (field_name == "internal_force") {
     return mesh.createStridedNodalField(internal_force, group_name, n, 0,
                                         padding_size);
   }
 
   if (field_name == "internal_momentum") {
     return mesh.createStridedNodalField(
         internal_force, group_name, nb_degree_of_freedom - n, n, 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*/) {
   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*/) {}
 
 /// callback to assemble the residual StructuralMechanicsModel::(rhs)
 void StructuralMechanicsModel::assembleResidual() {
   AKANTU_DEBUG_IN();
 
   auto & dof_manager = getDOFManager();
 
   assembleInternalForce();
 
   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();
 
   if ("external" == residual_part) {
     this->getDOFManager().assembleToResidual("displacement",
                                              *this->external_force, 1);
     AKANTU_DEBUG_OUT();
     return;
   }
 
   if ("internal" == residual_part) {
     this->assembleInternalForce();
     this->getDOFManager().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);
   }
   default:
     return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
   }
 }
 
 /* ------------------------------------------------------------------------ */
 ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions(
     const TimeStepSolverType & type) const {
   ModelSolverOptions options;
 
   switch (type) {
   case TimeStepSolverType::_static: {
     options.non_linear_solver_type = NonLinearSolverType::_linear;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_pseudo_time;
     options.solution_type["displacement"] = IntegrationScheme::_not_defined;
     break;
   }
   case TimeStepSolverType::_dynamic: {
     options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_trapezoidal_rule_2;
     options.solution_type["displacement"] = IntegrationScheme::_displacement;
     break;
   }
   default:
     AKANTU_EXCEPTION(type << " is not a valid time step solver type");
   }
 
   return options;
 }
 
 /* -------------------------------------------------------------------------- */
 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);
   BtSigma.resize(0);
 
   getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force,
                                                    type, gt, -1.);
 }
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu