diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index 0484633..9238cc7 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,554 +1,556 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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 "elastic_functional.hh"
 #include "field_container.hh"
 #include "functional.hh"
 #include "hooke.hh"
 #include "integral_operator.hh"
 #include "model_dumper.hh"
 #include "model_extensions.hh"
 #include "model_factory.hh"
 #include "model_type.hh"
 #include "numpy.hh"
 #include "residual.hh"
 #include "wrap.hh"
 #include <boost/variant.hpp>
 #include <pybind11/attr.h>
 #include <pybind11/detail/common.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 namespace wrap {
 
 using namespace py::literals;
 
 struct model_operator_accessor {
   Model& m;
 
   decltype(auto) get(const std::string& name) {
     return m.getIntegralOperator(name);
   }
 };
 
 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");
   }
 }
 
 void wrapFieldContainer(py::module& mod) {
   py::class_<FieldContainer>(mod, "FieldContainer")
       .def(py::init<>())
       .def("__getitem__", &FieldContainer::at,
            py::return_value_policy::reference_internal)
       .def(
           "__setitem__",
           [](FieldContainer& f, const std::string& key, numpy<Real> grid) {
             f[key] = instanciateFromNumpy(grid);
           },
           py::keep_alive<1, 3>())
       .def(
           "__setitem__",
           [](FieldContainer& f, const std::string& key, numpy<Complex> grid) {
             f[key] = instanciateFromNumpy(grid);
           },
           py::keep_alive<1, 3>())
       .def(
           "__setitem__",
           [](FieldContainer& f, const std::string& key, numpy<Int> grid) {
             f[key] = instanciateFromNumpy(grid);
           },
           py::keep_alive<1, 3>())
       .def(
           "__contains__",
           [](const FieldContainer& m, const std::string& key) {
             const auto fields = m.fields();
             return std::find(fields.begin(), fields.end(), key) != fields.end();
           },
           "Test field existence")
       .def(
           "__iter__",
           [](const FieldContainer& m) {
             const auto& fields = m.fields_map();
             return py::make_key_iterator(fields.cbegin(), fields.cend());
           },
           py::keep_alive<0, 1>(), "Iterator on fields")
       /// Deprecated interface
       .def(
           "registerField",
           [](Model& m, const std::string& name, numpy<Real> field) {
             TAMAAS_DEPRECATE("registerField()", "the [] operator");
             m[name] = instanciateFromNumpy(field);
           },
           "field_name"_a, "field"_a, py::keep_alive<1, 3>())
       .def(
           "getField",
           [](const Model& m, const std::string& name) -> decltype(m.at(name)) {
             TAMAAS_DEPRECATE("getField()", "the [] operator");
             return m.at(name);
           },
           "field_name"_a, py::return_value_policy::reference_internal)
       .def(
           "getFields",
           [](const Model& m) {
             TAMAAS_DEPRECATE("getFields()", "list(model)");
             return m.fields();
           },
           "Return fields list");
 }
 
 /// Wrap functional classes
 void wrapFunctionals(py::module& mod) {
   py::class_<functional::Functional, std::shared_ptr<functional::Functional>,
              functional::wrap::PyFunctional>
       func(mod, "Functional");
   func.def(py::init<>())
       .def("computeF", &functional::Functional::computeF,
            "Compute functional value")
       .def("computeGradF", &functional::Functional::computeGradF,
            "Compute functional gradient");
 
   py::class_<functional::ElasticFunctionalPressure>(
       mod, "ElasticFunctionalPressure", func)
       .def(py::init<const IntegralOperator&, const GridBase<Real>&>());
   py::class_<functional::ElasticFunctionalGap>(mod, "ElasticFunctionalGap",
                                                func)
       .def(py::init<const IntegralOperator&, const GridBase<Real>&>());
 
   py::class_<functional::AdhesionFunctional> adh(mod, "AdhesionFunctional",
                                                  func);
   adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
                    &functional::AdhesionFunctional::setParameters,
                    "Parameters dictionary")
       .def("setParameters", [](functional::AdhesionFunctional& f,
                                const std::map<std::string, Real>& m) {
         TAMAAS_DEPRECATE("setParameters()", "the parameters property");
         f.setParameters(m);
       });
 
   py::class_<functional::ExponentialAdhesionFunctional>(
       mod, "ExponentialAdhesionFunctional", adh,
       "Potential of the form F = -γ·exp(-g/ρ)")
       .def(py::init<const GridBase<Real>&>(), "surface"_a);
   py::class_<functional::MaugisAdhesionFunctional>(
       mod, "MaugisAdhesionFunctional", adh,
       "Cohesive zone potential F = H(g - ρ)·γ/ρ")
       .def(py::init<const GridBase<Real>&>(), "surface"_a);
   py::class_<functional::SquaredExponentialAdhesionFunctional>(
       mod, "SquaredExponentialAdhesionFunctional", adh,
       "Potential of the form F = -γ·exp(-0.5·(g/ρ)²)")
       .def(py::init<const GridBase<Real>&>(), "surface"_a);
 }
 
 /// Wrap IntegralOperator
 void wrapIntegralOperator(py::module& mod) {
-  py::class_<IntegralOperator, FieldContainer>(mod, "IntegralOperator")
+  py::class_<IntegralOperator, PyIntegralOperator,
+             std::shared_ptr<IntegralOperator>>(mod, "IntegralOperator")
+      .def(py::init<Model*>())
       .def("apply",
            [](IntegralOperator& op, numpy<Real> input, numpy<Real> output) {
              TAMAAS_DEPRECATE("apply()", "the () operator");
              auto in = instanciateFromNumpy(input);
              auto out = instanciateFromNumpy(output);
 
              op.apply(*in, *out);
            })
       .def(TAMAAS_DEPRECATE_ACCESSOR(getModel, IntegralOperator, "model"),
            py::return_value_policy::reference)
       .def(TAMAAS_DEPRECATE_ACCESSOR(getKind, IntegralOperator, "kind"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getType, IntegralOperator, "type"))
 
       .def(
           "__call__",
           [](IntegralOperator& op, numpy<Real> input, numpy<Real> output) {
             auto in = instanciateFromNumpy(input);
             auto out = instanciateFromNumpy(output);
 
             op.apply(*in, *out);
           },
           "Apply the integral operator")
       .def("updateFromModel", &IntegralOperator::updateFromModel,
            "Resets internal persistent variables from the model")
 
       .def_property_readonly("kind", &IntegralOperator::getKind)
       .def_property_readonly("model", &IntegralOperator::getModel)
       .def_property_readonly("type", &IntegralOperator::getType)
       .def_property_readonly("shape", &IntegralOperator::matvecShape)
       .def(
           "matvec",
           [](const IntegralOperator& op, numpy<Real> X) -> GridBase<Real> {
             GridBaseNumpy<Real> x(X);
             auto y = op.matvec(x);
             return y;
           },
           py::return_value_policy::move);
 
   py::enum_<integration_method>(mod, "integration_method",
                                 "Integration method used for the computation "
                                 "of volumetric Fourier operators")
       .value("linear", integration_method::linear,
              "No approximation error, O(N₁·N₂·N₃) time complexity, may cause "
              "float overflow/underflow")
       .value("cutoff", integration_method::cutoff,
              "Approximation, O(sqrt(N₁²+N₂²)·N₃²) time complexity, no "
              "overflow/underflow risk");
 }
 
 /// Wrap BEEngine classes
 void wrapBEEngine(py::module& mod) {
   py::class_<BEEngine>(mod, "BEEngine")
       .def("solveNeumann", &BEEngine::solveNeumann)
       .def("solveDirichlet", &BEEngine::solveDirichlet)
       .def("registerNeumann", &BEEngine::registerNeumann)
       .def("registerDirichlet", &BEEngine::registerDirichlet)
       .def(TAMAAS_DEPRECATE_ACCESSOR(getModel, BEEngine, "model"),
            py::return_value_policy::reference)
 
       .def_property_readonly("model", &BEEngine::getModel);
 }
 
 template <model_type type>
 void wrapModelTypeTrait(py::module& mod) {
   using trait = model_type_traits<type>;
   py::class_<trait>(mod, trait::repr)
       .def_property_readonly_static(
           "dimension", [](const py::object&) { return trait::dimension; },
           "Dimension of computational domain")
       .def_property_readonly_static(
           "components", [](const py::object&) { return trait::components; },
           "Number of components of vector fields")
       .def_property_readonly_static(
           "boundary_dimension",
           [](const py::object&) { return trait::boundary_dimension; },
           "Dimension of boundary of computational domain")
       .def_property_readonly_static(
           "voigt", [](const py::object&) { return trait::voigt; },
           "Number of components of symmetrical tensor fields")
       .def_property_readonly_static(
           "indices", [](const py::object&) { return trait::indices; });
 }
 
 /// Wrap Models
 void wrapModelClass(py::module& mod) {
   py::enum_<model_type>(mod, "model_type")
       .value("basic_1d", model_type::basic_1d,
              "Normal contact with 1D interface")
       .value("basic_2d", model_type::basic_2d,
              "Normal contact with 2D interface")
       .value("surface_1d", model_type::surface_1d,
              "Normal & tangential contact with 1D interface")
       .value("surface_2d", model_type::surface_2d,
              "Normal & tangential contact with 2D interface")
       .value("volume_1d", model_type::volume_1d,
              "Contact with volumetric representation and 1D interface")
       .value("volume_2d", model_type::volume_2d,
              "Contact with volumetric representation and 2D interface");
 
   auto trait_mod = mod.def_submodule("_type_traits");
   trait_mod.doc() = "Defines type trait objects.";
   wrapModelTypeTrait<model_type::basic_1d>(trait_mod);
   wrapModelTypeTrait<model_type::surface_1d>(trait_mod);
   wrapModelTypeTrait<model_type::volume_1d>(trait_mod);
   wrapModelTypeTrait<model_type::basic_2d>(trait_mod);
   wrapModelTypeTrait<model_type::surface_2d>(trait_mod);
   wrapModelTypeTrait<model_type::volume_2d>(trait_mod);
 
   py::class_<model_operator_accessor>(mod, "_model_operator_acessor")
       .def(py::init<Model&>())
       .def(
           "__getitem__",
           [](model_operator_accessor& acc, std::string name) {
             try {
               return acc.get(name);
             } catch (std::out_of_range&) {
               throw py::key_error(name);
             }
           },
           py::return_value_policy::reference_internal)
       .def("__contains__",
            [](model_operator_accessor& acc, std::string key) {
              const auto ops = acc.m.getIntegralOperators();
              return std::find(ops.begin(), ops.end(), key) != ops.end();
            })
       .def(
           "__iter__",
           [](const model_operator_accessor& acc) {
             const auto& ops = acc.m.getIntegralOperatorsMap();
             return py::make_key_iterator(ops.cbegin(), ops.cend());
           },
           py::keep_alive<0, 1>());
 
   py::class_<Model, FieldContainer>(mod, "Model")
       .def(py::init(py::overload_cast<model_type, const std::vector<Real>&,
                                       const std::vector<UInt>&>(
                &ModelFactory::createModel)),
            R"-(Create a new model of a given type, physical size and *global*
 discretization.
 
 :param model_type: the type of desired model
 :param system_size: the physical size of the domain in each direction
 :param global_discretization: number of points in each direction)-")
       .def_property_readonly("type", &Model::getType)
       .def_property("E", &Model::getYoungModulus, &Model::setYoungModulus,
                     "Young's modulus")
       .def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio,
                     "Poisson's ratio")
       .def_property_readonly("mu", &Model::getShearModulus, "Shear modulus")
       .def_property_readonly("E_star", &Model::getHertzModulus,
                              "Contact (Hertz) modulus")
       .def_property_readonly("be_engine", &Model::getBEEngine,
                              "Boundary element engine")
 
       .def(
           "setElasticity",
           [](Model& m, Real E, Real nu) {
             TAMAAS_DEPRECATE("setElasticity()", "the E and nu properties");
             m.setElasticity(E, nu);
           },
           "E"_a, "nu"_a)
       .def(TAMAAS_DEPRECATE_ACCESSOR(getHertzModulus, Model, "E_star"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getYoungModulus, Model, "E"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getShearModulus, Model, "mu"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getPoissonRatio, Model, "nu"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getTraction, Model, "traction"),
            py::return_value_policy::reference_internal)
       .def(TAMAAS_DEPRECATE_ACCESSOR(getDisplacement, Model, "displacement"),
            py::return_value_policy::reference_internal)
       .def(TAMAAS_DEPRECATE_ACCESSOR(getSystemSize, Model, "system_size"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getDiscretization, Model, "shape"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getBoundarySystemSize, Model,
                                      "boundary_system_size"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getBoundaryDiscretization, Model,
                                      "boundary_shape"))
 
       .def("setIntegrationMethod", &Model::setIntegrationMethod)
       .def("solveNeumann", &Model::solveNeumann,
            "Solve surface tractions -> displacements")
       .def("solveDirichlet", &Model::solveDirichlet,
            "Solve surface displacemnts -> tractions")
       .def("dump", &Model::dump, "Write model data to registered dumpers")
       .def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>(),
            "Register a dumper")
 
       .def(
           "getBEEngine",
           [](Model& m) -> decltype(m.getBEEngine()) {
             TAMAAS_DEPRECATE("getBEEngine()", "the be_engine property");
             return m.getBEEngine();
           },
           py::return_value_policy::reference_internal)
       .def(
           "getIntegralOperator",
           [](const Model& m, std::string name) {
             TAMAAS_DEPRECATE("getIntegralOperator()", "the operators property");
             return m.getIntegralOperator(name);
           },
           "operator_name"_a, py::return_value_policy::reference_internal)
       .def(
           "applyElasticity",
           [](Model& model, numpy<Real> stress, numpy<Real> strain) {
             auto out = instanciateFromNumpy(stress);
             auto in = instanciateFromNumpy(strain);
 
             model.applyElasticity(*out, *in);
           },
           "Apply Hooke's law")
 
       // Python magic functions
       .def("__repr__",
            [](const Model& m) {
              std::stringstream ss;
              ss << m;
              return ss.str();
            })
       .def(
           "__copy__",
           [](const Model&) {
             throw std::runtime_error("__copy__ not implemented");
           },
           "Shallow copy of model. Not implemented.")
       .def(
           "__deepcopy__",
           [](const Model& m, py::dict) { return ModelFactory::createModel(m); },
           "Deep copy of model.")
       .def_property_readonly("boundary_fields", &Model::getBoundaryFields)
       .def_property_readonly(
           "operators", [](Model& m) { return model_operator_accessor{m}; },
           "Returns a dict-like object allowing access to the model's "
           "integral operators")
       // More python-like access to model properties
       .def_property_readonly("shape", &Model::getDiscretization,
                              "Discretization (local in MPI environment)")
       .def_property_readonly("global_shape", &Model::getGlobalDiscretization,
                              "Global discretization (in MPI environement)")
       .def_property_readonly("boundary_shape",
                              &Model::getBoundaryDiscretization,
                              "Number of points on boundary")
       .def_property_readonly("system_size", &Model::getSystemSize,
                              "Size of physical domain")
       .def_property_readonly("boundary_system_size",
                              &Model::getBoundarySystemSize,
                              "Physical size of surface")
       .def_property_readonly("traction",
                              (const GridBase<Real>& (Model::*)() const) &
                                  Model::getTraction,
                              "Surface traction field")
       .def_property_readonly("displacement",
                              (const GridBase<Real>& (Model::*)() const) &
                                  Model::getDisplacement,
                              "Displacement field");
 
   py::class_<ModelDumper, PyModelDumper, std::shared_ptr<ModelDumper>>(
       mod, "ModelDumper")
       .def(py::init<>())
       .def("dump", &ModelDumper::dump, "model"_a, "Dump model")
       .def(
           "__lshift__",
           [](ModelDumper& dumper, Model& model) { dumper << model; },
           "Dump model");
 }
 
 /// Wrap factory for models
 void wrapModelFactory(py::module& mod) {
   py::class_<ModelFactory>(mod, "ModelFactory")
       .def_static(
           "createModel",
           py::overload_cast<model_type, const std::vector<Real>&,
                             const std::vector<UInt>&>(
               &ModelFactory::createModel),
           "model_type"_a, "system_size"_a, "global_discretization"_a,
           R"-(Create a new model of a given type, physical size and *global*
 discretization.
 
 :param model_type: the type of desired model
 :param system_size: the physical size of the domain in each direction
 :param global_discretization: number of points in each direction)-")
       .def_static("createModel",
                   py::overload_cast<const Model&>(&ModelFactory::createModel),
                   "model"_a, "Create a deep copy of a model.")
       .def_static(
           "createResidual",
           [](Model& model, Real sigma_y, Real hardening) {
             TAMAAS_DEPRECATE("createResidual", "Residual constructor");
             return ModelFactory::createResidual(model, sigma_y, hardening);
           },
           "model"_a, "sigma_y"_a, "hardening"_a = 0.,
           R"-(Create an isotropic linear hardening residual.
 :param model: the model on which to define the residual
 :param sigma_y: the (von Mises) yield stress
 :param hardening: the hardening modulus)-")
       .def_static("registerVolumeOperators",
                   &ModelFactory::registerVolumeOperators, "model"_a,
                   "Register Boussinesq and Mindlin operators to model.")
       .def_static("registerNonPeriodic", &ModelFactory::registerNonPeriodic,
                   "model"_a, "name"_a,
                   "Register non-periodic Boussinesq operator")
       .def_static("setIntegrationMethod", &ModelFactory::setIntegrationMethod,
                   "operator"_a, "method"_a, "cutoff"_a,
                   "Set the integration method (linear or cutoff) for a volume "
                   "integral operator")
       .def_static(
           "registerHookeField",
           [](Model& m, const std::string& name) {
             model_type_dispatch(
                 [&name, &m](auto&& type) {
                   constexpr auto mtype = std::decay_t<decltype(type)>::value;
                   return m.registerIntegralOperator<HookeField<mtype>>(name);
                 },
                 m.getType());
           },
           "model"_a, "name"_a, "Register a HookeField operator");
 }
 
 /// Wrap residual class
 void wrapResidual(py::module& mod) {
   py::class_<Residual, PyResidual>(mod, "Residual")
       .def(py::init<Model&, std::shared_ptr<Material>>(),
            py::keep_alive<1, 2>(), "model"_a, "material"_a,
            R"-(Create a residual object with a material.
 Defines the following residual equation:
     ɛ - ∇N[σ(ɛ)] - ∇M[t] = 0
 
 Where σ(ɛ) is the *eigenstress* associated with the constitutive law.
 
 :param model: the model on which to define the residual
 :param material: material object which defines a constitutive law
 )-")
       .def("computeResidual", &Residual::computeResidual)
       .def("updateState", &Residual::updateState)
       .def("computeResidualDisplacement",
            &Residual::computeResidualDisplacement)
       .def(TAMAAS_DEPRECATE_ACCESSOR(getVector, Residual, "vector"),
            py::return_value_policy::reference_internal)
       .def(
           "getStress",
           [](const Residual& res) {
             TAMAAS_DEPRECATE("getStress", "model[\"stress\"]");
             const auto& stress = res.getModel().field<Real>("stress");
             Grid<Real, 3> view(res.getModel().getDiscretization(),
                                stress.getNbComponents(), stress.view());
             return view;
           },
           py::return_value_policy::reference_internal)
       .def("setIntegrationMethod",
            [](Residual& residual, integration_method method, Real cutoff) {
              TAMAAS_DEPRECATE("setIntegrationMethod",
                               "model.setIntegrationMethod()");
              residual.getModel().setIntegrationMethod(method, cutoff);
            })
       .def_property_readonly("vector", &Residual::getVector)
       .def_property_readonly(
           "model",
           static_cast<const Model& (Residual::*)() const>(&Residual::getModel))
       .def_property_readonly("material", &Residual::getMaterial);
 }
 
 void wrapModel(py::module& mod) {
   wrapFieldContainer(mod);
   wrapBEEngine(mod);
   wrapModelClass(mod);
   wrapModelFactory(mod);
   wrapFunctionals(mod);
   wrapResidual(mod);
   wrapIntegralOperator(mod);
 }
 
 }  // namespace wrap
 
 }  // namespace tamaas
