diff --git a/python/cast.hh b/python/cast.hh
index 77ed63a..beaa489 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,208 +1,210 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __CAST_HH__
 #define __CAST_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid_base.hh"
-#include "surface.hh"
 #include "numpy.hh"
+#include "surface.hh"
 /* -------------------------------------------------------------------------- */
-#include <pybind11/cast.h>
 #include <boost/preprocessor/seq.hpp>
+#include <pybind11/cast.h>
 /* -------------------------------------------------------------------------- */
 
 namespace pybind11 {
 
 // Format descriptor necessary for correct wrap of tamaas complex type
 template <typename T>
 struct format_descriptor<
     tamaas::complex<T>, detail::enable_if_t<std::is_floating_point<T>::value>> {
   static constexpr const char c = format_descriptor<T>::c;
   static constexpr const char value[3] = {'Z', c, '\0'};
   static std::string format() { return std::string(value); }
 };
 
 #ifndef PYBIND11_CPP17
 
 template <typename T>
 constexpr const char format_descriptor<
     tamaas::complex<T>,
     detail::enable_if_t<std::is_floating_point<T>::value>>::value[3];
 
 #endif
 
 namespace detail {
 
 // declare tamaas complex as a complex type for pybind11
 template <typename T>
 struct is_complex<tamaas::complex<T>> : std::true_type {};
 
 template <typename T>
 struct is_fmt_numeric<tamaas::complex<T>,
                       detail::enable_if_t<std::is_floating_point<T>::value>>
     : std::true_type {
   static constexpr int index = is_fmt_numeric<T>::index + 3;
 };
 
+static handle policy_switch(return_value_policy policy, handle parent) {
+  switch (policy) {
+  case return_value_policy::copy:
+  case return_value_policy::move:
+    return handle();
+  case return_value_policy::reference:
+    return none();
+  case return_value_policy::reference_internal:
+    return parent;
+  default:
+    TAMAAS_EXCEPTION("Policy is not handled");
+  }
+}
+
+template <class array, typename T, unsigned int dim>
+handle grid_to_python(const tamaas::Grid<T, dim>& grid,
+                      return_value_policy policy, handle parent) {
+  parent = policy_switch(policy, parent);  // reusing variable
+  std::vector<unsigned int> sizes(dim);
+  std::copy(grid.sizes().begin(), grid.sizes().end(), sizes.begin());
+
+  if (grid.getNbComponents() != 1)
+    sizes.push_back(grid.getNbComponents());
+
+  return array(sizes, grid.getInternalData(), parent).release();
+}
 
 /**
  * Type caster for grid classes
  * inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
  * and pybind11/eigen.h
  */
 template <template <typename, unsigned int> class G, typename T,
           unsigned int dim>
 struct type_caster<G<T, dim>> {
   using type = G<T, dim>;
-  template <typename U>
-  using array_type = array_t<U, array::c_style | array::forcecast>;
+  using array_type =
+      array_t<typename type::value_type, array::c_style | array::forcecast>;
 
 public:
   PYBIND11_TYPE_CASTER(type, _("GridWrap<T, dim>"));
 
   /**
    * Conversion part 1 (Python->C++): convert a PyObject into a grid
    * 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))
+    if (!convert && !array_type::check_(src))
       return false;
 
-    auto buf = array_type<typename type::value_type>::ensure(src);
+    auto buf = array_type::ensure(src);
     value.move(tamaas::wrap::GridNumpy<G<T, dim>>(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*/) {
-    if (src.getNbComponents() == 1) {
-      // none() passed as parent to get a correct nocopy
-      auto a = array_type<typename type::value_type>(
-          src.sizes(), src.getInternalData(), none());
-      return a.release();
-    } else {
-      std::array<unsigned int, dim+1> sizes;
-      std::copy_n(src.sizes().begin(), dim, sizes.begin());
-      sizes.back() = src.getNbComponents();
-      // none() passed as parent to get a correct nocopy
-      auto a = array_type<typename type::value_type>(
-          sizes, src.getInternalData(), none());
-      return a.release();
-    }
+  static handle cast(const type& src, return_value_policy policy,
+                     handle parent) {
+    return grid_to_python<array_type>(src, policy, parent);
   }
 };
 
 /**
  * Type caster for GridBase classes
  */
 template <typename T>
 struct type_caster<tamaas::GridBase<T>> {
   using type = tamaas::GridBase<T>;
-  template <typename U>
-  using array_type = array_t<U, array::c_style | array::forcecast>;
+  using array_type =
+      array_t<typename type::value_type, array::c_style | array::forcecast>;
 
 public:
   PYBIND11_TYPE_CASTER(type, _("GridBaseWrap<T>"));
 
   bool load(handle src, bool convert) {
-    if (!convert && !array_type<T>::check_(src))
+    if (!convert && !array_type::check_(src))
       return false;
 
-    auto buf = array_type<T>::ensure(src);
+    auto buf = array_type::ensure(src);
     value.move(tamaas::wrap::GridBaseNumpy<T>(buf));
 
     return true;
   }
 
-  static handle cast(const type& src, return_value_policy /* policy */,
-                     handle /*parent*/) {
+  static handle cast(const type& src, return_value_policy policy,
+                     handle parent) {
 #define GRID_BASE_CASE(unused1, unused2, dim)                                  \
   case dim: {                                                                  \
     const tamaas::Grid<T, dim>& conv = dynamic_cast<decltype(conv)>(src);      \
-    if (src.getNbComponents() != 1) {                                          \
-      std::array<unsigned int, dim + 1> sizes;                                 \
-      std::copy_n(conv.sizes().begin(), dim, sizes.begin());                   \
-      sizes.back() = src.getNbComponents();                                    \
-      return array_type<T>(sizes, src.getInternalData(), none()).release();    \
-    } else {                                                                   \
-      return array_type<T>(conv.sizes(), src.getInternalData(), none())        \
-          .release();                                                          \
-    }                                                                          \
+    return grid_to_python<array_type>(conv, policy, parent);                   \
   }
 
     switch (src.getDimension()) {
       BOOST_PP_SEQ_FOR_EACH(GRID_BASE_CASE, ~, (1)(2)(3));
     default:
       return nullptr;
     }
 #undef GRID_BASE_CASE
   }
 };
 
 /**
  * Type caster for surface class
  */
 template <typename T>
 struct type_caster<tamaas::Surface<T>> {
   using type = tamaas::Surface<T>;
-  template <typename U>
-  using array_type = array_t<U, array::c_style | array::forcecast>;
+  using array_type =
+      array_t<typename type::value_type, array::c_style | array::forcecast>;
 
 public:
   PYBIND11_TYPE_CASTER(type, _("SurfaceWrap<T, dim>"));
 
   bool load(handle src, bool convert) {
-    if (!convert && !array_type<T>::check_(src))
+    if (!convert && !array_type::check_(src))
       return false;
 
-    auto buf = array_type<T>::ensure(src);
+    auto buf = array_type::ensure(src);
     value.move(tamaas::wrap::SurfaceNumpy<type>(buf));
 
     return true;
   }
 
-  static handle cast(const type& src, return_value_policy /* policy */,
-                     handle /*parent*/) {
-    auto a = array_type<T>(src.sizes(), src.getInternalData(), none());
-    return a.release();
+  static handle cast(const type& src, return_value_policy policy,
+                     handle parent) {
+    return grid_to_python<array_type>(src, policy, parent);
   }
 };
 }  // namespace detail
 }  // namespace pybind11
 
-#endif // __CAST_HH__
+#endif  // __CAST_HH__
diff --git a/python/wrap/bem.cpp b/python/wrap/bem.cpp
index 815c14c..7933fe6 100644
--- a/python/wrap/bem.cpp
+++ b/python/wrap/bem.cpp
@@ -1,111 +1,119 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
-#include "wrap.hh"
 #include "bem_fft_base.hh"
-#include "bem_polonski.hh"
-#include "bem_kato.hh"
 #include "bem_gigipol.hh"
 #include "bem_grid.hh"
-#include "bem_grid_polonski.hh"
 #include "bem_grid_condat.hh"
+#include "bem_grid_polonski.hh"
+#include "bem_kato.hh"
+#include "bem_polonski.hh"
+#include "wrap.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace wrap {
 
 using namespace py::literals;
 
 /* -------------------------------------------------------------------------- */
 void wrapBEM(py::module& mod) {
   py::class_<BemFFTBase>(mod, "BemFFTBase")
       .def("computeDisplacementsFromTractions",
            &BemFFTBase::computeDisplacementsFromTractions)
       .def("computeTractionsFromDisplacements",
            &BemFFTBase::computeTractionsFromDisplacements)
       .def("applyInfluenceFunctions", &BemFFTBase::applyInfluenceFunctions)
       .def("applyInverseInfluenceFunctions",
            &BemFFTBase::applyInverseInfluenceFunctions)
       .def("computeSpectralInfluenceOverDisplacements",
            &BemFFTBase::computeSpectralInfluenceOverDisplacement)
       .def("computeTrueDisplacements", &BemFFTBase::computeTrueDisplacements)
       .def("computeGaps", &BemFFTBase::computeGaps)
       .def("getSpectralInfluenceOverDisplacement",
            &BemFFTBase::getSpectralInfluenceOverDisplacement)
-      .def("getTractions", &BemFFTBase::getTractions)
-      .def("getDisplacements", &BemFFTBase::getDisplacements)
-      .def("getTrueDisplacements", &BemFFTBase::getTrueDisplacements)
-      .def("getGap", &BemFFTBase::getGap)
+      .def("getTractions", &BemFFTBase::getTractions,
+           py::return_value_policy::reference_internal)
+      .def("getDisplacements", &BemFFTBase::getDisplacements,
+           py::return_value_policy::reference_internal)
+      .def("getTrueDisplacements", &BemFFTBase::getTrueDisplacements,
+           py::return_value_policy::reference_internal)
+      .def("getGap", &BemFFTBase::getGap,
+           py::return_value_policy::reference_internal)
       .def("setEffectiveModulus", &BemFFTBase::setEffectiveModulus, "E_star"_a)
       .def("getEffectiveModulus", &BemFFTBase::getEffectiveModulus)
       .def("setMaxIterations", &BemFFTBase::setMaxIterations)
       .def("getMaxIterations", &BemFFTBase::getMaxIterations)
       .def("setDumpFreq", &BemFFTBase::setDumpFreq)
       .def("getDumpFreq", &BemFFTBase::getDumpFreq);
 
   py::class_<BemPolonski, BemFFTBase>(mod, "BemPolonski")
       .def(py::init<Surface<Real>&>())
       .def("computeEquilibrium", &BemPolonski::computeEquilibrium)
       .def("computeEquilibriuminit", &BemPolonski::computeEquilibriuminit)
       .def("computeMeanGapsInContact", &BemPolonski::computeMeanGapsInContact)
       .def("setMaxPressure", &BemPolonski::setMaxPressure)
       .def("getNbIterations", &BemPolonski::getNbIterations);
 
   py::class_<BemKato, BemFFTBase>(mod, "BemKato")
       .def(py::init<Surface<Real>&>())
       .def("computeEquilibrium", &BemKato::computeEquilibrium)
       .def("setMaxPressure", &BemKato::setMaxPressure);
 
   py::class_<BemGigipol, BemFFTBase>(mod, "BemGigiPol")
       .def(py::init<Surface<Real>&>())
       .def("computeEquilibrium", &BemGigipol::computeEquilibrium)
       .def("computeEquilibrium2", &BemGigipol::computeEquilibrium2)
       .def("computeEquilibriuminit", &BemGigipol::computeEquilibriuminit)
       .def("computeEquilibrium2init", &BemGigipol::computeEquilibrium2init)
       .def("computeTractionsFromDisplacements",
            &BemGigipol::computeTractionsFromDisplacements);
 
   py::class_<BemGrid>(mod, "BemGrid")
-    .def("computeDisplacementsFromTractions", &BemGrid::computeDisplacementsFromTractions)
-    .def("computeSurfaceNormals", &BemGrid::computeSurfaceNormals)
-    .def("computeInfluence", &BemGrid::computeInfluence)
-    .def("setElasticity", &BemGrid::setElasticity)
-    .def("getTractions", &BemGrid::getTractions)
-    .def("getDisplacements", &BemGrid::getDisplacements)
-    .def("getSurfaceNormals", &BemGrid::getSurfaceNormals);
+      .def("computeDisplacementsFromTractions",
+           &BemGrid::computeDisplacementsFromTractions)
+      .def("computeSurfaceNormals", &BemGrid::computeSurfaceNormals)
+      .def("computeInfluence", &BemGrid::computeInfluence)
+      .def("setElasticity", &BemGrid::setElasticity)
+      .def("getTractions", &BemGrid::getTractions,
+           py::return_value_policy::reference_internal)
+      .def("getDisplacements", &BemGrid::getDisplacements,
+           py::return_value_policy::reference_internal)
+      .def("getSurfaceNormals", &BemGrid::getSurfaceNormals,
+           py::return_value_policy::reference_internal);
 
   py::class_<BemGridPolonski, BemGrid>(mod, "BemGridPolonski")
-    .def(py::init<Surface<Real>&>());
+      .def(py::init<Surface<Real>&>());
 
   py::class_<BemGridCondat, BemGrid>(mod, "BemGridCondat")
-    .def(py::init<Surface<Real>&>())
-    .def("computeEquilibrium", &BemGridCondat::computeEquilibrium);
+      .def(py::init<Surface<Real>&>())
+      .def("computeEquilibrium", &BemGridCondat::computeEquilibrium);
 }
 
-} // namespace wrap
+}  // namespace wrap
 
 __END_TAMAAS__
diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index d636ceb..74018b0 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,236 +1,242 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "model.hh"
 #include "adhesion_functional.hh"
 #include "functional.hh"
 #include "integral_operator.hh"
 #include "model_dumper.hh"
 #include "model_extensions.hh"
 #include "model_factory.hh"
 #include "numpy.hh"
 #include "residual.hh"
 #include "wrap.hh"
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace wrap {
 
 using namespace py::literals;
 
 /// Wrap functional classes
 void wrapFunctionals(py::module& mod) {
   py::class_<functional::Functional, functional::wrap::PyFunctional> func(
       mod, "Functional");
   func.def(py::init<>())
       .def("computeF", &functional::Functional::computeF)
       .def("computeGradF", &functional::Functional::computeGradF);
 
   py::class_<functional::AdhesionFunctional> adh(mod, "AdhesionFunctional",
                                                  func);
   adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
                    &functional::AdhesionFunctional::setParameters)
       // legacy wrapper code
       .def("setParameters", &functional::AdhesionFunctional::setParameters);
 
   py::class_<functional::ExponentialAdhesionFunctional>(
       mod, "ExponentialAdhesionFunctional", adh)
       .def(py::init<const GridBase<Real>&>(), "surface"_a);
   py::class_<functional::MaugisAdhesionFunctional>(
       mod, "MaugisAdhesionFunctional", adh)
       .def(py::init<const GridBase<Real>&>(), "surface"_a);
   py::class_<functional::SquaredExponentialAdhesionFunctional>(
       mod, "SquaredExponentialAdhesionFunctional", adh)
       .def(py::init<const GridBase<Real>&>(), "surface"_a);
 }
 
 template <typename T>
 std::unique_ptr<GridBase<T>> instanciateFromNumpy(numpy<T>& num) {
   std::unique_ptr<GridBase<T>> result = nullptr;
 
   switch (num.ndim()) {
   case 2:
     result = std::make_unique<GridNumpy<Grid<T, 1>>>(num);
     return result;
   case 3:
     result = std::make_unique<GridNumpy<Grid<T, 2>>>(num);
     return result;
   case 4:
     result = std::make_unique<GridNumpy<Grid<T, 3>>>(num);
     return result;
   default:
     TAMAAS_EXCEPTION("instanciateFromNumpy expects the last dimension of numpy "
                      "array to be the number of components");
   }
 }
 
 /// Wrap IntegralOperator
 void wrapIntegralOperator(py::module& mod) {
   py::class_<IntegralOperator>(mod, "IntegralOperator")
       .def("apply",
            [](IntegralOperator& op, numpy<Real> input, numpy<Real> output) {
              auto in = instanciateFromNumpy<Real>(input);
              auto out = instanciateFromNumpy<Real>(output);
 
              op.apply(*in, *out);
            })
       .def("updateFromModel", &IntegralOperator::updateFromModel)
       .def("getModel", &IntegralOperator::getModel,
            py::return_value_policy::reference)
       .def("getKind", &IntegralOperator::getKind)
       .def("getType", &IntegralOperator::getType);
 }
 
 /// Wrap BEEngine classes
 void wrapBEEngine(py::module& mod) {
   py::class_<BEEngine>(mod, "BEEngine")
       .def("solveNeumann", &BEEngine::solveNeumann)
       .def("solveDirichlet", &BEEngine::solveDirichlet)
       .def("getModel", &BEEngine::getModel, py::return_value_policy::reference);
 }
 
 /// Wrap Models
 void wrapModelClass(py::module& mod) {
   py::enum_<model_type>(mod, "model_type")
       .value("basic_1d", model_type::basic_1d)
       .value("basic_2d", model_type::basic_2d)
       .value("surface_1d", model_type::surface_1d)
       .value("surface_2d", model_type::surface_2d)
       .value("volume_1d", model_type::volume_1d)
       .value("volume_2d", model_type::volume_2d);
 
   py::class_<Model>(mod, "Model")
       .def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a)
       .def_property("E", &Model::getYoungModulus, &Model::setYoungModulus)
       .def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio)
       .def("getHertzModulus", &Model::getHertzModulus)
       .def("getYoungModulus", &Model::getYoungModulus)
       .def("getShearModulus", &Model::getShearModulus)
       .def("getPoissonRatio", &Model::getPoissonRatio)
