diff --git a/examples/python/plate-hole/material.dat b/examples/python/plate-hole/material.dat
index 8f20aa1e1..bf0627a5a 100644
--- a/examples/python/plate-hole/material.dat
+++ b/examples/python/plate-hole/material.dat
@@ -1,6 +1,6 @@
 material elastic [
-         name = steel
-	 rho = 7800     # density
-	 E   = 2.1e11   # young's modulus
-	 nu  = 0.3      # poisson's ratio
+   name = steel
+   rho = 7800     # density
+   E   = 2.1e11   # young's modulus
+   nu  = 0.3      # poisson's ratio
 ]
diff --git a/python/py_aka_array.hh b/python/py_aka_array.hh
index a2fd3956a..c3a594dd3 100644
--- a/python/py_aka_array.hh
+++ b/python/py_aka_array.hh
@@ -1,275 +1,224 @@
 /**
  * Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * This file is part of Akantu
  *
  * 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/eigen.h>
 #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(Int /*size*/, Int /*nb_component*/)
-    final {}
+    void allocate(Int /*size*/, Int /*nb_component*/) final {}
 
     // allocate and initialize the memory
     void allocate(Int /*size*/, Int /*nb_component*/,
                   const T & /*value*/) final {}
 
   public:
     ArrayProxy(T * data, Int size, Int nb_component) {
       this->values = data;
       this->size_ = size;
       this->nb_component = nb_component;
     }
 
     ArrayProxy(const Array<T> & src) {
       this->values = src.data();
       this->size_ = src.size();
       this->nb_component = src.getNbComponent();
     }
 
     ~ArrayProxy() override { this->values = nullptr; }
 
     void resize(Int 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(Int new_size) final {
       if (new_size != this->size()) {
         AKANTU_EXCEPTION("cannot resize a temporary array");
       }
     }
 
-    void reserve(Int /*size*/, Int /*new_size*/)
-    final { AKANTU_EXCEPTION("cannot resize a temporary array"); }
+    void reserve(Int /*size*/, Int /*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*/) {
     return std::make_unique<_aka::detail::ProxyType_t<_aka::Array<T>>>(
         ref.mutable_data(), ref.shape(0), ref.shape(1));
   }
 
   /* ------------------------------------------------------------------------ */
   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.data(), 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.data(), 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.data(),
-  //                                     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_akantu.cc b/python/py_akantu.cc
index 4471660fe..08613e3c2 100644
--- a/python/py_akantu.cc
+++ b/python/py_akantu.cc
@@ -1,159 +1,166 @@
 /**
  * Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * This file is part of Akantu
  *
  * 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"
+// for NLSNotConvergedException
+#include "non_linear_solver.hh"
 /* -------------------------------------------------------------------------- */
 #include "py_aka_common.hh"
 #include "py_aka_error.hh"
 #include "py_boundary_conditions.hh"
 #include "py_dof_manager.hh"
 #include "py_dumpable.hh"
 #include "py_fe_engine.hh"
 #include "py_group_manager.hh"
 #include "py_integration_scheme.hh"
 #include "py_mesh.hh"
 #include "py_model.hh"
 #include "py_parser.hh"
 #include "py_solver.hh"
 
 #if defined(AKANTU_SOLID_MECHANICS)
 #include "py_material.hh"
 #include "py_material_selector.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_fragment_manager.hh"
 #include "py_solid_mechanics_model_cohesive.hh"
 #endif
 
 #if defined(AKANTU_CONTACT_MECHANICS)
 #include "py_contact_mechanics_model.hh"
 #include "py_model_couplers.hh"
 #endif
 
 #if defined(AKANTU_PHASE_FIELD)
 #include "py_phase_field_model.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);
   register_dumpable(mod);
   register_mesh(mod);
 
   register_fe_engine(mod);
 
   register_integration_schemes(mod);
   register_dof_manager(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);
   register_material_selector(mod);
 #endif
 
 #if defined(AKANTU_COHESIVE_ELEMENT)
   register_solid_mechanics_model_cohesive(mod);
   register_fragment_manager(mod);
 #endif
 
 #if defined(AKANTU_STRUCTURAL_MECHANICS)
   register_structural_mechanics_model(mod);
 #endif
 
 #if defined(AKANTU_CONTACT_MECHANICS)
   register_contact_mechanics_model(mod);
   register_model_couplers(mod);
 #endif
 
 #if defined(AKANTU_PHASE_FIELD)
   register_phase_field_model(mod);
   register_phase_field_coupler(mod);
 #endif
 }
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 PYBIND11_MODULE(py11_akantu, mod) {
   mod.doc() = "Akantu python interface";
 
   static py::exception<akantu::debug::Exception> akantu_exception(mod,
                                                                   "Exception");
 
+  static py::exception<akantu::debug::NLSNotConvergedException>
+      akantu_exception_nls_not_converged(mod, "NLSNotConvergedException");
+
   py::register_exception_translator([](std::exception_ptr ptr) {
     try {
       if (ptr) {
         std::rethrow_exception(ptr);
       }
+    } catch (akantu::debug::NLSNotConvergedException & e) {
+      akantu_exception_nls_not_converged(e.info().c_str());
     } 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;
+            return false;
 #endif
           })
       .def("getVersion", &akantu::getVersion);
 
 } // Module akantu
