diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index 61c2097..0c35641 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,261 +1,261 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://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>
 /* -------------------------------------------------------------------------- */
 
 namespace 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(input);
              auto out = instanciateFromNumpy(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)
       .def("registerNeumann", &BEEngine::registerNeumann)
       .def("registerDirichlet", &BEEngine::registerDirichlet);
 }
 
 /// 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_property_readonly("type", &Model::getType)
       .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,
            py::return_value_policy::reference_internal)
       .def("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("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>())
       .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_internal)
       .def("getIntegralOperator", &Model::getIntegralOperator,
            "operator_name"_a, py::return_value_policy::reference_internal)
       .def("registerField",
            [](Model& m, std::string name, numpy<Real> field) {
              auto f = instanciateFromNumpy(field);
              m.registerField(name, std::move(f));
            },
            "field_name"_a, "field"_a, py::keep_alive<1, 3>())
       .def("getField",
            (GridBase<Real> & (Model::*)(const std::string&)) & Model::getField,
            "field_name"_a, py::return_value_policy::reference_internal)
       .def("getFields", &Model::getFields)
       .def("applyElasticity",
            [](Model& model, numpy<Real> stress, numpy<Real> strain) {
              auto out = instanciateFromNumpy(stress);
              auto in = instanciateFromNumpy(strain);
 
              model.applyElasticity(*out, *in);
            })
       .def("__repr__",
            [](const Model& m) {
              std::stringstream ss;
              ss << m;
              return ss.str();
            })
       .def("__getitem__",
            [](const Model& m, std::string key) -> const GridBase<Real>& {
              try {
                return m.getField(key);
              } catch (std::out_of_range&) {
                throw py::key_error();
              }
            },
            py::return_value_policy::reference_internal);
 
   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)
       .def_static("registerVolumeOperators",
                   &ModelFactory::registerVolumeOperators, "model"_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,
            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
 
 }  // namespace tamaas