-      .def("getTraction", (GridBase<Real> & (Model::*)()) & Model::getTraction)
+      .def("getTraction", (GridBase<Real> & (Model::*)()) & Model::getTraction,
+           py::return_value_policy::reference_internal)
       .def("getDisplacement",
-           (GridBase<Real> & (Model::*)()) & Model::getDisplacement)
+           (GridBase<Real> & (Model::*)()) & Model::getDisplacement,
+           py::return_value_policy::reference_internal)
       .def("getSystemSize", &Model::getSystemSize)
       .def("getDiscretization", &Model::getDiscretization)
       .def("getBoundarySystemSize", &Model::getBoundarySystemSize)
       .def("getBoundaryDiscretization", &Model::getBoundaryDiscretization)
       .def("solveNeumann", &Model::solveNeumann)
       .def("solveDirichlet", &Model::solveDirichlet)
       .def("dump", &Model::dump)
       .def("addDumper", &Model::addDumper, "dumper"_a)
       .def("setDumper",
            [](Model& m, const ModelDumper&) {
              TAMAAS_EXCEPTION("setDumper() is not a member of Model; use "
                               "addDumper() instead");
            })
       .def("getBEEngine", &Model::getBEEngine,
-           py::return_value_policy::reference)
+           py::return_value_policy::reference_internal)
       .def("getIntegralOperator", &Model::getIntegralOperator,
-           "operator_name"_a, py::return_value_policy::reference)
-      .def("getField", &Model::getField, "field_name"_a)
+           "operator_name"_a, py::return_value_policy::reference_internal)
+      .def("getField", &Model::getField, "field_name"_a,
+           py::return_value_policy::reference_internal)
       .def("applyElasticity",
            [](Model& model, numpy<Real> stress, numpy<Real> strain) {
              auto out = instanciateFromNumpy<Real>(stress);
              auto in = instanciateFromNumpy<Real>(strain);
 
              model.applyElasticity(*out, *in);
            })
       .def("__repr__", [](const Model& m) {
         std::stringstream ss;
         ss << m;
         return ss.str();
       });
 
   py::class_<ModelDumper, PyModelDumper, std::shared_ptr<ModelDumper>>(
       mod, "ModelDumper")
       .def(py::init<>())
       .def("dump", &ModelDumper::dump, "model"_a)
       .def("__lshift__",
            [](ModelDumper& dumper, Model& model) { dumper << model; });
 }
 
 /// Wrap factory for models
 void wrapModelFactory(py::module& mod) {
   py::class_<ModelFactory>(mod, "ModelFactory")
       .def_static("createModel", &ModelFactory::createModel, "model_type"_a,
                   "system_size"_a, "discretization"_a)
       .def_static("createResidual", &ModelFactory::createResidual,
                   "model_type"_a, "hardening"_a, "sigma_y"_a);
 }
 
 /// Wrap residual class
 void wrapResidual(py::module& mod) {
   // TODO adapt to n-dim
   py::class_<Residual, PyResidual>(mod, "Residual")
       .def(py::init<Model*>())
       .def("computeResidual",
            [](Residual& res, numpy<Real>& x) {
              auto in = instanciateFromNumpy(x);
              res.computeResidual(*in);
            })
       .def("computeStress",
            [](Residual& res, numpy<Real>& x) {
              auto in = instanciateFromNumpy(x);
              res.computeStress(*in);
            })
       .def("updateState",
            [](Residual& res, numpy<Real>& x) {
              auto in = instanciateFromNumpy(x);
              res.updateState(*in);
            })
       .def("computeResidualDisplacement",
            [](Residual& res, numpy<Real>& x) {
              auto in = instanciateFromNumpy(x);
              res.computeResidualDisplacement(*in);
            })
       .def("applyTangent",
            [](Residual& res, numpy<Real>& output, numpy<Real>& input,
               numpy<Real>& current_strain_inc) {
              auto out = instanciateFromNumpy(output);
              auto in = instanciateFromNumpy(input);
              auto inc = instanciateFromNumpy(current_strain_inc);
              res.applyTangent(*out, *in, *inc);
            },
            "output"_a, "input"_a, "current_strain_increment"_a)
-      .def("getVector", &Residual::getVector)
-      .def("getPlasticStrain", &Residual::getPlasticStrain)
-      .def("getStress", &Residual::getStress);
+      .def("getVector", &Residual::getVector,
+           py::return_value_policy::reference_internal)
+      .def("getPlasticStrain", &Residual::getPlasticStrain,
+           py::return_value_policy::reference_internal)
+      .def("getStress", &Residual::getStress,
+           py::return_value_policy::reference_internal);
 }
 
 void wrapModel(py::module& mod) {
   wrapBEEngine(mod);
   wrapModelClass(mod);
   wrapModelFactory(mod);
   wrapFunctionals(mod);
   wrapResidual(mod);
   wrapIntegralOperator(mod);
 }
 
 }  // namespace wrap
 
 __END_TAMAAS__
diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp
index 0d7591c..695b1d4 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,134 +1,135 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
-#include "wrap.hh"
-#include "contact_solver.hh"
-#include "polonsky_keer_rey.hh"
-#include "kato.hh"
 #include "beck_teboulle.hh"
 #include "condat.hh"
-#include "polonsky_keer_tan.hh"
-#include "epic.hh"
+#include "contact_solver.hh"
 #include "ep_solver.hh"
+#include "epic.hh"
+#include "kato.hh"
+#include "polonsky_keer_rey.hh"
+#include "polonsky_keer_tan.hh"
+#include "wrap.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace wrap {
 
 /* -------------------------------------------------------------------------- */
 using namespace py::literals;
 
 class PyEPSolver : public EPSolver {
 public:
   using EPSolver::EPSolver;
-  void solve() override {
-    PYBIND11_OVERLOAD_PURE(void, EPSolver, solve);
-  }
+  void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); }
 
   void updateState() override {
     PYBIND11_OVERLOAD(void, EPSolver, updateState);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 void wrapSolvers(py::module& mod) {
   py::class_<ContactSolver>(mod, "ContactSolver")
       // .def(py::init<Model&, const GridBase<Real>&, Real>(), "model"_a,
       //      "surface"_a, "tolerance"_a)
       .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a)
       .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a)
       .def("addFunctionalTerm", &ContactSolver::addFunctionalTerm)
       .def("solve", py::overload_cast<std::vector<Real>>(&ContactSolver::solve),
            "target_force"_a)
       .def("solve", py::overload_cast<Real>(&ContactSolver::solve),
            "target_normal_pressure"_a);
 
   py::class_<PolonskyKeerRey, ContactSolver> pkr(mod, "PolonskyKeerRey");
   pkr.def(py::init<Model&, const GridBase<Real>&, Real, PolonskyKeerRey::type,
                    PolonskyKeerRey::type>(),
           "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a,
           "constraint_type"_a, py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
       .def("computeError", &PolonskyKeerRey::computeError);
 
   py::enum_<PolonskyKeerRey::type>(pkr, "type")
       .value("gap", PolonskyKeerRey::gap)
       .value("pressure", PolonskyKeerRey::pressure)
       .export_values();
 
   py::class_<Kato, ContactSolver> kato(mod, "Kato");
-  kato.def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
-          "model"_a, "surface"_a, "tolerance"_a, "mu"_a,
-          py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
-      .def("solve", &Kato::solve, "p0"_a, "proj_iter"_a=50)
+  kato.def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
+           "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
+           py::keep_alive<1, 3>())
+      .def("solve", &Kato::solve, "p0"_a, "proj_iter"_a = 50)
       .def("solveRelaxed", &Kato::solveRelaxed, "g0"_a)