diff --git a/python/wrap/model_extensions.hh b/python/wrap/model_extensions.hh
index d82b65a..50ac281 100644
--- a/python/wrap/model_extensions.hh
+++ b/python/wrap/model_extensions.hh
@@ -1,95 +1,104 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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 "functional.hh"
 #include "model.hh"
 #include "model_dumper.hh"
 #include "residual.hh"
 #include "wrap.hh"
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 namespace functional {
 namespace wrap {
 /// Class for extension of Functional in python
 class PyFunctional : public Functional {
 public:
   using Functional::Functional;
 
   // Overriding pure virtual functions
   Real computeF(GridBase<Real>& variable, GridBase<Real>& dual) const override {
     // NOLINTNEXTLINE(readability-else-after-return)
     PYBIND11_OVERLOAD_PURE(Real, Functional, computeF, variable, dual);
   }
 
   void computeGradF(GridBase<Real>& variable,
                     GridBase<Real>& gradient) const override {
     // NOLINTNEXTLINE(readability-else-after-return)
     PYBIND11_OVERLOAD_PURE(void, Functional, computeGradF, variable, gradient);
   }
 };
 }  // namespace wrap
 }  // namespace functional
 /* -------------------------------------------------------------------------- */
 namespace wrap {
 
 /* -------------------------------------------------------------------------- */
 class PyModelDumper : public ModelDumper {
 public:
   using ModelDumper::ModelDumper;
 
   void dump(const Model& model) override {
     // NOLINTNEXTLINE(readability-else-after-return)
     PYBIND11_OVERLOAD_PURE(void, ModelDumper, dump, model);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 class PyResidual : public Residual {
 public:
   using Residual::Residual;
 
   void computeResidual(GridBase<Real>& strain_increment) override {
     // NOLINTNEXTLINE(readability-else-after-return)
     PYBIND11_OVERLOAD(void, Residual, computeResidual, strain_increment);
   }
 
   void updateState(GridBase<Real>& converged_strain_increment) override {
     // NOLINTNEXTLINE(readability-else-after-return)
     PYBIND11_OVERLOAD(void, Residual, updateState,
                            converged_strain_increment);
   }
 
   void computeResidualDisplacement(GridBase<Real>& strain_increment) override {
     // NOLINTNEXTLINE(readability-else-after-return)
     PYBIND11_OVERLOAD(void, Residual, computeResidualDisplacement,
                            strain_increment);
   }
 };
 
+class PyIntegralOperator : public IntegralOperator {
+public:
+  using IntegralOperator::IntegralOperator;
+
+  void apply(GridBase<Real>& input, GridBase<Real>& output) const override {
+    PYBIND11_OVERLOAD_PURE(void, IntegralOperator, apply, input, output);
+  }
+};
+
 /* -------------------------------------------------------------------------- */
 
 }  // namespace wrap
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/model/be_engine.cpp b/src/model/be_engine.cpp
index 81f5921..5e5054a 100644
--- a/src/model/be_engine.cpp
+++ b/src/model/be_engine.cpp
@@ -1,89 +1,91 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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 "be_engine.hh"
 #include "logger.hh"
 #include "model.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 void solve(IntegralOperator::kind kind,
-           const std::map<IntegralOperator::kind, IntegralOperator*>& operators,
+           const std::map<IntegralOperator::kind,
+                          std::shared_ptr<IntegralOperator>>& operators,
            GridBase<Real>& input, GridBase<Real>& output) try {
   operators.at(kind)->apply(input, output);
 } catch (std::out_of_range& e) {
   Logger().get(LogLevel::warning)
       << "Operator ("
       << "Westergaard::" << kind << ") not registered\n";
   throw e;
 }
 
 template <model_type mtype, IntegralOperator::kind otype>
 void registerWestergaardOperator(
-    std::map<IntegralOperator::kind, IntegralOperator*>& operators,
+    std::map<IntegralOperator::kind, std::shared_ptr<IntegralOperator>>&
+        operators,
     Model& model) {
   std::stringstream sstr;
   sstr << "Westergaard::" << otype;
   if (operators.find(otype) == operators.end())
     operators[otype] =
         model.template registerIntegralOperator<Westergaard<mtype, otype>>(
             sstr.str());
 }
 
 template <model_type type>
 void BEEngineTmpl<type>::solveNeumann(GridBase<Real>& neumann,
                                       GridBase<Real>& dirichlet) const {
   solve(IntegralOperator::neumann, this->operators, neumann, dirichlet);
 }
 
 template <model_type type>
 void BEEngineTmpl<type>::solveDirichlet(GridBase<Real>& dirichlet,
                                         GridBase<Real>& neumann) const {
   solve(IntegralOperator::dirichlet, this->operators, dirichlet, neumann);
 }
 
 template <model_type type>
 void BEEngineTmpl<type>::registerNeumann() {
   Logger().get(LogLevel::debug)
       << TAMAAS_DEBUG_MSG("Registering Neumann Westergaard");
   registerWestergaardOperator<type, IntegralOperator::neumann>(this->operators,
                                                                *this->model);
 }
 /// Register dirichlet operator
 template <model_type type>
 void BEEngineTmpl<type>::registerDirichlet() {
   Logger().get(LogLevel::debug)
       << TAMAAS_DEBUG_MSG("Registering Dirichlet Westergaard");
   registerWestergaardOperator<type, IntegralOperator::dirichlet>(
       this->operators, *this->model);
 }
 
 template class BEEngineTmpl<model_type::basic_1d>;
 template class BEEngineTmpl<model_type::basic_2d>;
 template class BEEngineTmpl<model_type::surface_1d>;
 template class BEEngineTmpl<model_type::surface_2d>;
 template class BEEngineTmpl<model_type::volume_1d>;
 template class BEEngineTmpl<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/model/be_engine.hh b/src/model/be_engine.hh
index f13be4a..4c94907 100644
--- a/src/model/be_engine.hh
+++ b/src/model/be_engine.hh
@@ -1,87 +1,86 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef BE_ENGINE_HH
 #define BE_ENGINE_HH
 /* -------------------------------------------------------------------------- */
-#include "model_type.hh"
 #include "integral_operator.hh"
+#include "model_type.hh"
 #include "westergaard.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /**
  * @brief Boundary equation engine class. Solves the Neumann/Dirichlet problem
  * This class should be dimension and model-type agnostic.
  */
 class BEEngine {
 public:
   BEEngine(Model* model) : model(model) {}
   /// Destructor
   virtual ~BEEngine() = default;
 
 public:
   /// Solve Neumann problem (expects boundary data)
   virtual void solveNeumann(GridBase<Real>& neumann,
                             GridBase<Real>& dirichlet) const = 0;
   /// Solve Dirichlet problem (expects boundary data)
   virtual void solveDirichlet(GridBase<Real>& dirichlet,
                               GridBase<Real>& neumann) const = 0;
 
   /// Register neumann operator
   virtual void registerNeumann() = 0;
   /// Register dirichlet operator
   virtual void registerDirichlet() = 0;
 
   /// Get model
   const Model& getModel() const { return *model; }
 
   /// Compute L_2 norm of influence functions
   Real getNeumannNorm() {
     return operators[IntegralOperator::neumann]->getOperatorNorm();
   }
 
 protected:
   Model* model;
-  std::map<IntegralOperator::kind, IntegralOperator*> operators;
+  std::map<IntegralOperator::kind, std::shared_ptr<IntegralOperator>> operators;
 };
 
 template <model_type type>
 class BEEngineTmpl : public BEEngine {
 public:
-  BEEngineTmpl(Model* model): BEEngine(model) {}
+  BEEngineTmpl(Model* model) : BEEngine(model) {}
 
   void solveNeumann(GridBase<Real>& neumann,
-		    GridBase<Real>& dirichlet) const override;
+                    GridBase<Real>& dirichlet) const override;
   void solveDirichlet(GridBase<Real>& dirichlet,
                       GridBase<Real>& neumann) const override;
   /// Register neumann operator
   void registerNeumann() override;
   /// Register dirichlet operator
   void registerDirichlet() override;
-
 };
 
 }  // namespace tamaas
 
 #endif  // BE_ENGINE_HH
