diff --git a/python/py_aka_array.hh b/python/py_aka_array.hh
index 6da82a07c..d70615d06 100644
--- a/python/py_aka_array.hh
+++ b/python/py_aka_array.hh
@@ -1,282 +1,287 @@
 /**
  * @file   py_aka_array.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Oct 31 2018
  * @date last modification: Fri Nov 13 2020
  *
  * @brief  pybind11 interface to akantu 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 <aka_array.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 
 namespace py = pybind11;
 namespace _aka = akantu;
 
 namespace akantu {
 namespace detail {
   template <class T> struct is_array_type : public std::false_type {};
   template <class T> struct is_array_type<Vector<T>> : public std::true_type {};
   template <class T> struct is_array_type<Matrix<T>> : public std::true_type {};
   template <class T> struct is_array_type<Array<T>> : public std::true_type {};
 
   /* ------------------------------------------------------------------------ */
   template <typename T> class ArrayProxy : public Array<T> {
   protected:
     // deallocate the memory
     void deallocate() final {}
 
     // allocate the memory
     void allocate(UInt /*size*/, UInt /*nb_component*/) final {}
 
     // allocate and initialize the memory
     void allocate(UInt /*size*/, UInt /*nb_component*/,
                   const T & /*value*/) final {}
 
   public:
     ArrayProxy(T * data, UInt size, UInt nb_component) {
       this->values = data;
       this->size_ = size;
       this->nb_component = nb_component;
     }
 
     ArrayProxy(const Array<T> & src) {
       this->values = src.storage();
       this->size_ = src.size();
       this->nb_component = src.getNbComponent();
     }
 
     ~ArrayProxy() override { this->values = nullptr; }
 
     void resize(UInt size, const T & /*val*/) final {
       if (size != this->size()) {
         AKANTU_EXCEPTION("cannot resize a temporary array");
       }
       // std::fill(this->begin(), this->end(), val);
     }
 
     void resize(UInt new_size) final {
       if (new_size != this->size()) {
         AKANTU_EXCEPTION("cannot resize a temporary array");
       }
     }
 
     void reserve(UInt /*size*/, UInt /*new_size*/) final {
       AKANTU_EXCEPTION("cannot resize a temporary array");
     }
   };
 
   /* ------------------------------------------------------------------------ */
   template <typename T> struct ProxyType {};
   template <typename T> struct ProxyType<Vector<T>> { using type = Vector<T>; };
   template <typename T> struct ProxyType<Matrix<T>> { using type = Matrix<T>; };
   template <typename T> struct ProxyType<Array<T>> {
     using type = ArrayProxy<T>;
   };
   template <typename array> using ProxyType_t = typename ProxyType<array>::type;
 } // namespace detail
 } // namespace akantu
 
 namespace pybind11 {
 namespace detail {
 
   template <typename T> struct AkaArrayType {
     using type =
         array_t<typename T::value_type, array::c_style | array::forcecast>;
   };
 
   template <typename T> struct AkaArrayType<_aka::Vector<T>> {
     using type = array_t<T, array::f_style | array::forcecast>;
   };
   template <typename T> struct AkaArrayType<_aka::Matrix<T>> {
     using type = array_t<T, array::f_style | array::forcecast>;
   };
 
   template <typename U> using array_type_t = typename AkaArrayType<U>::type;
 
   /* ------------------------------------------------------------------------ */
   template <typename T>
   decltype(auto) create_proxy(array_type_t<_aka::Vector<T>> & ref,
                               const _aka::Vector<T> * /*unused*/) {
     return std::make_unique<_aka::detail::ProxyType_t<_aka::Vector<T>>>(
         ref.mutable_data(), ref.shape(0));
   }
 
   template <typename T>
   decltype(auto) create_proxy(array_type_t<_aka::Matrix<T>> & ref,
                               const _aka::Matrix<T> * /*unused*/) {
     return std::make_unique<_aka::detail::ProxyType_t<_aka::Matrix<T>>>(
         ref.mutable_data(), ref.shape(0), ref.shape(1));
   }
 
   template <typename T>
   decltype(auto) create_proxy(array_type_t<_aka::Array<T>> & ref,
                               const _aka::Array<T> * /*unused*/) {
+    int size = ref.shape(0);
+    int nb_components = 1;
+    if (ref.ndim() > 1) {
+      nb_components = ref.shape(1);
+    }
     return std::make_unique<_aka::detail::ProxyType_t<_aka::Array<T>>>(
-        ref.mutable_data(), ref.shape(0), ref.shape(1));
+        ref.mutable_data(), size, nb_components);
   }
 
   /* ------------------------------------------------------------------------ */
   template <typename T>
   py::handle aka_array_cast(const _aka::Array<T> & src,
                             py::handle base = handle(), bool writeable = true) {
     array a;
     a = array_type_t<_aka::Array<T>>({src.size(), src.getNbComponent()},
                                      src.storage(), base);
 
     if (not writeable) {
       array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_;
     }
 
     return a.release();
   }
 
   template <typename T>
   py::handle aka_array_cast(const _aka::Vector<T> & src,
                             py::handle base = handle(), bool writeable = true) {
     array a;
     a = array_type_t<_aka::Vector<T>>({src.size()}, src.storage(), base);
 
     if (not writeable) {
       array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_;
     }
 
     return a.release();
   }
 
   template <typename T>
   py::handle aka_array_cast(const _aka::Matrix<T> & src,
                             py::handle base = handle(), bool writeable = true) {
     array a;
     a = array_type_t<_aka::Matrix<T>>({src.size(0), src.size(1)}, src.storage(),
                                       base);
 
     if (not writeable) {
       array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_;
     }
 
     return a.release();
   }
 
   /* ------------------------------------------------------------------------ */
   template <typename AkaArrayType>
   class type_caster<
       AkaArrayType,
       std::enable_if_t<_aka::detail::is_array_type<AkaArrayType>::value>> {
   protected:
     using T = typename AkaArrayType::value_type;
     using type = AkaArrayType;
     using proxy_type = _aka::detail::ProxyType_t<AkaArrayType>;
     using array_type = array_type_t<AkaArrayType>;
 
     std::unique_ptr<proxy_type> array_proxy;
     array_type_t<AkaArrayType> copy_or_ref;
 
   public:
 #if PYBIND11_VERSION_MAJOR >= 2 && PYBIND11_VERSION_MINOR >= 3
     static constexpr auto name = _("AkaArray");
     operator type &&() && { return std::move(*array_proxy); }
     template <typename T_>
     using cast_op_type = pybind11::detail::movable_cast_op_type<T_>;
 #else
     static PYBIND11_DESCR name() { return type_descr(_("AkaArray")); };
     template <typename _T>
     using cast_op_type = pybind11::detail::cast_op_type<_T>;
 #endif
     operator type *() { return array_proxy.get(); }
     operator type &() { return *array_proxy; }
 
     /**
      * Conversion part 1 (Python->C++)
      */
     bool load(handle src, bool convert) {
       bool need_copy = not isinstance<array_type>(src);
 
       auto && fits = [&](auto && aref) {
         auto && dims = aref.ndim();
         if (dims < 1 || dims > 2) {
           return false;
         }
 
         return true;
       };
 
       if (not need_copy) {
         // We don't need a converting copy, but we also need to check whether
         // the strides are compatible with the Ref's stride requirements
         auto aref = py::cast<array_type>(src);
 
         if (not fits(aref)) {
           return false;
         }
         copy_or_ref = std::move(aref);
       } else {
         if (not convert) {
           return false;
         }
 
         auto copy = array_type::ensure(src);
         if (not copy) {
           return false;
         }
 
         if (not fits(copy)) {
           return false;
         }
         copy_or_ref = std::move(array_type::ensure(src));
         loader_life_support::add_patient(copy_or_ref);
       }
 
       AkaArrayType * dispatch = nullptr; // cannot detect T from the expression
       array_proxy = create_proxy(copy_or_ref, dispatch);
       return true;
     }
 
     /**
      * Conversion part 2 (C++ -> Python)
      */
     static handle cast(const type & src, return_value_policy policy,
                        handle parent) {
       switch (policy) {
       case return_value_policy::copy:
         return aka_array_cast<T>(src);
       case return_value_policy::reference_internal:
         return aka_array_cast<T>(src, parent);
       case return_value_policy::reference:
       case return_value_policy::automatic:
       case return_value_policy::automatic_reference:
         return aka_array_cast<T>(src, none());
       default:
         pybind11_fail("Invalid return_value_policy for ArrayProxy type");
       }
     }
   };
 } // namespace detail
 } // namespace pybind11
diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc
index 1a073da25..cf0aff4d9 100644
--- a/python/py_dof_manager.cc
+++ b/python/py_dof_manager.cc
@@ -1,216 +1,228 @@
 /*
  * 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"))
-      .def("zeroResidual", &DOFManager::zeroResidual);
+      .def("zeroResidual", &DOFManager::zeroResidual)
+      .def(
+          "assembleElementalArrayLocalArray",
+          [](DOFManager & self, const Array<Real> & elementary_vect,
+             Array<Real> & array_assembeled, ElementType type,
+             GhostType ghost_type, Real scale_factor) {
+            self.assembleElementalArrayLocalArray(elementary_vect,
+                                                  array_assembeled, type,
+                                                  ghost_type, scale_factor);
+          },
+          py::arg("elementary_vect"), py::arg("array_assembeled"),
+          py::arg("type"), py::arg("ghost_type") = _not_ghost,
+          py::arg("scale_factor") = 1.);
 
   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_fe_engine.cc b/python/py_fe_engine.cc
index 11af53f0a..19f9be213 100644
--- a/python/py_fe_engine.cc
+++ b/python/py_fe_engine.cc
@@ -1,160 +1,266 @@
 /**
  * @file   py_fe_engine.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Nov 27 2019
  * @date last modification: Sat Dec 12 2020
  *
  * @brief  pybind11 interface to FEEngine
  *
  *
  * @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 "py_aka_common.hh"
 /* -------------------------------------------------------------------------- */
 #include <element.hh>
 #include <fe_engine.hh>
 #include <integration_point.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/functional.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 void register_fe_engine(py::module & mod) {
   py::class_<Element>(mod, "Element")
       .def(py::init([](ElementType type, UInt id) {
         return new Element{type, id, _not_ghost};
       }))
       .def(py::init([](ElementType type, UInt id, GhostType ghost_type) {
         return new Element{type, id, ghost_type};
       }))
       .def("__lt__",
            [](Element & self, const Element & other) { return (self < other); })
       .def("__repr__", [](Element & self) { return std::to_string(self); });
 
   mod.attr("ElementNull") = ElementNull;
 
   py::class_<FEEngine>(mod, "FEEngine")
       .def(
           "getNbIntegrationPoints",
           [](FEEngine & fem, const ElementType & type,
              const GhostType & ghost_type) {
             return fem.getNbIntegrationPoints(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def(
           "gradientOnIntegrationPoints",
           [](FEEngine & fem, const Array<Real> & u, Array<Real> & nablauq,
              UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type,
              const Array<UInt> * filter_elements) {
             if (filter_elements == nullptr) {
               // This is due to the ArrayProxy that looses the
               // empty_filter information
               filter_elements = &empty_filter;
             }
             fem.gradientOnIntegrationPoints(u, nablauq, nb_degree_of_freedom,
                                             type, ghost_type, *filter_elements);
           },
           py::arg("u"), py::arg("nablauq"), py::arg("nb_degree_of_freedom"),
           py::arg("type"), py::arg("ghost_type") = _not_ghost,
           py::arg("filter_elements") = nullptr)
+      .def(
+          "gradientOnIntegrationPoints",
+          [](FEEngine & fem, const Array<Real> & u, Array<Real> & nablauq,
+             UInt nb_degree_of_freedom, ElementType type,
+             GhostType ghost_type) {
+            fem.gradientOnIntegrationPoints(u, nablauq, nb_degree_of_freedom,
+                                            type, ghost_type);
+          },
+          py::arg("u"), py::arg("nablauq"), py::arg("nb_degree_of_freedom"),
+          py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def(
           "interpolateOnIntegrationPoints",
-          [](FEEngine & self, const Array<Real> & u, Array<Real> & uq,
+          [](FEEngine & self, Array<Real> & u, Array<Real> & uq,
+             UInt nb_degree_of_freedom, ElementType type,
+             GhostType ghost_type) {
+            self.interpolateOnIntegrationPoints(u, uq, nb_degree_of_freedom,
+                                                type, ghost_type);
+          },
+          py::arg("u"), py::arg("uq"), py::arg("nb_degree_of_freedom"),
+          py::arg("type"), py::arg("ghost_type") = _not_ghost)
+      .def(
+          "interpolateOnIntegrationPoints",
+          [](FEEngine & self, Array<Real> & u, Array<Real> & uq,
              UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type,
              const Array<UInt> * filter_elements) {
             if (filter_elements == nullptr) {
               // This is due to the ArrayProxy that looses the
               // empty_filter information
               filter_elements = &empty_filter;
             }
 
             self.interpolateOnIntegrationPoints(u, uq, nb_degree_of_freedom,
                                                 type, ghost_type,
                                                 *filter_elements);
           },
           py::arg("u"), py::arg("uq"), py::arg("nb_degree_of_freedom"),
           py::arg("type"), py::arg("ghost_type") = _not_ghost,
           py::arg("filter_elements") = nullptr)
       .def(
           "interpolateOnIntegrationPoints",
           [](FEEngine & self, const Array<Real> & u,
              ElementTypeMapArray<Real> & uq,
              const ElementTypeMapArray<UInt> * filter_elements) {
             self.interpolateOnIntegrationPoints(u, uq, filter_elements);
           },
           py::arg("u"), py::arg("uq"), py::arg("filter_elements") = nullptr)
       .def(
           "computeIntegrationPointsCoordinates",
           [](FEEngine & self, ElementTypeMapArray<Real> & coordinates,
              const ElementTypeMapArray<UInt> * filter_elements)
               -> decltype(auto) {
             return self.computeIntegrationPointsCoordinates(coordinates,
                                                             filter_elements);
           },
           py::arg("coordinates"), py::arg("filter_elements") = nullptr)
       .def(
           "assembleFieldLumped",
           [](FEEngine & fem,
              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) {
             fem.assembleFieldLumped(field_funct, matrix_id, dof_id, dof_manager,
                                     type, ghost_type);
           },
           py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"),
           py::arg("dof_manager"), py::arg("type"),
           py::arg("ghost_type") = _not_ghost)
       .def(
           "assembleFieldMatrix",
           [](FEEngine & fem,
              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 = _not_ghost) {
             fem.assembleFieldMatrix(field_funct, matrix_id, dof_id, dof_manager,
                                     type, ghost_type);
           },
           py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"),
           py::arg("dof_manager"), py::arg("type"),
           py::arg("ghost_type") = _not_ghost)
       .def("getElementInradius",
            [](FEEngine & self, const Element & element) {
              return self.getElementInradius(element);
            })
       .def("getNormalsOnIntegrationPoints",
            &FEEngine::getNormalsOnIntegrationPoints, py::arg("type"),
            py::arg("ghost_type") = _not_ghost,
-           py::return_value_policy::reference);
+           py::return_value_policy::reference)
+      .def(
+          "computeNtb",
+          [](FEEngine & fem, const Array<Real> & bs, Array<Real> & Ntbs,
+             ElementType type, GhostType ghost_type) {
+            fem.computeNtb(bs, Ntbs, type, ghost_type);
+          },
+          py::arg("bs"), py::arg("Ntbs"), py::arg("type"),
+          py::arg("ghost_type") = _not_ghost)
+      .def(
+          "computeNtb",
+          [](FEEngine & fem, const Array<Real> & bs, Array<Real> & Ntbs,
+             ElementType type, GhostType ghost_type,
+             const Array<UInt> * filter_elements) {
+            if (filter_elements == nullptr) {
+              // This is due to the ArrayProxy that looses the
+              // empty_filter information
+              filter_elements = &empty_filter;
+            }
+            fem.computeNtb(bs, Ntbs, type, ghost_type, *filter_elements);
+          },
+          py::arg("Ds"), py::arg("BtDs"), py::arg("type"),
+          py::arg("ghost_type") = _not_ghost,
+          py::arg("filter_elements") = nullptr)
+      .def(
+          "computeBtD",
+          [](FEEngine & fem, const Array<Real> & Ds, Array<Real> & BtDs,
+             ElementType type, GhostType ghost_type) {
+            fem.computeBtD(Ds, BtDs, type, ghost_type);
+          },
+          py::arg("Ds"), py::arg("BtDs"), py::arg("type"),
+          py::arg("ghost_type") = _not_ghost)
+      .def(
+          "computeBtD",
+          [](FEEngine & fem, const Array<Real> & Ds, Array<Real> & BtDs,
+             ElementType type, GhostType ghost_type,
+             const Array<UInt> * filter_elements) {
+            if (filter_elements == nullptr) {
+              // This is due to the ArrayProxy that looses the
+              // empty_filter information
+              filter_elements = &empty_filter;
+            }
+            fem.computeBtD(Ds, BtDs, type, ghost_type, *filter_elements);
+          },
+          py::arg("Ds"), py::arg("BtDs"), py::arg("type"),
+          py::arg("ghost_type") = _not_ghost,
+          py::arg("filter_elements") = nullptr)
+      .def(
+          "getShapes",
+          [](FEEngine & fem, const ElementType & type,
+             const GhostType & ghost_type,
+             UInt id) { return fem.getShapes(type, ghost_type, id); },
+          py::arg("type"), py::arg("ghost_type") = _not_ghost,
+          py::arg("id") = 0)
+      .def(
+          "getShapesDerivatives",
+          [](FEEngine & fem, const ElementType & type,
+             const GhostType & ghost_type, UInt id) {
+            return fem.getShapesDerivatives(type, ghost_type, id);
+          },
+          py::arg("type"), py::arg("ghost_type") = _not_ghost,
+          py::arg("id") = 0)
+      .def(
+          "integrate",
+          [](FEEngine & fem, const Array<Real> & f, Array<Real> & intf,
+             UInt nb_degree_of_freedom, ElementType type,
+             GhostType ghost_type) {
+            fem.integrate(f, intf, nb_degree_of_freedom, type, ghost_type);
+          },
+          py::arg("f"), py::arg("intf"), py::arg("nb_degree_of_feedom"),
+          py::arg("type"), py::arg("ghost_type") = _not_ghost)
+      .def(
+          "integrate",
+          [](FEEngine & fem, const Array<Real> & f, Array<Real> & intf,
+             UInt nb_degree_of_freedom, ElementType type, GhostType ghost_type,
+             const Array<UInt> * filter_elements) {
+            if (filter_elements == nullptr) {
+              // This is due to the ArrayProxy that looses the
+              // empty_filter information
+              filter_elements = &empty_filter;
+            }
+            fem.integrate(f, intf, nb_degree_of_freedom, type, ghost_type,
+                          *filter_elements);
+          },
+          py::arg("f"), py::arg("intf"), py::arg("nb_degree_of_feedom"),
+          py::arg("type"), py::arg("ghost_type") = _not_ghost,
+          py::arg("filter_elements") = nullptr);
 
   py::class_<IntegrationPoint>(mod, "IntegrationPoint");
 }
 } // namespace akantu
diff --git a/python/py_heat_transfer_model.cc b/python/py_heat_transfer_model.cc
index fac0a0a6b..9994cbe46 100644
--- a/python/py_heat_transfer_model.cc
+++ b/python/py_heat_transfer_model.cc
@@ -1,102 +1,115 @@
 /**
  * @file   py_heat_transfer_model.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jun 16 2019
  * @date last modification: Sun Jun 16 2019
  *
  * @brief  pybind11 interface to HeatTransferModel
  *
  *
  * @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 <heat_transfer_model.hh>
 #include <non_linear_solver.hh>
 /* -------------------------------------------------------------------------- */
 //#include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 //#include <pybind11/stl.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,                                                              \
       [](HeatTransferModel & self) -> decltype(auto) {                         \
         return self.func_name();                                               \
       },                                                                       \
       py::return_value_policy::reference)
 
 #define def_function(func_name)                                                \
   def(#func_name, [](HeatTransferModel & self) -> decltype(auto) {             \
     return self.func_name();                                                   \
   })
 /* -------------------------------------------------------------------------- */
 
 void register_heat_transfer_model(py::module & mod) {
   py::class_<HeatTransferModelOptions>(mod, "HeatTransferModelOptions")
       .def(py::init<AnalysisMethod>(),
            py::arg("analysis_method") = _explicit_lumped_mass);
 
   py::class_<HeatTransferModel, Model>(mod, "HeatTransferModel",
                                        py::multiple_inheritance())
       .def(py::init<Mesh &, UInt, const ID &>(), py::arg("mesh"),
            py::arg("spatial_dimension") = _all_dimensions,
            py::arg("id") = "heat_transfer_model")
       .def(
           "initFull",
           [](HeatTransferModel & self,
              const HeatTransferModelOptions & options) {
             self.initFull(options);
           },
           py::arg("_analysis_method") = HeatTransferModelOptions())
       .def(
           "initFull",
           [](HeatTransferModel & self,
              const AnalysisMethod & _analysis_method) {
             self.initFull(HeatTransferModelOptions(_analysis_method));
           },
           py::arg("_analysis_method"))
       .def("setTimeStep", &HeatTransferModel::setTimeStep, py::arg("time_step"),
            py::arg("solver_id") = "")
       .def_function(getStableTimeStep)
       .def_function_nocopy(getTemperature)
       .def_function_nocopy(getBlockedDOFs)
+      .def_function_nocopy(getExternalHeatRate)
+      .def_function_nocopy(getBlockedDOFs)
+      .def_function_nocopy(getMesh)
       .def("getTemperatureGradient", &HeatTransferModel::getTemperatureGradient,
            py::arg("el_type"), py::arg("ghost_type") = _not_ghost,
            py::return_value_policy::reference)
       .def("getKgradT", &HeatTransferModel::getKgradT, py::arg("el_type"),
            py::arg("ghost_type") = _not_ghost,
-           py::return_value_policy::reference);
+           py::return_value_policy::reference)
+      .def("applyBC",
+           [](HeatTransferModel & self, BC::Dirichlet::DirichletFunctor & func,
+              const std::string & element_group) {
+             self.applyBC(func, element_group);
+           })
+      .def("applyBC",
+           [](HeatTransferModel & self, BC::Neumann::NeumannFunctor & func,
+              const std::string & element_group) {
+             self.applyBC(func, element_group);
+           });
 }
 
 } // namespace akantu