-      .def("solveRegularized", &Kato::solveRegularized, "p0"_a, "r"_a=0.01)
-      .def("computeCost", &Kato::computeCost, "use_tresca"_a=false);
+      .def("solveRegularized", &Kato::solveRegularized, "p0"_a, "r"_a = 0.01)
+      .def("computeCost", &Kato::computeCost, "use_tresca"_a = false);
 
   py::class_<BeckTeboulle, ContactSolver> bt(mod, "BeckTeboulle");
-  bt.def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
-          "model"_a, "surface"_a, "tolerance"_a, "mu"_a,
-          py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
+  bt.def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
+         "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
+         py::keep_alive<1, 3>())
       .def("solve", &BeckTeboulle::solve, "p0"_a)
       .def("computeCost", &BeckTeboulle::computeCost);
 
   py::class_<Condat, ContactSolver> cd(mod, "Condat");
-  cd.def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
-          "model"_a, "surface"_a, "tolerance"_a, "mu"_a,
-          py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
-      .def("solve", &Condat::solve, "p0"_a, "grad_step"_a=0.9)
+  cd.def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
+         "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
+         py::keep_alive<1, 3>())
+      .def("solve", &Condat::solve, "p0"_a, "grad_step"_a = 0.9)
       .def("computeCost", &Condat::computeCost);
 
   py::class_<PolonskyKeerTan, ContactSolver> pkt(mod, "PolonskyKeerTan");
