diff --git a/python/pybind11/py_aka_array.cc b/python/pybind11/py_aka_array.cc
index fcf089a89..3c5ee7a13 100644
--- a/python/pybind11/py_aka_array.cc
+++ b/python/pybind11/py_aka_array.cc
@@ -1,163 +1,215 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 
 namespace py = pybind11;
 namespace _aka = akantu;
 
 namespace akantu {
 template <typename T> class ArrayProxy : public Array<T> {
 protected:
   // deallocate the memory
   void deallocate() override final {}
 
   // allocate the memory
-  void allocate(UInt size, UInt nb_component) override final {}
+  void allocate(__attribute__((unused)) UInt size,
+                __attribute__((unused)) UInt nb_component) override final {}
 
   // allocate and initialize the memory
-  void allocate(UInt size, UInt nb_component, const T & value) override final {}
+  void allocate(__attribute__((unused)) UInt size,
+                __attribute__((unused)) UInt nb_component,
+                __attribute__((unused)) const T & value) override 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();
   }
 
-  void resize(UInt size, const T & val) override final {
+  void resize(__attribute__((unused)) UInt size,
+              __attribute__((unused)) const T & val) override final {
     AKANTU_EXCEPTION("cannot resize a temporary array");
   }
 
-  void resize(UInt new_size) override final {
+  void resize(__attribute__((unused)) UInt new_size) override final {
     AKANTU_EXCEPTION("cannot resize a temporary array");
   }
 
-  void reserve(UInt new_size) override final {
+  void reserve(__attribute__((unused)) UInt new_size) override final {
     AKANTU_EXCEPTION("cannot resize a temporary array");
   }
 };
 
 } // namespace akantu
 
 namespace pybind11 {
 namespace detail {
   /* ------------------------------------------------------------------------ */
   template <typename T>
   py::handle aka_array_cast(const _aka::Array<T> & src,
                             py::handle base = handle(), bool writeable = true) {
     constexpr ssize_t elem_size = sizeof(T);
     array a;
     a = array({src.size(), src.getNbComponent()},
               {elem_size * src.getNbComponent(), elem_size}, src.storage(),
               base);
 
     if (not writeable)
       array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_;
 
     return a.release();
   }
 
   template <typename T> struct type_caster<_aka::Array<T>> {
   protected:
     using py_array = array_t<T, array::forcecast | array::c_style>;
+    using type = _aka::Array<T>;
 
   public:
-    static constexpr auto name =
-        _("numpy.ndarray[") + npy_format_descriptor<T>::name +
-        _(", flags.writeable") + _(", flags.c_contiguous") + _("]");
+    PYBIND11_TYPE_CASTER(type, _("Array<T>"));
+
+    // static constexpr auto name =
+    //     _("numpy.ndarray[") + npy_format_descriptor<T>::name +
+    //     _(", flags.writeable") + _(", flags.c_contiguous") + _("]");
 
     /**
      * Conversion part 1 (Python->C++)
      */
     bool load(handle src, bool convert) {
       bool need_copy = not isinstance<py_array>(src);
 
       auto && fits = [&](auto && a) {
         auto && dims = copy_or_ref.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
         py_array aref = reinterpret_borrow<py_array>(src);
 
         if (not fits(aref)) {
           return false;
         }
 
         copy_or_ref = std::move(aref);
       } else {
         if (not convert) {
           return false;
         }
 
         auto copy = py_array::ensure(src);
         if (not copy) {
           return false;
         }
 
         if (not fits(copy)) {
           return false;
         }
         copy_or_ref = std::move(py_array::ensure(src));
         loader_life_support::add_patient(copy_or_ref);
       }
 
       array_proxy.reset(new _aka::ArrayProxy<T>> (copy_or_ref.mutable_data(),
                                                   copy_or_ref.shape(0),
                                                   copy_or_ref.shape(1)));
       return true;
     }
 
-    operator _aka::Array<T> *() { return array_proxy.get(); }
-    operator _aka::Array<T> &() { return *array_proxy; }
-    template <typename _T>
-    using cast_op_type = pybind11::detail::cast_op_type<_T>;
+    // operator _aka::Array<T> *() { return array_proxy.get(); }
+    // operator _aka::Array<T> &() { return *array_proxy; }
+    // template <typename _T>
+    // using cast_op_type = pybind11::detail::cast_op_type<_T>;
 
     /**
      * Conversion part 2 (C++ -> Python)
      */
     static handle cast(const _aka::Array<T> & 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");
       }
     }
 
   protected:
     std::unique_ptr<_aka::Array<T>> array_proxy;
     py_array copy_or_ref;
   };
+
+  /**
+   * Type caster for Vector classes
+   */
+  template <template <typename> class V, typename T> struct type_caster<V<T>> {
+    using type = V<T>;
+    template <typename U>
+    using array_type = array_t<U, array::c_style | array::forcecast>;
+
+  public:
+    PYBIND11_TYPE_CASTER(type, _("Vector<T>"));
+
+    /**
+     * Conversion part 1 (Python->C++): convert a PyObject into a Vector
+     * instance or return false upon failure. The second argument
+     * indicates whether implicit conversions should be applied.
+     */
+    bool load(handle src, bool convert) {
+      if (!convert && !array_type<typename type::value_type>::check_(src))
+        return false;
+
+      auto buf = array_type<typename type::value_type>::ensure(src);
+      value.move(_aka::VectorProxy<T>(buf));
+
+      return true;
+    }
+
+    /**
+     * Conversion part 2 (C++ -> Python): convert a grid instance into
+     * a Python object. The second and third arguments are used to
+     * indicate the return value policy and parent object (for
+     * ``return_value_policy::reference_internal``) and are generally
+     * ignored by implicit casters.
+     *
+     * TODO: do not ignore policy (see pybind11/eigen.h)
+     */
+    static handle cast(const type & src, return_value_policy /* policy */,
+                       handle /*parent*/) {
+      // none() passed as parent to get a correct nocopy
+      auto a = array_type<typename type::value_type>(src.size(), src.storage(),
+                                                     none());
+      return a.release();
+    }
+  };
+
 } // namespace detail
 } // namespace pybind11
 
-namespace{
-  py::module & register_arrays(py::module & mod) {
+namespace {
+py::module & register_arrays(py::module & mod) {
 
-    py::class_<_aka::Vector<bool>>(mod, "VectorBool");
-    py::class_<_aka::Vector<double>>(mod, "VectorReal");    
-    return mod;
-  }
+  // py::class_<_aka::Vector<bool>>(mod, "VectorBool");
+  // py::class_<_aka::Vector<double>>(mod, "VectorReal");
+  return mod;
 }
+} // namespace
diff --git a/python/pybind11/py_aka_boundary_conditions.cc b/python/pybind11/py_aka_boundary_conditions.cc
index 1001666b8..119f6bf63 100644
--- a/python/pybind11/py_aka_boundary_conditions.cc
+++ b/python/pybind11/py_aka_boundary_conditions.cc
@@ -1,45 +1,66 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 #include "boundary_condition_python_functor.hh"
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 namespace _aka = akantu;
 
 namespace {
 
 /* -------------------------------------------------------------------------- */
 
 class PyDirichletFunctor : public _aka::BC::DirichletFunctor {
 public:
   /* Inherit the constructors */
   using _aka::BC::DirichletFunctor::DirichletFunctor;
 
   /* Trampoline (need one for each virtual function) */
   void operator()(_aka::UInt node, _aka::Vector<bool> & flags,
                   _aka::Vector<_aka::Real> & primal,
                   const _aka::Vector<_aka::Real> & coord) const override {
 
     PYBIND11_OVERLOAD_NAME(void, _aka::BC::DirichletFunctor,
                            "__call__", operator(), node, flags, primal, coord);
   }
 };
 /* -------------------------------------------------------------------------- */
 
+class PyNeumannFunctor : public _aka::BC::NeumannFunctor {
+public:
+  /* Inherit the constructors */
+  using _aka::BC::NeumannFunctor::NeumannFunctor;
+
+  /* Trampoline (need one for each virtual function) */
+  void operator()(const _aka::IntegrationPoint & quad_point,
+                  _aka::Vector<_aka::Real> & dual,
+                  const _aka::Vector<_aka::Real> & coord,
+                  const _aka::Vector<_aka::Real> & normals) const override {
+
+    PYBIND11_OVERLOAD_PURE_NAME(void, _aka::BC::NeumannFunctor,
+                           "__call__", operator(), quad_point, dual, coord,
+                           normals);
+  }
+};
+/* -------------------------------------------------------------------------- */
+
 py::module & register_boundary_conditions(py::module & mod) {
 
   py::class_<_aka::BC::Functor>(mod, "BCFunctor");
   py::class_<_aka::BC::DirichletFunctor, PyDirichletFunctor, _aka::BC::Functor>(
       mod, "DirichletFunctor")
       .def(py::init())
       .def(py::init<_aka::SpatialDirection>());
 
+  py::class_<_aka::BC::NeumannFunctor, PyNeumannFunctor, _aka::BC::Functor>(
+      mod, "NeumannFunctor")
+      .def(py::init());
 
   return mod;
 } // namespace
 
 } // namespace
diff --git a/python/pybind11/py_aka_common.cc b/python/pybind11/py_aka_common.cc
index 458d205ac..0729b7505 100644
--- a/python/pybind11/py_aka_common.cc
+++ b/python/pybind11/py_aka_common.cc
@@ -1,115 +1,118 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
+#include "integration_point.hh"
 #include "mesh.hh"
 #include "parser.hh"
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 namespace _aka = akantu;
 
 namespace {
 
 /* -------------------------------------------------------------------------- */
 
 py::module & register_enums(py::module & mod) {
   py::enum_<_aka::SpatialDirection>(mod, "SpatialDirection")
       .value("_x", _aka::_x)
       .value("_y", _aka::_y)
       .value("_z", _aka::_z)
       .export_values();
 
   py::enum_<_aka::AnalysisMethod>(mod, "AnalysisMethod")
       .value("_static", _aka::_static)
       .value("_implicit_dynamic", _aka::_implicit_dynamic)
       .value("_explicit_lumped_mass", _aka::_explicit_lumped_mass)
       .value("_explicit_lumped_capacity", _aka::_explicit_lumped_capacity)
       .value("_explicit_consistent_mass", _aka::_explicit_consistent_mass)
       .export_values();
 
   py::enum_<_aka::NonLinearSolverType>(mod, "NonLinearSolverType")
       .value("_nls_linear", _aka::_nls_linear)
       .value("_nls_newton_raphson", _aka::_nls_newton_raphson)
       .value("_nls_newton_raphson_modified", _aka::_nls_newton_raphson_modified)
       .value("_nls_lumped", _aka::_nls_lumped)
       .value("_nls_auto", _aka::_nls_auto)
       .export_values();
 
   py::enum_<_aka::TimeStepSolverType>(mod, "TimeStepSolverType")
       .value("_tsst_static", _aka::_tsst_static)
       .value("_tsst_dynamic", _aka::_tsst_dynamic)
       .value("_tsst_dynamic_lumped", _aka::_tsst_dynamic_lumped)
       .value("_tsst_not_defined", _aka::_tsst_not_defined)
       .export_values();
 
   py::enum_<_aka::IntegrationSchemeType>(mod, "IntegrationSchemeType")
       .value("_ist_pseudo_time", _aka::_ist_pseudo_time)
       .value("_ist_forward_euler", _aka::_ist_forward_euler)
       .value("_ist_trapezoidal_rule_1", _aka::_ist_trapezoidal_rule_1)
       .value("_ist_backward_euler", _aka::_ist_backward_euler)
       .value("_ist_central_difference", _aka::_ist_central_difference)
       .value("_ist_fox_goodwin", _aka::_ist_fox_goodwin)
       .value("_ist_trapezoidal_rule_2", _aka::_ist_trapezoidal_rule_2)
       .value("_ist_linear_acceleration", _aka::_ist_linear_acceleration)
       .value("_ist_newmark_beta", _aka::_ist_newmark_beta)
       .value("_ist_generalized_trapezoidal", _aka::_ist_generalized_trapezoidal)
       .export_values();
 
   py::enum_<_aka::SolveConvergenceCriteria>(mod, "SolveConvergenceCriteria")
       .value("_scc_residual", _aka::_scc_residual)
       .value("_scc_solution", _aka::_scc_solution)
       .value("_scc_residual_mass_wgh", _aka::_scc_residual_mass_wgh)
       .export_values();
 
   py::enum_<_aka::CohesiveMethod>(mod, "CohesiveMethod")
       .value("_intrinsic", _aka::_intrinsic)
       .value("_extrinsic", _aka::_extrinsic)
       .export_values();
 
   py::enum_<_aka::GhostType>(mod, "GhostType")
       .value("_not_ghost", _aka::_not_ghost)
       .value("_ghost", _aka::_ghost)
       .value("_casper", _aka::_casper)
       .export_values();
 
   py::enum_<_aka::MeshIOType>(mod, "MeshIOType")
       .value("_miot_auto", _aka::_miot_auto)
       .value("_miot_gmsh", _aka::_miot_gmsh)
       .value("_miot_gmsh_struct", _aka::_miot_gmsh_struct)
       .value("_miot_diana", _aka::_miot_diana)
       .value("_miot_abaqus", _aka::_miot_abaqus)
       .export_values();
 
   py::enum_<_aka::ModelType>(mod, "ModelType")
 
       .value("_model", _aka::ModelType::_model)
       .value("_solid_mechanics_model", _aka::ModelType::_solid_mechanics_model)
       .value("_solid_mechanics_model_cohesive",
              _aka::ModelType::_solid_mechanics_model_cohesive)
       .value("_heat_transfer_model", _aka::ModelType::_heat_transfer_model)
       .value("_structural_mechanics_model",
              _aka::ModelType::_structural_mechanics_model)
       .value("_embedded_model", _aka::ModelType::_embedded_model)
       .export_values();
 
   mod.def("parseInput",
           [](const std::string & input_file) {
             _aka::getStaticParser().parse(input_file);
           },
           "Parse an Akantu input file");
 
   py::class_<_aka::Mesh>(mod, "Mesh")
       .def(py::init<_aka::UInt, const _aka::ID &, const _aka::MemoryID &>(),
            py::arg("spatial_dimension"), py::arg("id") = "mesh",
            py::arg("memory_id") = 0)
       .def("read", &_aka::Mesh::read, py::arg("filename"),
            py::arg("mesh_io_type") = _aka::_miot_auto,
            "read the mesh from a file");
 
+  py::class_<_aka::IntegrationPoint>(mod, "IntegrationPoint");
+
   return mod;
 } // namespace
 
 } // namespace
diff --git a/python/pybind11/py_aka_solid_mechanics_model.cc b/python/pybind11/py_aka_solid_mechanics_model.cc
index 54887bb10..289ed4f6a 100644
--- a/python/pybind11/py_aka_solid_mechanics_model.cc
+++ b/python/pybind11/py_aka_solid_mechanics_model.cc
@@ -1,51 +1,82 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
+#include "non_linear_solver.hh"
 #include "solid_mechanics_model.hh"
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 namespace _aka = akantu;
 
 namespace {
 
 /* -------------------------------------------------------------------------- */
 
 py::module & register_solid_mechanics_models(py::module & mod) {
 
   py::class_<_aka::SolidMechanicsModelOptions>(mod,
                                                "SolidMechanicsModelOptions")
       .def(py::init<_aka::AnalysisMethod>(),
            py::arg("analysis_method") = _aka::_explicit_lumped_mass);
 
   py::class_<_aka::Model>(mod, "Model")
       .def("setBaseName", &_aka::Model::setBaseName)
       .def("addDumpFieldVector", &_aka::Model::addDumpFieldVector)
-      .def("addDumpField", &_aka::Model::addDumpField);
+      .def("addDumpField", &_aka::Model::addDumpField)
+      .def("dump", &_aka::Model::dump);
 
-  py::class_<_aka::SolidMechanicsModel, _aka::Model>(mod, "SolidMechanicsModel")
+  py::class_<_aka::NonLinearSolver>(mod, "NonLinearSolver")
+      .def("set",
+           [](_aka::NonLinearSolver & self, const std::string & id,
+              const _aka::Real & val) {
+             if (id == "max_iterations")
+               self.set(id, int(val));
+             else if (id == "convergence_type")
+               self.set(id, akantu::SolveConvergenceCriteria(_aka::UInt(val)));
+             else
+               self.set(id, val);
+           })
+      .def("set", &_aka::NonLinearSolver::set<_aka::SolveConvergenceCriteria>);
+
+  py::class_<_aka::ModelSolver>(mod, "ModelSolver")
+      .def("getNonLinearSolver",
+           (_aka::NonLinearSolver & (_aka::ModelSolver::*)(const _aka::ID &)) &
+               _aka::ModelSolver::getNonLinearSolver,
+           py::arg("solver_id") = "", py::return_value_policy::reference)
+      .def("solveStep", &_aka::ModelSolver::solveStep,
+           py::arg("solver_id") = "");
+
+  py::class_<_aka::SolidMechanicsModel, _aka::Model, _aka::ModelSolver>(
+      mod, "SolidMechanicsModel")
       .def(py::init<_aka::Mesh &, _aka::UInt, const _aka::ID &,
                     const _aka::MemoryID &, const _aka::ModelType>(),
            py::arg("mesh"),
            py::arg("spatial_dimension") = _aka::_all_dimensions,
            py::arg("id") = "solid_mechanics_model", py::arg("memory_id") = 0,
            py::arg("model_type") = _aka::ModelType::_solid_mechanics_model)
       .def("initFull",
            [](_aka::SolidMechanicsModel & self,
               const _aka::SolidMechanicsModelOptions & options) {
              self.initFull(options);
            },
            py::arg("options") = _aka::SolidMechanicsModelOptions())
-      .def("applyBC", [](_aka::SolidMechanicsModel & self,
-                         _aka::BC::DirichletFunctor & func,
-                         const std::string & element_group) {
-        self.applyBC(func, element_group);
-      });
+      .def("applyBC",
+           [](_aka::SolidMechanicsModel & self,
+              _aka::BC::DirichletFunctor & func,
+              const std::string & element_group) {
+             self.applyBC(func, element_group);
+           })
+      .def("applyBC",
+           [](_aka::SolidMechanicsModel & self, _aka::BC::NeumannFunctor & func,
+              const std::string & element_group) {
+             self.applyBC(func, element_group);
+           })
+      .def("getExternalForce", &_aka::SolidMechanicsModel::getExternalForce);
 
   return mod;
 } // namespace
 
 } // namespace
diff --git a/python/pybind11/py_akantu.cc b/python/pybind11/py_akantu.cc
index 4582e93d6..188a60d3a 100644
--- a/python/pybind11/py_akantu.cc
+++ b/python/pybind11/py_akantu.cc
@@ -1,17 +1,19 @@
 #include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
 
 namespace py = pybind11;
 
 #include "py_aka_common.cc"
 #include "py_aka_array.cc"
 #include "py_aka_boundary_conditions.cc"
 #include "py_aka_solid_mechanics_model.cc"
 
 PYBIND11_MODULE(py11_akantu, mod) {
   mod.doc() = "Akantu python interface";
 
   register_enums(mod);
   register_arrays(mod);
   register_boundary_conditions(mod);
   register_solid_mechanics_models(mod);
+
 } // Module akantu