diff --git a/python/py_mesh.cc b/python/py_mesh.cc
index 3c4de188a..b658e1dae 100644
--- a/python/py_mesh.cc
+++ b/python/py_mesh.cc
@@ -1,254 +1,271 @@
 /**
  * @file   py_mesh.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Philip Mueller <philip.paul.mueller@bluemail.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jun 16 2019
  * @date last modification: Mon Mar 15 2021
  *
  * @brief  pybind11 interface to Mesh
  *
  *
  * @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_config.hh"
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <mesh.hh>
 #include <mesh_accessor.hh>
 #include <mesh_utils.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 namespace {
   /* ------------------------------------------------------------------------ */
   template <typename T>
   decltype(auto) register_element_type_map_array(py::module & mod,
                                                  const std::string & name) {
     return py::class_<ElementTypeMapArray<T>,
                       std::shared_ptr<ElementTypeMapArray<T>>>(
                mod, ("ElementTypeMapArray" + name).c_str())
         .def(py::init<const ID &, const ID &>(),
              py::arg("id") = "by_element_type_array",
              py::arg("parent_id") = "no_parent")
         .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, py::keep_alive<0, 1>())
         .def(
             "elementTypes",
             [](ElementTypeMapArray<T> & self, UInt _dim, GhostType _ghost_type,
                ElementKind _kind) -> std::vector<ElementType> {
               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)
         .def(
             "initialize",
             [](ElementTypeMapArray<T> & self, const Mesh & mesh,
                GhostType ghost_type = _casper, UInt nb_component = 1,
                UInt spatial_dimension = UInt(-2),
                ElementKind element_kind = _ek_not_defined,
                bool with_nb_element = false,
                bool with_nb_nodes_per_element = false, T default_value = T(),
                bool do_not_default = false) {
               self.initialize(
                   mesh, _ghost_type = ghost_type, _nb_component = nb_component,
                   _spatial_dimension = (spatial_dimension == UInt(-2)
                                             ? mesh.getSpatialDimension()
                                             : spatial_dimension),
                   _element_kind = element_kind,
                   _with_nb_element = with_nb_element,
                   _with_nb_nodes_per_element = with_nb_nodes_per_element,
                   _default_value = default_value,
                   _do_not_default = do_not_default);
             },
             py::arg("mesh"), py::arg("ghost_type") = _casper,
             py::arg("nb_component") = 1,
             py::arg("spatial_dimension") = UInt(-2),
             py::arg("element_kind") = _ek_not_defined,
             py::arg("with_nb_element") = false,
             py::arg("with_nb_nodes_per_element") = false,
             py::arg("default_value") = T(), py::arg("do_not_default") = false);
   }
 } // namespace
 /* -------------------------------------------------------------------------- */
 void register_mesh(py::module & mod) {
   py::class_<Mesh::PeriodicSlaves>(mod, "PeriodicSlaves")
       .def(
           "__iter__",
           [](Mesh::PeriodicSlaves & _this) {
             py::make_iterator(_this.begin(), _this.end());
           },
           py::keep_alive<0, 1>());
 
   py::class_<MeshData>(mod, "MeshData")
       .def(
           "getElementalDataUInt",
           [](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 &>(), py::arg("spatial_dimension"),
            py::arg("id") = "mesh")
       .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(
           "getConnectivities",
           [](Mesh & self) -> decltype(auto) {
             return self.getConnectivities();
           },
           py::return_value_policy::reference)
       .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("isDistributed",
            [](const Mesh & self) { return self.isDistributed(); })
       .def("fillNodesToElements", &Mesh::fillNodesToElements,
            py::arg("dimension") = _all_dimensions)
       .def("getAssociatedElements",
            [](Mesh & self, const UInt & node, py::list list) {
              Array<Element> elements;
              self.getAssociatedElements(node, elements);
              for (auto && element : elements) {
                list.append(element);
              }
            })
       .def("makePeriodic",
            [](Mesh & self, const SpatialDirection & direction) {
              self.makePeriodic(direction);
            })
       .def(
           "getNbElement",
           [](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) {
             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(
+          "getNbNodesPerElement",
+          [](ElementType & type) { return Mesh::getNbNodesPerElement(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)
       .def(
           "hasDataReal",
           [](Mesh & _this, const ID & name, ElementType type,
              GhostType ghost_type) -> bool {
             return _this.hasData<Real>(name, type, ghost_type);
           },
           py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def("isPeriodic", [](const Mesh & _this) { return _this.isPeriodic(); })
       .def("getPeriodicMaster", &Mesh::getPeriodicMaster)
       .def("getPeriodicSlaves", &Mesh::getPeriodicSlaves)
       .def("isPeriodicSlave", &Mesh::isPeriodicSlave)
       .def("isPeriodicMaster", &Mesh::isPeriodicMaster)
       .def(
           "getMeshFacets",
           [](const Mesh & self) -> const Mesh & {
             return self.getMeshFacets();
           },
           py::return_value_policy::reference)
       .def("initMeshFacets", &Mesh::initMeshFacets,
-           py::arg("id") = "mesh_facets", py::return_value_policy::reference);
+           py::arg("id") = "mesh_facets", py::return_value_policy::reference)
+      .def(
+          "elementTypes",
+          [](Mesh & self, UInt spatial_dimension, GhostType ghost_type,
+             ElementKind kind) -> std::vector<ElementType> {
+            auto types = self.elementTypes(spatial_dimension, ghost_type, kind);
+            std::vector<ElementType> _types;
+            for (auto && t : types) {
+              _types.push_back(t);
+            }
+            return _types;
+          },
+          py::arg("spatial_dimension") = _all_dimensions,
+          py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular);
 
-  /* ------------------------------------------------------------------------ */
+  /* ------------------------------------------------------------------------
+   */
   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(
           "resizeNodes",
           [](MeshAccessor & self, UInt new_size) -> void {
             self.resizeNodes(new_size);
           },
           py::arg("new_size"))
       .def("makeReady", &MeshAccessor::makeReady);
 
   register_element_type_map_array<Real>(mod, "Real");
   register_element_type_map_array<UInt>(mod, "UInt");
   register_element_type_map_array<bool>(mod, "bool");
   // register_element_type_map_array<std::string>(mod, "String");
 }
 } // namespace akantu
diff --git a/python/py_model.cc b/python/py_model.cc
index 4f0a9312c..96b091cb9 100644
--- a/python/py_model.cc
+++ b/python/py_model.cc
@@ -1,149 +1,150 @@
 /**
  * @file   py_model.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Emil Gallyamov <emil.gallyamov@epfl.ch>
  * @author Philip Mueller <philip.paul.mueller@bluemail.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jun 16 2019
  * @date last modification: Sat Mar 13 2021
  *
  * @brief  pybind11 interface to Model and parent classes
  *
  *
  * @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 <model.hh>
 #include <non_linear_solver.hh>
 #include <solver_callback.hh>
 #include <sparse_matrix_aij.hh>
 #include <time_step_solver.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_<ModelSolver, SolverCallback, Parsable>(mod, "ModelSolver",
                                                     py::multiple_inheritance())
       .def(
           "getNonLinearSolver",
           [](ModelSolver & self, const ID & solver_id) -> NonLinearSolver & {
             return self.getNonLinearSolver(solver_id);
           },
           py::arg("solver_id") = "", py::return_value_policy::reference)
       .def(
           "getTimeStepSolver",
           [](ModelSolver & self, const ID & solver_id) -> TimeStepSolver & {
             return self.getTimeStepSolver(solver_id);
           },
           py::arg("solver_id") = "", py::return_value_policy::reference)
       .def(
           "solveStep",
           [](ModelSolver & self, const ID & solver_id) {
             self.solveStep(solver_id);
           },
           py::arg("solver_id") = "")
       .def(
           "solveStep",
           [](ModelSolver & self, SolverCallback & callback,
              const ID & solver_id) { self.solveStep(callback, solver_id); },
           py::arg("callback"), py::arg("solver_id") = "");
 
   py::class_<Model, ModelSolver>(mod, "Model", py::multiple_inheritance())
+      .def("getSpatialDimension", &Model::getSpatialDimension)
       .def("setBaseName", &Model::setBaseName)
       .def("setDirectory", &Model::setDirectory)
       .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", [](Model & self) { self.dump(); })
       .def(
           "dump", [](Model & self, UInt step) { self.dump(step); },
           py::arg("step"))
       .def(
           "dump",
           [](Model & self, Real time, UInt step) { self.dump(time, step); },
           py::arg("time"), py::arg("step"))
       .def(
           "dump",
           [](Model & self, const std::string & dumper) { self.dump(dumper); },
           py::arg("dumper_name"))
       .def(
           "dump",
           [](Model & self, const std::string & dumper, UInt step) {
             self.dump(dumper, step);
           },
           py::arg("dumper_name"), py::arg("step"))
       .def(
           "dump",
           [](Model & self, const std::string & dumper, Real time, UInt step) {
             self.dump(dumper, time, step);
           },
           py::arg("dumper_name"), py::arg("time"), py::arg("step"))
       .def("initNewSolver", &Model::initNewSolver)
       .def(
           "getNewSolver",
           [](Model & self, const std::string id,
              const TimeStepSolverType & time,
              const NonLinearSolverType & type) {
             self.getNewSolver(id, time, type);
           },
           py::return_value_policy::reference)
       .def(
           "setIntegrationScheme",
           [](Model & self, const std::string id, const std::string primal,
              const IntegrationSchemeType & scheme_type,
              IntegrationScheme::SolutionType solution_type) {
             self.setIntegrationScheme(id, primal, scheme_type, solution_type);
           },
           py::arg("id"), py::arg("primal"), py::arg("scheme_type"),
           py::arg("solution_type") =
               IntegrationScheme::SolutionType::_not_defined)
       // .def("setIntegrationScheme",
       //      [](Model & self, const std::string id, const std::string primal,
       //         std::unique_ptr<IntegrationScheme> & scheme,
       //         IntegrationScheme::SolutionType solution_type) {
       //        self.setIntegrationScheme(id, primal, scheme, solution_type);
       //      })
 
       .def("getDOFManager", &Model::getDOFManager,
            py::return_value_policy::reference)
       .def("assembleMatrix", &Model::assembleMatrix);
 }
 
 } // namespace akantu
diff --git a/python/py_solid_mechanics_model.cc b/python/py_solid_mechanics_model.cc
index bd8d801fc..ad8309338 100644
--- a/python/py_solid_mechanics_model.cc
+++ b/python/py_solid_mechanics_model.cc
@@ -1,186 +1,185 @@
 /**
  * @file   py_solid_mechanics_model.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jun 16 2019
  * @date last modification: Sat Mar 13 2021
  *
  * @brief  pybind11 interface to SolidMechanicsModel
  *
  *
  * @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 <non_linear_solver.hh>
 #include <solid_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,                                                              \
       [](SolidMechanicsModel & self) -> decltype(auto) {                       \
         return self.func_name();                                               \
       },                                                                       \
       py::return_value_policy::reference)
 
 #define def_function(func_name)                                                \
   def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) {           \
     return self.func_name();                                                   \
   })
 
 /* -------------------------------------------------------------------------- */
 void register_solid_mechanics_model(py::module & mod) {
 
   py::class_<SolidMechanicsModelOptions>(mod, "SolidMechanicsModelOptions")
       .def(py::init<AnalysisMethod>(),
            py::arg("_analysis_method") = _explicit_lumped_mass);
 
   py::class_<SolidMechanicsModel, Model>(mod, "SolidMechanicsModel",
                                          py::multiple_inheritance())
       .def(py::init<Mesh &, UInt, const ID &, std::shared_ptr<DOFManager>,
                     const ModelType>(),
            py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions,
            py::arg("id") = "solid_mechanics_model",
            py::arg("dof_manager") = nullptr,
            py::arg("model_type") = ModelType::_solid_mechanics_model)
       .def(
           "initFull",
           [](SolidMechanicsModel & self,
              const SolidMechanicsModelOptions & options) {
             self.initFull(options);
           },
           py::arg("option") = SolidMechanicsModelOptions())
       .def(
           "initFull",
           [](SolidMechanicsModel & self,
              const AnalysisMethod & analysis_method) {
             self.initFull(_analysis_method = analysis_method);
           },
           py::arg("_analysis_method"))
       .def_deprecated("applyDirichletBC", "Deprecated: use applyBC")
       .def("applyBC",
            [](SolidMechanicsModel & self,
               BC::Dirichlet::DirichletFunctor & func,
               const std::string & element_group) {
              self.applyBC(func, element_group);
            })
       .def("applyBC",
            [](SolidMechanicsModel & self, BC::Neumann::NeumannFunctor & func,
               const std::string & element_group) {
              self.applyBC(func, element_group);
            })
       .def("setTimeStep", &SolidMechanicsModel::setTimeStep,
            py::arg("time_step"), py::arg("solver_id") = "")
       .def(
           "getEnergy",
           [](SolidMechanicsModel & self, const std::string & energy_id) {
             return self.getEnergy(energy_id);
           },
           py::arg("energy_id"))
       .def(
           "getEnergy",
           [](SolidMechanicsModel & self, const std::string & energy_id,
              const std::string & group_id) {
             return self.getEnergy(energy_id, group_id);
           },
           py::arg("energy_id"), py::arg("group_id"))
 
       .def_function(assembleStiffnessMatrix)
       .def_function(assembleInternalForces)
       .def_function(assembleMass)
       .def_function(assembleMassLumped)
       .def_function(getStableTimeStep)
       .def_function_nocopy(getExternalForce)
       .def_function_nocopy(getDisplacement)
       .def_function_nocopy(getPreviousDisplacement)
       .def_function_nocopy(getCurrentPosition)
       .def_function_nocopy(getIncrement)
-      .def_function_nocopy(getInternalForce)
       .def_function_nocopy(getMass)
       .def_function_nocopy(getVelocity)
       .def_function_nocopy(getAcceleration)
       .def_function_nocopy(getInternalForce)
       .def_function_nocopy(getBlockedDOFs)
       .def_function_nocopy(getMesh)
       .def(
           "getMaterial",
           [](SolidMechanicsModel & self, UInt material_id) -> decltype(auto) {
             return self.getMaterial(material_id);
           },
           py::arg("material_id"), py::return_value_policy::reference)
       .def(
           "getMaterial",
           [](SolidMechanicsModel & self, const ID & material_name)
               -> decltype(auto) { return self.getMaterial(material_name); },
           py::arg("material_name"), py::return_value_policy::reference)
       .def(
           "getMaterial",
           [](SolidMechanicsModel & self, const Element & element)
               -> decltype(auto) { return self.getMaterial(element); },
           py::arg("element"), py::return_value_policy::reference)
 
       .def("getNbMaterials", &SolidMechanicsModel::getNbMaterials)
       .def("getMaterialIndex", &SolidMechanicsModel::getMaterialIndex)
       .def("setMaterialSelector",
            [](SolidMechanicsModel & self,
               std::shared_ptr<MaterialSelector> material_selector) {
              std::cout << (*material_selector)(ElementNull) << std::endl;
              self.setMaterialSelector(material_selector);
            })
       .def("getMaterialSelector", &SolidMechanicsModel::getMaterialSelector)
       .def(
           "getMaterialByElement",
           [](const SolidMechanicsModel & self) -> decltype(auto) {
             return self.getMaterialByElement();
           },
           py::return_value_policy::reference, py::keep_alive<0, 1>())
       .def("reassignMaterial", &SolidMechanicsModel::reassignMaterial)
       .def(
           "registerNewMaterial",
           [](SolidMechanicsModel & self, const ID & mat_name,
              const ID & mat_type, const ID & opt_param) -> decltype(auto) {
             return self.registerNewMaterial(mat_name, mat_type, opt_param);
           },
           py::arg("material_name"), py::arg("material_type"),
           py::arg("option") = "", py::return_value_policy::reference)
       .def("initMaterials", &SolidMechanicsModel::initMaterials)
       .def("flattenInternal", &SolidMechanicsModel::flattenInternal,
            py::return_value_policy::reference)
       .def("inflateInternal", &SolidMechanicsModel::inflateInternal,
            py::return_value_policy::reference);
 }
 
 } // namespace akantu