diff --git a/python/py_fe_engine.cc b/python/py_fe_engine.cc
index eaa70207c..b0c6f03d0 100644
--- a/python/py_fe_engine.cc
+++ b/python/py_fe_engine.cc
@@ -1,135 +1,138 @@
 /**
  * Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * This file is part of Akantu
  *
  * 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>
 /* -------------------------------------------------------------------------- */
+#include <memory>
+/* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 void register_fe_engine(py::module & mod) {
   py::class_<Element>(mod, "Element")
-      .def(py::init([](ElementType type, Int id) {
-        return new Element{type, id, _not_ghost};
-      }))
       .def(py::init([](ElementType type, Int id, GhostType ghost_type) {
-        return new Element{type, id, ghost_type};
-      }))
+             return std::make_unique<Element>(Element{type, id, ghost_type});
+           }),
+           py::arg("type"), py::arg("ghost_type"),
+           py::arg("ghost_type") = _not_ghost)
       .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, ElementType type, GhostType ghost_type) {
             return fem.getNbIntegrationPoints(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost)
+      .def("initShapeFunctions", &FEEngine::initShapeFunctions,
+           py::arg("ghost_type") = _not_ghost)
       .def(
           "gradientOnIntegrationPoints",
           [](FEEngine & fem, const Array<Real> & u, Array<Real> & nablauq,
              const Int nb_degree_of_freedom, ElementType type,
              GhostType ghost_type, const Array<Idx> & filter_elements) {
             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") = empty_filter)
       .def(
           "interpolateOnIntegrationPoints",
           [](FEEngine & self, const Array<Real> & u, Array<Real> & uq,
              Int nb_degree_of_freedom, ElementType type, GhostType ghost_type,
              const Array<Idx> & filter_elements) {
             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") = empty_filter)
       .def(
           "interpolateOnIntegrationPoints",
           [](FEEngine & self, const Array<Real> & u,
              ElementTypeMapArray<Real> & uq,
              const ElementTypeMapArray<Idx> * 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<Idx> * 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::class_<IntegrationPoint>(mod, "IntegrationPoint");
 }
 } // namespace akantu
diff --git a/python/py_mesh.cc b/python/py_mesh.cc
index d25b203cd..7bde07c20 100644
--- a/python/py_mesh.cc
+++ b/python/py_mesh.cc
@@ -1,245 +1,247 @@
 /**
  * Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * This file is part of Akantu
  *
  * 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>
   auto register_element_type_map_array(py::module & mod,
                                        const std::string & name) {
     auto element_type_map_class =
         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, Int nb_component,
                    Int spatial_dimension, ElementKind element_kind,
                    bool with_nb_element, bool with_nb_nodes_per_element,
                    T default_value, bool do_not_default) {
                   self.initialize(
                       mesh, _ghost_type = ghost_type,
                       _nb_component = nb_component,
                       _spatial_dimension =
                           (spatial_dimension == -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") = -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);
 
     return element_type_map_class;
   }
 } // namespace
 /* -------------------------------------------------------------------------- */
 void register_mesh(py::module & mod) {
   py::class_<Mesh::PeriodicSlaves>(mod, "PeriodicSlaves")
       .def(
           "__iter__",
           [](Mesh::PeriodicSlaves & _this) {
             return 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<Int, 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 Idx & 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 Int 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(
           "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("getUpperBounds", &Mesh::getUpperBounds)
+      .def("getLowerBounds", &Mesh::getLowerBounds)
       .def("initMeshFacets", &Mesh::initMeshFacets,
            py::arg("id") = "mesh_facets", 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, Int 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, Int 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<Int>(mod, "Int");
   register_element_type_map_array<bool>(mod, "bool");
   // register_element_type_map_array<std::string>(mod, "String");
 }
 } // namespace akantu