diff --git a/python/py_aka_array.hh b/python/py_aka_array.hh
index 82ab5de1e..3fa7956b6 100644
--- a/python/py_aka_array.hh
+++ b/python/py_aka_array.hh
@@ -1,240 +1,244 @@
 /* -------------------------------------------------------------------------- */
 #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() override final {}
 
     // allocate the memory
     void allocate(UInt /*size*/, UInt /*nb_component*/) override final {}
 
     // allocate and initialize the memory
     void allocate(UInt /*size*/, UInt /*nb_component*/,
                   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();
     }
 
-    ~ArrayProxy() {
-      this->values = nullptr;
-    }
+    ~ArrayProxy() { this->values = nullptr; }
 
-    void resize(UInt /*size*/, const T & /*val */) override final {
-      AKANTU_EXCEPTION("cannot resize a temporary array");
+    void resize(UInt size, const T & val) override final {
+      if (size != this->size())
+        AKANTU_EXCEPTION("cannot resize a temporary array");
+      std::fill(this->begin(), this->end(), val);
     }
 
-    void resize(UInt /*new_size*/) override final {
-      AKANTU_EXCEPTION("cannot resize a temporary array");
+    void resize(UInt new_size) override final {
+      if (new_size != this->size())
+        AKANTU_EXCEPTION("cannot resize a temporary array");
     }
 
     void reserve(UInt /*size*/, UInt /*new_size*/) override 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> *) {
+  decltype(auto) create_proxy(array_type_t<_aka::Vector<T>> & ref,
+                              const _aka::Vector<T> *) {
     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> *) {
+  decltype(auto) create_proxy(array_type_t<_aka::Matrix<T>> & ref,
+                              const _aka::Matrix<T> *) {
     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> *) {
+  decltype(auto) create_proxy(array_type_t<_aka::Array<T>> & ref,
+                              const _aka::Array<T> *) {
     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.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 [[gnu::visibility("default")]] 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
+      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_mesh.cc b/python/py_mesh.cc
index fd0a6096f..7a3ff9783 100644
--- a/python/py_mesh.cc
+++ b/python/py_mesh.cc
@@ -1,68 +1,72 @@
 /* -------------------------------------------------------------------------- */
 #include "aka_config.hh"
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <mesh.hh>
 #include <mesh_utils.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 void register_mesh(py::module & mod) {
   py::class_<MeshData>(mod, "MeshData")
       .def(
           "getElementalDataUInt",
           [](MeshData & _this, const ID & name) -> ElementTypeMapArray<UInt> & {
             return _this.getElementalData<UInt>(name);
           },
           py::return_value_policy::reference);
 
   py::class_<Mesh, GroupManager, Dumpable, MeshData>(mod, "Mesh",
                                                      py::multiple_inheritance())
       .def(py::init<UInt, const ID &, const MemoryID &>(),
            py::arg("spatial_dimension"), py::arg("id") = "mesh",
            py::arg("memory_id") = 0)
       .def("read", &Mesh::read, py::arg("filename"),
            py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file")
       .def(
           "getNodes",
           [](Mesh & self) -> decltype(auto) { return self.getNodes(); },
           py::return_value_policy::reference)
       .def("getNbNodes", &Mesh::getNbNodes)
       .def(
           "getConnectivity",
           [](Mesh & self, const ElementType & type) -> decltype(auto) {
             return self.getConnectivity(type);
           },
           py::return_value_policy::reference)
       .def("distribute", [](Mesh & self) { self.distribute(); })
+      .def("makePeriodic",
+           [](Mesh & self, const SpatialDirection & direction) {
+             self.makePeriodic(direction);
+           })
       .def(
           "getNbElement",
           [](Mesh & self, const UInt spatial_dimension,
              const GhostType & ghost_type, const 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, const ElementType & type,
              const 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);
       });
 
   /* ------------------------------------------------------------------------ */
   py::class_<MeshUtils>(mod, "MeshUtils")
       .def_static("buildFacets", &MeshUtils::buildFacets);
 }
 } // namespace akantu
diff --git a/python/py_model.cc b/python/py_model.cc
index ed0dff0cb..45a820343 100644
--- a/python/py_model.cc
+++ b/python/py_model.cc
@@ -1,96 +1,106 @@
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <model.hh>
 #include <non_linear_solver.hh>
 #include <sparse_matrix_aij.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 
 void register_model(py::module & mod) {
   py::class_<SparseMatrix>(mod, "SparseMatrix")
       .def("getMatrixType", &SparseMatrix::getMatrixType)
       .def("size", &SparseMatrix::size)
       .def("clear", &SparseMatrix::clear)
       .def("saveProfile", &SparseMatrix::saveProfile)
       .def("saveMatrix", &SparseMatrix::saveMatrix)
       .def(
           "add", [](SparseMatrix & self, UInt i, UInt j) { self.add(i, j); },
           "Add entry in the profile")
       .def(
           "add",
           [](SparseMatrix & self, UInt i, UInt j, Real value) {
             self.add(i, j, value);
           },
           "Add the value to the matrix")
       .def("__call__", [](const SparseMatrix & self, UInt i, UInt j) {
         return self(i, j);
       });
 
   py::class_<SparseMatrixAIJ, SparseMatrix>(mod, "SparseMatrixAIJ")
       .def("getIRN", &SparseMatrixAIJ::getIRN)
       .def("getJCN", &SparseMatrixAIJ::getJCN)
       .def("getA", &SparseMatrixAIJ::getA);
 
+  py::class_<SolverVector>(mod, "SolverVector");
+
   py::class_<DOFManager>(mod, "DOFManager")
       .def("getMatrix", &DOFManager::getMatrix,
            py::return_value_policy::reference)
       .def(
           "getNewMatrix",
           [](DOFManager & self, const std::string & name,
              const std::string & matrix_to_copy_id) -> decltype(auto) {
             return self.getNewMatrix(name, matrix_to_copy_id);
           },
-          py::return_value_policy::reference);
+          py::return_value_policy::reference)
+      .def(
+          "getResidual",
+          [](DOFManager & self) -> decltype(auto) {
+            return self.getResidual();
+          },
+          py::return_value_policy::reference)
+      .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs)
+      .def("assembleToResidual", &DOFManager::assembleToResidual);
 
   py::class_<NonLinearSolver>(mod, "NonLinearSolver")
       .def(
           "set",
           [](NonLinearSolver & self, const std::string & id, const Real & val) {
             if (id == "max_iterations")
               self.set(id, int(val));
             else
               self.set(id, val);
           })
       .def("set",
            [](NonLinearSolver & self, const std::string & id,
               const SolveConvergenceCriteria & val) { self.set(id, val); });
 
   py::class_<ModelSolver, Parsable>(mod, "ModelSolver",
                                     py::multiple_inheritance())
       .def("getNonLinearSolver",
            (NonLinearSolver & (ModelSolver::*)(const ID &)) &
                ModelSolver::getNonLinearSolver,
            py::arg("solver_id") = "", py::return_value_policy::reference)
       .def("solveStep", [](ModelSolver & self) { self.solveStep(); })
       .def("solveStep", [](ModelSolver & self, const ID & solver_id) {
         self.solveStep(solver_id);
       });
 
   py::class_<Model, ModelSolver>(mod, "Model", py::multiple_inheritance())
       .def("setBaseName", &Model::setBaseName)
       .def("getFEEngine", &Model::getFEEngine, py::arg("name") = "",
            py::return_value_policy::reference)
       .def("getFEEngineBoundary", &Model::getFEEngine, py::arg("name") = "",
            py::return_value_policy::reference)
       .def("addDumpFieldVector", &Model::addDumpFieldVector)
       .def("addDumpField", &Model::addDumpField)
       .def("setBaseNameToDumper", &Model::setBaseNameToDumper)
       .def("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper)
       .def("addDumpFieldToDumper", &Model::addDumpFieldToDumper)
       .def("dump", &Model::dump)
       .def("initNewSolver", &Model::initNewSolver)
       .def("getDOFManager", &Model::getDOFManager,
            py::return_value_policy::reference);
 }
 
 } // namespace akantu