-  pkt.def(py::init<Model&, const GridBase<Real>&, Real, Real>(),
-          "model"_a, "surface"_a, "tolerance"_a, "mu"_a,
-          py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
+  pkt.def(py::init<Model&, const GridBase<Real>&, Real, Real>(), "model"_a,
+          "surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
+          py::keep_alive<1, 3>())
       .def("solve", &PolonskyKeerTan::solve, "p0"_a)
       .def("solveTresca", &PolonskyKeerTan::solveTresca, "p0"_a)
-      .def("computeCost", &PolonskyKeerTan::computeCost, "use_tresca"_a=false);
+      .def("computeCost", &PolonskyKeerTan::computeCost,
+           "use_tresca"_a = false);
 
   py::class_<EPSolver, PyEPSolver>(mod, "EPSolver")
       .def(py::init<Residual&>())
       .def("solve", &EPSolver::solve)
-      .def("getStrainIncrement", &EPSolver::getStrainIncrement)
-      .def("getResidual", &EPSolver::getResidual)
+      .def("getStrainIncrement", &EPSolver::getStrainIncrement,
+           py::return_value_policy::reference_internal)
+      .def("getResidual", &EPSolver::getResidual,
+           py::return_value_policy::reference_internal)
       .def("updateState", &EPSolver::updateState);
 
   py::class_<EPICSolver>(mod, "EPICSolver")
       .def(py::init<ContactSolver&, EPSolver&, Real, Real>(),
            "contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10,
            "relaxation"_a = 0.3)
       .def("solve", &EPICSolver::solve, "force"_a)
       .def("solve",
            [](EPICSolver& solver, Real pressure) {
              return solver.solve(std::vector<Real>{pressure});
            },
            "normal_pressure"_a);
 }
 
 }  // namespace wrap
 
 __END_TAMAAS__