diff --git a/src/model/model.cpp b/src/model/model.cpp
index 1941f94..a9a87d1 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,168 +1,175 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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 "be_engine.hh"
 #include "logger.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 void Model::setElasticity(Real E_, Real nu_) {
   setYoungModulus(E_);
   setPoissonRatio(nu_);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::applyElasticity(GridBase<Real>& stress,
                             const GridBase<Real>& strain) const {
   operators.at("hooke")->apply(const_cast<GridBase<Real>&>(strain), stress);
 }
 
 /* -------------------------------------------------------------------------- */
 GridBase<Real>& Model::getTraction() { return this->field<Real>("traction"); }
 const GridBase<Real>& Model::getTraction() const {
   return this->field<Real>("traction");
 }
 
 /* -------------------------------------------------------------------------- */
 GridBase<Real>& Model::getDisplacement() {
   return this->field<Real>("displacement");
 }
 const GridBase<Real>& Model::getDisplacement() const {
   return this->field<Real>("displacement");
 }
 
 /* -------------------------------------------------------------------------- */
 const std::vector<Real>& Model::getSystemSize() const { return system_size; }
 
 /* -------------------------------------------------------------------------- */
 const std::vector<UInt>& Model::getDiscretization() const {
   return discretization;
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::solveNeumann() {
   engine->registerNeumann();
   engine->solveNeumann(getTraction(), getDisplacement());
 }
 
 void Model::solveDirichlet() {
   engine->registerDirichlet();
   engine->solveDirichlet(getDisplacement(), getTraction());
 }
 
 /* -------------------------------------------------------------------------- */
-IntegralOperator* Model::getIntegralOperator(const std::string& name) const {
-  return operators.at(name).get();
+void Model::registerIntegralOperator(const std::string& name,
+                                     std::shared_ptr<IntegralOperator> op) {
+  operators[name] = std::move(op);
+}
+
+/* -------------------------------------------------------------------------- */
+std::shared_ptr<IntegralOperator>
+Model::getIntegralOperator(const std::string& name) const {
+  return operators.at(name);
 }
 
 /* -------------------------------------------------------------------------- */
 std::vector<std::string> Model::getIntegralOperators() const {
   std::vector<std::string> keys;
   keys.reserve(operators.size());
 
   std::transform(operators.begin(), operators.end(), std::back_inserter(keys),
                  [](auto&& pair) { return pair.first; });
   return keys;
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::updateOperators() {
   for (auto& op : operators)
     op.second->updateFromModel();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumper(std::shared_ptr<ModelDumper> dumper) {
   this->dumpers.push_back(std::move(dumper));
 }
 
 void Model::dump() const {
   if (dumpers.empty())
     Logger().get(LogLevel::warning) << "dumper list is empty\n";
 
   for (const auto& dumper : dumpers)
     if (dumper)
       dumper->dump(*this);
 }
 
 /* -------------------------------------------------------------------------- */
 std::vector<std::string> Model::getBoundaryFields() const {
   std::vector<std::string> boundary_fields;
 
   for (auto&& pair : this->fields_map()) {
     const auto boundary = boost::apply_visitor(
         [&](auto&& grid_ptr) { return this->isBoundaryField(*grid_ptr); },
         pair.second);
 
     if (boundary)
       boundary_fields.push_back(pair.first);
   }
 
   return boundary_fields;
 }
 
 /* -------------------------------------------------------------------------- */
 std::ostream& operator<<(std::ostream& o, const Model& _this) {
   o << "Model<" << _this.getType() << "> (E = " << _this.getYoungModulus()
     << ", nu = " << _this.getPoissonRatio() << ")\n";
 
   auto out_collec = [&o](auto&& collec) {
     std::for_each(collec.begin(), collec.end() - 1,
                   [&o](const auto& x) { o << x << ", "; });
     o << collec.back();
   };
 
   // Printing domain size
   o << "  - domain = [";
   out_collec(_this.getSystemSize());
   o << "]\n";
 
   // Printing discretization
   o << "  - discretization = [";
   out_collec(_this.getDiscretization());
   o << "]\n";
 
   if (mpi::size() > 1) {
     o << "  - global discretization = [";
     out_collec(_this.getGlobalDiscretization());
     o << "]\n";
   }
 
   // Print fields
   o << "  - registered fields = [";
   out_collec(_this.fields());
   o << "]\n";
 
   o << "  - registered operators = [";
   out_collec(_this.getIntegralOperators());
   o << "]";
 
   if (not _this.dumpers.empty())
     o << "\n  - " << _this.dumpers.size() << " registered dumpers";
 
   return o;
 }
 
 /* --------------------------------------------------------------------------
  */
 }  // namespace tamaas
diff --git a/src/model/model.hh b/src/model/model.hh
index 6da6a7e..168129c 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,202 +1,208 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef MODEL_HH
 #define MODEL_HH
 /* -------------------------------------------------------------------------- */
 #include "be_engine.hh"
 #include "field_container.hh"
 #include "grid_base.hh"
 #include "integral_operator.hh"
 #include "model_dumper.hh"
 #include "model_type.hh"
 #include "tamaas.hh"
 #include <algorithm>
 #include <memory>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 
 /// Forward declaration
 enum class integration_method;
 
 /**
  * @brief Model containing pressure and displacement
  * This class is a container for the model fields. It is supposed to be
  * dimension agnostic, hence the GridBase members.
  */
 class Model : public FieldContainer {
 protected:
   /// Constructor
   Model(std::vector<Real> system_size, std::vector<UInt> discretization)
       : FieldContainer(), system_size(std::move(system_size)),
         discretization(std::move(discretization)) {}
 
 public:
   /// Set elasticity parameters
   void setElasticity(Real E, Real nu);
   /// Get Hertz contact modulus
   Real getHertzModulus() const { return E / (1 - nu * nu); }
   /// Get Young's modulus
   Real getYoungModulus() const { return E; }
   /// Get Poisson's ratio
   Real getPoissonRatio() const { return nu; }
   /// Get shear modulus
   Real getShearModulus() const { return E / (2 * (1 + nu)); }
   /// Set Young's modulus
   void setYoungModulus(Real E_) {
     if (E_ < 0)
       TAMAAS_EXCEPTION("Elastic modulus should be positive");
     this->E = E_;
     updateOperators();
   }
 
   /// Set Poisson's ratio
   void setPoissonRatio(Real nu_) {
     if (nu_ > 0.5 or nu_ <= -1)
       TAMAAS_EXCEPTION("Poisson ratio should be in ]-1, 0.5]");
     this->nu = nu_;
     updateOperators();
   }
 
   // Model info accessors
 public:
   /// Get model type
   virtual model_type getType() const = 0;
   /// Get system physical size
   const std::vector<Real>& getSystemSize() const;
   /// Get boundary system physical size
   virtual std::vector<Real> getBoundarySystemSize() const = 0;
   /// Get discretization
   const std::vector<UInt>& getDiscretization() const;
   /// Get discretization of global MPI system
   virtual std::vector<UInt> getGlobalDiscretization() const = 0;
   /// Get boundary discretization
   virtual std::vector<UInt> getBoundaryDiscretization() const = 0;
   /// Get boundary element engine
   BEEngine& getBEEngine() {
     TAMAAS_ASSERT(engine, "BEEngine was not initialized");
     return *engine;
   }
 
   // Exposing some common operators
 public:
   /// Apply Hooke's law
   void applyElasticity(GridBase<Real>& stress,
                        const GridBase<Real>& strain) const;
 
   /// Solve Neumann problem using default neumann operator
   void solveNeumann();
   /// Solve Dirichlet problem using default dirichlet operator
   void solveDirichlet();
 
 public:
   /// Register a new integral operator
   template <typename Operator>
-  IntegralOperator* registerIntegralOperator(const std::string& name);
+  std::shared_ptr<IntegralOperator>
+  registerIntegralOperator(const std::string& name);
+  /// Register external operator
+  void registerIntegralOperator(const std::string& name,
+                                std::shared_ptr<IntegralOperator> op);
   /// Get a registerd integral operator
-  IntegralOperator* getIntegralOperator(const std::string& name) const;
+  std::shared_ptr<IntegralOperator>
+  getIntegralOperator(const std::string& name) const;
   /// Get list of integral operators
   std::vector<std::string> getIntegralOperators() const;
   /// Get operators mapcar
   const auto& getIntegralOperatorsMap() const { return operators; }
   /// Tell operators to update their cache
   void updateOperators();
   /// Set integration method for registered volume operators
   virtual void setIntegrationMethod(integration_method method, Real cutoff) = 0;
 
 public:
   /// Get pressure
   GridBase<Real>& getTraction();
   /// Get pressure
   const GridBase<Real>& getTraction() const;
   /// Get displacement
   GridBase<Real>& getDisplacement();
   /// Get displacement
   const GridBase<Real>& getDisplacement() const;
   /// Determine if a field is defined on boundary
   template <class T>
   bool isBoundaryField(const GridBase<T>& field) const {
     return field.getDimension() == getBoundaryDiscretization().size();
   }
   /// Return list of fields defined on boundary
   std::vector<std::string> getBoundaryFields() const;
 
 public:
   /// Set the dumper object
   void addDumper(std::shared_ptr<ModelDumper> dumper);
   /// Dump the model
   void dump() const;
 
   friend std::ostream& operator<<(std::ostream& o, const Model& _this);
 
 protected:
   Real E = 1, nu = 0;
 
   std::vector<Real> system_size;
   std::vector<UInt> discretization;
   std::unique_ptr<BEEngine> engine = nullptr;
   std::unordered_map<std::string, std::shared_ptr<IntegralOperator>> operators;
   std::vector<std::shared_ptr<ModelDumper>> dumpers;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template functions */
 /* -------------------------------------------------------------------------- */
 template <typename Operator>
-IntegralOperator* Model::registerIntegralOperator(const std::string& name) {
+std::shared_ptr<IntegralOperator>
+Model::registerIntegralOperator(const std::string& name) {
   Logger().get(LogLevel::debug)
       << TAMAAS_DEBUG_MSG("registering operator " + name);
   operators[name] = std::make_unique<Operator>(this);
-  return operators[name].get();
+  return operators[name];
 }
 
 /* -------------------------------------------------------------------------- */
 /* Output model to stream */
 /* -------------------------------------------------------------------------- */
 std::ostream& operator<<(std::ostream& o, const Model& _this);
 
 /* -------------------------------------------------------------------------- */
 /* Simpler grid allocation */
 /* -------------------------------------------------------------------------- */
 template <bool boundary, typename T>
 std::unique_ptr<GridBase<T>> allocateGrid(Model& model) {
   return allocateGrid<boundary, T>(
       model.getType(), (boundary) ? model.getBoundaryDiscretization()
                                   : model.getDiscretization());
 }
 
 template <bool boundary, typename T>
 std::unique_ptr<GridBase<T>> allocateGrid(Model& model, UInt nc) {
   return allocateGrid<boundary, T>(model.getType(),
                                    (boundary)
                                        ? model.getBoundaryDiscretization()
                                        : model.getDiscretization(),
                                    nc);
 }
 
 }  // namespace tamaas
 
 #endif  // MODEL_HH
diff --git a/src/solvers/polonsky_keer_rey.hh b/src/solvers/polonsky_keer_rey.hh
index e712df3..892ff7b 100644
--- a/src/solvers/polonsky_keer_rey.hh
+++ b/src/solvers/polonsky_keer_rey.hh
@@ -1,137 +1,137 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef POLONSKY_KEER_REY_HH
 #define POLONSKY_KEER_REY_HH
 /* -------------------------------------------------------------------------- */
 #include "contact_solver.hh"
 #include "grid_view.hh"
 #include "meta_functional.hh"
 #include "westergaard.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 class PolonskyKeerRey : public ContactSolver {
 public:
   /// Types of algorithm (primal/dual) or constraint
   enum type { gap, pressure };
 
 public:
   /// Constructor
   PolonskyKeerRey(Model& model, const GridBase<Real>& surface, Real tolerance,
                   type variable_type, type constraint_type);
   ~PolonskyKeerRey() override = default;
 
 public:
   /// \cond DO_NOT_DOCUMENT
   using ContactSolver::solve;
   ///\endcond
 
   /// Solve
   Real solve(std::vector<Real> target) override;
   /// Mean on unsaturated constraint zone
   virtual Real meanOnUnsaturated(const GridBase<Real>& field) const;
   /// Compute squared norm
   virtual Real computeSquaredNorm(const GridBase<Real>& field) const;
   /// Update search direction
   virtual void updateSearchDirection(Real factor);
   /// Compute critical step
   virtual Real computeCriticalStep(Real target = 0);
   /// Update primal and check non-admissible state
   virtual bool updatePrimal(Real step);
   /// Compute error/stopping criterion
   virtual Real computeError() const;
   /// Enforce contact constraints on final state
   virtual void enforceAdmissibleState();
   /// Enfore mean value constraint
   virtual void enforceMeanValue(Real mean);
 
   /// Set integral operator for gradient computation
   void setIntegralOperator(std::string name);
 
 private:
   /// Set correct views on normal traction and gap
   template <model_type type>
   void setViews();
 
 protected:
   type variable_type, constraint_type;
   model_type operation_type;
   GridBase<Real>* primal = nullptr;  ///< non-owning pointer for primal varaible
   GridBase<Real>* dual = nullptr;    ///< non-owning pointer for dual variable
 
   /// CG search direction
   std::unique_ptr<GridBase<Real>> search_direction = nullptr;
   /// Projected CG search direction
   std::unique_ptr<GridBase<Real>> projected_search_direction = nullptr;
   /// View on normal pressure, gap, displacement
   std::unique_ptr<GridBase<Real>> pressure_view, gap_view,
       displacement_view = nullptr;
 
   /// Integral operator for gradient computation
-  IntegralOperator* integral_op = nullptr;
+  std::shared_ptr<IntegralOperator> integral_op = nullptr;
 };
 
 /* -------------------------------------------------------------------------- */
 template <model_type mtype>
 void PolonskyKeerRey::setViews() {
   constexpr UInt dim = model_type_traits<mtype>::dimension;
   constexpr UInt bdim = model_type_traits<mtype>::boundary_dimension;
   constexpr UInt comp = model_type_traits<mtype>::components;
 
   pressure_view = std::unique_ptr<GridBase<Real>>{
       new GridView<Grid, Real, bdim, bdim>(model.getTraction(), {}, comp - 1)};
   gap_view = std::unique_ptr<GridBase<Real>>{
       new GridView<Grid, Real, bdim, bdim>(*this->_gap, {}, comp - 1)};
   displacement_view =
       std::unique_ptr<GridBase<Real>>{new GridView<Grid, Real, dim, bdim>(
           model.getDisplacement(), model_type_traits<mtype>::indices,
           comp - 1)};
 
 #define WESTERGAARD(type, kind, desc)                                          \
   this->model.registerIntegralOperator<                                        \
       Westergaard<model_type::type, IntegralOperator::kind>>(desc)
 
-  IntegralOperator *dirichlet, *neumann;
+  std::shared_ptr<IntegralOperator> dirichlet, neumann;
 
   // I noz is fugly - TODO factory?
   if (bdim == 1) {
     operation_type = model_type::basic_1d;
     dirichlet = WESTERGAARD(basic_1d, dirichlet, "westergaard_dirichlet");
     neumann = WESTERGAARD(basic_1d, neumann, "westergaard_neumann");
   } else if (bdim == 2) {
     operation_type = model_type::basic_2d;
     dirichlet = WESTERGAARD(basic_2d, dirichlet, "westergaard_dirichlet");
     neumann = WESTERGAARD(basic_2d, neumann, "westergaard_neumann");
   }
 
   if (variable_type == gap)
-    integral_op = dirichlet;
+    integral_op = std::move(dirichlet);
   else
-    integral_op = neumann;
+    integral_op = std::move(neumann);
 
 #undef WESTERGAARD
 }
 
 }  // namespace tamaas
 
 #endif  // POLONSKY_KEER_REY_HH