diff --git a/python/wrap/surface.cpp b/python/wrap/surface.cpp
index 5975865..68b8e4e 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,184 +1,186 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "isopowerlaw.hh"
 #include "surface_generator_filter.hh"
 #include "surface_generator_filter_fft.hh"
 #include "surface_generator_random_phase.hh"
 #include "wrap.hh"
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /* -------------------------------------------------------------------------- */
 namespace wrap {
 
 using namespace py::literals;
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 class PyFilter : public Filter<dim> {
 public:
   using Filter<dim>::Filter;
 
   // Overriding pure virtual functions
   void
   computeFilter(GridHermitian<Real, dim>& filter_coefficients) const override {
     PYBIND11_OVERLOAD_PURE(void, Filter<dim>, computeFilter,
                            filter_coefficients);
   }
 };
 
 template <UInt dim>
 void wrapFilter(py::module& mod) {
   auto name = makeDimensionName("Filter", dim);
   py::class_<Filter<dim>, PyFilter<dim>>(mod, name.c_str())
       .def(py::init<>())
       .def("computeFilter", &Filter<dim>::computeFilter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 void wrapIsopowerlaw(py::module& mod) {
   std::string name = makeDimensionName("Isopowerlaw", dim);
   py::class_<Isopowerlaw<dim>, Filter<dim>>(mod, name.c_str())
       .def(py::init<>())
       .def_property("q0", &Isopowerlaw<dim>::getQ0, &Isopowerlaw<dim>::setQ0)
       .def_property("q1", &Isopowerlaw<dim>::getQ1, &Isopowerlaw<dim>::setQ1)
       .def_property("q2", &Isopowerlaw<dim>::getQ2, &Isopowerlaw<dim>::setQ2)
       .def_property("hurst", &Isopowerlaw<dim>::getHurst,
                     &Isopowerlaw<dim>::setHurst)
       .def("rmsHeights", &Isopowerlaw<dim>::rmsHeights)
       .def("moments", &Isopowerlaw<dim>::moments)
       .def("alpha", &Isopowerlaw<dim>::alpha)
       .def("rmsSlopes", &Isopowerlaw<dim>::rmsSlopes)
       // legacy wrapper code
       .def("getQ0", &Isopowerlaw<dim>::getQ0,
            py::return_value_policy::reference)
       .def("getQ1", &Isopowerlaw<dim>::getQ1,
            py::return_value_policy::reference)
       .def("getQ2", &Isopowerlaw<dim>::getQ2,
            py::return_value_policy::reference)
       .def("setQ0", &Isopowerlaw<dim>::setQ0,
            py::return_value_policy::reference)
       .def("setQ1", &Isopowerlaw<dim>::setQ1,
            py::return_value_policy::reference)
       .def("setQ2", &Isopowerlaw<dim>::setQ2,
            py::return_value_policy::reference)
       .def("setHurst", &Isopowerlaw<dim>::setHurst)
       .def("getHurst", &Isopowerlaw<dim>::getHurst,
            py::return_value_policy::reference);
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 void wrapSurfaceGenerators(py::module& mod) {
   std::string generator_name = makeDimensionName("SurfaceGenerator", dim);
   py::class_<SurfaceGenerator<dim>>(mod, generator_name.c_str())
-      .def("buildSurface", &SurfaceGenerator<dim>::buildSurface)
+      .def("buildSurface", &SurfaceGenerator<dim>::buildSurface,
+           py::return_value_policy::reference_internal)
       .def("setSizes",
            py::overload_cast<const std::array<UInt, dim>&>(
                &SurfaceGenerator<dim>::setSizes),
            "surface_sizes"_a)
       .def_property("random_seed", &SurfaceGenerator<dim>::getRandomSeed,
                     &SurfaceGenerator<dim>::setRandomSeed);
 
   std::string filter_name = makeDimensionName("SurfaceGeneratorFilter", dim);
   py::class_<SurfaceGeneratorFilter<dim>, SurfaceGenerator<dim>>(
       mod, filter_name.c_str())
       .def(py::init<>())
       .def("setFilter", &SurfaceGeneratorFilter<dim>::setFilter,
            "Set PSD filter", "filter"_a)
       .def("setSpectrum", &SurfaceGeneratorFilter<dim>::setFilter,
            "Set PSD spectrum", "filter"_a)
       // legacy wrapper code
       .def("setRandomSeed", &SurfaceGeneratorFilter<dim>::setRandomSeed,
            "[[Deprecated: use property]]");
 
   std::string random_name =
       makeDimensionName("SurfaceGeneratorRandomPhase", dim);
   py::class_<SurfaceGeneratorRandomPhase<dim>, SurfaceGenerator<dim>>(
       mod, random_name.c_str())
       .def(py::init<>())
       .def("setSpectrum", &SurfaceGeneratorRandomPhase<dim>::setSpectrum,
            "Set PSD Spectrum", "spectrum"_a);
 }
 
 /* -------------------------------------------------------------------------- */
 // Legacy wrap
 void wrapSurfaceGeneratorFilterFFT(py::module& mod) {
   py::class_<SurfaceGeneratorFilterFFT>(mod, "SurfaceGeneratorFilterFFT")
       .def(py::init<>())
       .def("getQ0",
            [](SurfaceGeneratorFilterFFT& f) {
              return smart_pointer<UInt>(&f.getQ0());
            })
       .def("getQ1",
            [](SurfaceGeneratorFilterFFT& f) {
              return smart_pointer<UInt>(&f.getQ1());
            })
       .def("getQ2",
            [](SurfaceGeneratorFilterFFT& f) {
              return smart_pointer<UInt>(&f.getQ2());
            })
       .def("getHurst",
            [](SurfaceGeneratorFilterFFT& f) {
              return smart_pointer<Real>(&f.getHurst());
            })
       .def("analyticalRMS", &SurfaceGeneratorFilterFFT::analyticalRMS)
-      .def("buildSurface", &SurfaceGeneratorFilterFFT::buildSurface)
+      .def("buildSurface", &SurfaceGeneratorFilterFFT::buildSurface,
+           py::return_value_policy::reference_internal)
       .def("getGridSize",
            [](SurfaceGeneratorFilterFFT& f) {
              return smart_pointer<UInt>(&f.getGridSize());
            })
       .def("getRandomSeed",
            [](SurfaceGeneratorFilterFFT& f) {
              return smart_pointer<long>(&f.getRandomSeed());
            })
       .def("getRMS",
            [](SurfaceGeneratorFilterFFT& f) {
              return smart_pointer<Real>(&f.getRMS());
            })
       .def("Init", &SurfaceGeneratorFilterFFT::Init);
 }
 
 /* -------------------------------------------------------------------------- */
 void wrapSurface(py::module& mod) {
   wrapFilter<1>(mod);
   wrapFilter<2>(mod);
   wrapIsopowerlaw<1>(mod);
   wrapIsopowerlaw<2>(mod);
   wrapSurfaceGenerators<1>(mod);
   wrapSurfaceGenerators<2>(mod);
 
   // legacy wrap
   wrapSurfaceGeneratorFilterFFT(mod);
 }
 
 }  // namespace wrap
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__