diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3306b7a6f02afc544e4ce2b28ba3f1da41d689f5..0b55368662b1a287089ed5393a119b9205954567 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -28,7 +28,7 @@ IF(NOT CMAKE_BUILD_TYPE)
          "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel."
          FORCE)
 ENDIF(NOT CMAKE_BUILD_TYPE)
-MESSAGE(STATUS "Build type: ${CMAKE_BUILD_TYPE}")
+MESSAGE(STATUS "CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE}")
 
 LIST(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/CMake")
 
diff --git a/flow/_src/floww.i b/flow/_src/floww.i
index 9a6b59d52066dfb1e5bb2c51a614f45dc70e2892..9bfe0570846b07bfaac4c9bcda9aef3e33adcb89 100644
--- a/flow/_src/floww.i
+++ b/flow/_src/floww.i
@@ -19,16 +19,13 @@
 %feature("autodoc","1");
 
 %module(docstring=
-"'floww' module: full potential solver 2018-2019
+"'floww' module: full potential solver 2018-2020
 (c) ULg - A&M",
-directors="1",
+directors="0",
 threads="1"
 ) floww
 %{
 
-#include <string>
-#include <sstream>
-#include <typeinfo>
 #include "flow.h"
 
 #include "fwkw.h"
@@ -74,29 +71,21 @@ threads="1"
 %include "wWakeElement.h"
 %include "wKuttaElement.h"
 
-%feature("director") flow::F0PsC;
-%pythonappend flow::F0PsC::F0PsC "self.__disown__()"    // only for directors
-%feature("director") flow::F0PsPhiInf;
-%pythonappend flow::F0PsPhiInf::F0PsPhiInf "self.__disown__()"    // only for directors
-%nodefault flow::F0Ps;     // is pure virtual but swig doesn't know it.
+%shared_ptr(flow::F0Ps);
+%shared_ptr(flow::F0PsC);
+%shared_ptr(flow::F0PsPhiInf);
+%nodefault flow::F0Ps; // is pure virtual but swig doesn't know it.
 %include "wF0Ps.h"
 
-%feature("director") flow::F0ElRho;
-%pythonappend flow::F0ElRho::F0ElRho "self.__disown__()"    // only for directors
-%feature("director") flow::F0ElRhoL;
-%pythonappend flow::F0ElRhoL::F0ElRhoL "self.__disown__()"    // only for directors
-%feature("director") flow::F0ElMach;
-%pythonappend flow::F0ElMach::F0ElMach "self.__disown__()"    // only for directors
-%feature("director") flow::F0ElMachL;
-%pythonappend flow::F0ElMachL::F0ElMachL "self.__disown__()"    // only for directors
-%feature("director") flow::F0ElCp;
-%pythonappend flow::F0ElCp::F0ElCp "self.__disown__()"    // only for directors
-%feature("director") flow::F0ElCpL;
-%pythonappend flow::F0ElCpL::F0ElCpL "self.__disown__()"    // only for directors
+%shared_ptr(flow::F0ElRho);
+%shared_ptr(flow::F0ElRhoL);
+%shared_ptr(flow::F0ElMach);
+%shared_ptr(flow::F0ElMachL);
+%shared_ptr(flow::F0ElCp);
+%shared_ptr(flow::F0ElCpL);
 %include "wF0El.h"
 
-%feature("director") flow::F1ElVi;
-%pythonappend flow::F1ElVi::F1ElVi "self.__disown__()"    // only for directors
+%shared_ptr(flow::F1ElVi);
 %include "wF1El.h"
 
 %immutable flow::Problem::msh; // avoids the creation of the setter method
diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index cb04c1c55541f2051a44f9d6a4960e85def2d9e7..a31bed4133b978636d51012f646f1694cf5947ed 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -135,12 +135,12 @@ void Adjoint::buildDR(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR)
         Element *eU = p.second[0];
 
         // Subsonic contribution (drho*grad_phi*grad_psi + rho*grad_phi*grad_psi)
-        Eigen::MatrixXd Ae1 = e->buildDK(sol->phi, fluid->rho);
-        Eigen::MatrixXd Ae2 = e->buildK(sol->phi, fluid->rho);
+        Eigen::MatrixXd Ae1 = e->buildDK(sol->phi, *fluid->rho);
+        Eigen::MatrixXd Ae2 = e->buildK(sol->phi, *fluid->rho);
         // Supersonic contribution
         double mCOv = 0.95; // solver should be converged, so we set the lowest numerical viscosity
         double muCv = 1.0;
-        double mach = fluid->mach.eval(*e, sol->phi, 0);
+        double mach = fluid->mach->eval(*e, sol->phi, 0);
         if (mach > mCOv)
         {
             // switching function and derivative
@@ -152,13 +152,13 @@ void Adjoint::buildDR(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR)
             Ae2 *= 1 - mu;
 
             // 3rd term (mu*drhoU*grad_phi*grad_psi)
-            Eigen::MatrixXd Ae3 = mu * fluid->rho.evalD(*eU, sol->phi, 0) * e->buildDs(sol->phi, Fct0C(1.)) * eU->computeGrad(sol->phi, 0).transpose() * eU->getVCache().getDsf(0);
+            Eigen::MatrixXd Ae3 = mu * fluid->rho->evalD(*eU, sol->phi, 0) * e->buildDs(sol->phi, Fct0C(1.)) * eU->computeGrad(sol->phi, 0).transpose() * eU->getVCache().getDsf(0);
 
             // 4th term (mu*rhoU*grad_phi*grad_psi)
-            Eigen::MatrixXd Ae4 = mu * fluid->rho.eval(*eU, sol->phi, 0) * e->buildK(sol->phi, Fct0C(1.));
+            Eigen::MatrixXd Ae4 = mu * fluid->rho->eval(*eU, sol->phi, 0) * e->buildK(sol->phi, Fct0C(1.));
 
             // 5th term dmu*(rhoU-rho)*grad_phi*grad_psi
-            Eigen::MatrixXd Ae5 = dmu * (fluid->rho.eval(*eU, sol->phi, 0) - fluid->rho.eval(*e, sol->phi, 0)) * e->buildDK(sol->phi, fluid->mach);
+            Eigen::MatrixXd Ae5 = dmu * (fluid->rho->eval(*eU, sol->phi, 0) - fluid->rho->eval(*e, sol->phi, 0)) * e->buildDK(sol->phi, *fluid->mach);
 
             // Assembly (supersonic)
             tbb::spin_mutex::scoped_lock lock(mutex);
@@ -292,10 +292,10 @@ void Adjoint::buildDI(Eigen::Map<Eigen::VectorXd> &dL, Eigen::Map<Eigen::VectorX
             Eigen::Vector3d Vi(0., 0., 0.);
             Vi(0) = -sin(sol->pbl->alpha);
             Vi(sol->pbl->nDim - 1) = cos(sol->pbl->alpha);
-            double factorL = -1 / sol->pbl->S_ref * sol->pbl->medium->cP.evalD(*eV, sol->phi, 0) * e->normal().dot(Vi);
+            double factorL = -1 / sol->pbl->S_ref * sol->pbl->medium->cP->evalD(*eV, sol->phi, 0) * e->normal().dot(Vi);
             Vi(0) = cos(sol->pbl->alpha);
             Vi(sol->pbl->nDim - 1) = sin(sol->pbl->alpha);
-            double factorD = -1 / sol->pbl->S_ref * sol->pbl->medium->cP.evalD(*eV, sol->phi, 0) * e->normal().dot(Vi);
+            double factorD = -1 / sol->pbl->S_ref * sol->pbl->medium->cP->evalD(*eV, sol->phi, 0) * e->normal().dot(Vi);
 
             // Assembly
             tbb::spin_mutex::scoped_lock lock(mutex);
diff --git a/flow/src/wAssign.cpp b/flow/src/wAssign.cpp
index 919faf2651abac035ba5410c4c86b5ea39188799..f97a4ebe1bf83c7d9b6194adac37b63d4bc5be41 100644
--- a/flow/src/wAssign.cpp
+++ b/flow/src/wAssign.cpp
@@ -23,11 +23,11 @@
 using namespace tbox;
 using namespace flow;
 
-Assign::Assign(std::shared_ptr<MshData> _msh, int no, F0Ps &_f) : Group(_msh, no), f(_f)
+Assign::Assign(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<F0Ps> _f) : Group(_msh, no), f(_f)
 {
     getNodes();
 }
-Assign::Assign(std::shared_ptr<MshData> _msh, std::string const &name, F0Ps &_f) : Group(_msh, name), f(_f)
+Assign::Assign(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<F0Ps> _f) : Group(_msh, name), f(_f)
 {
     getNodes();
 }
@@ -41,7 +41,7 @@ Assign::Assign(std::shared_ptr<MshData> _msh, std::string const &name, F0Ps &_f)
 void Assign::apply(std::vector<double> &vec)
 {
     for (auto n : nodes)
-        vec[n->row] = f.eval(n->pos);
+        vec[n->row] = f->eval(n->pos);
 }
 
 /**
@@ -58,10 +58,10 @@ void Assign::getNodes()
     nodes.resize(std::distance(nodes.begin(), it));
 }
 
-Initial::Initial(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f) : Assign(_msh, no, _f)
+Initial::Initial(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<F0Ps> _f) : Assign(_msh, no, _f)
 {
 }
-Initial::Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f) : Assign(_msh, name, _f)
+Initial::Initial(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<F0Ps> _f) : Assign(_msh, name, _f)
 {
 }
 void Initial::write(std::ostream &out) const
@@ -69,14 +69,14 @@ void Initial::write(std::ostream &out) const
     out << "flow::Initial condition on " << *tag << ")\n";
     for (auto n : nodes)
     {
-        out << '\t' << *n << ": (val=" << f.eval(n->pos) << '\n';
+        out << '\t' << *n << ": (val=" << f->eval(n->pos) << '\n';
     }
 }
 
-Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f) : Assign(_msh, no, _f)
+Dirichlet::Dirichlet(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<F0Ps> _f) : Assign(_msh, no, _f)
 {
 }
-Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f) : Assign(_msh, name, _f)
+Dirichlet::Dirichlet(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<F0Ps> _f) : Assign(_msh, name, _f)
 {
 }
 void Dirichlet::write(std::ostream &out) const
@@ -84,6 +84,6 @@ void Dirichlet::write(std::ostream &out) const
     out << "flow::Dirichlet boundary condition on " << *tag << ")\n";
     for (auto n : nodes)
     {
-        out << '\t' << *n << ": (val=" << f.eval(n->pos) << '\n';
+        out << '\t' << *n << ": (val=" << f->eval(n->pos) << '\n';
     }
 }
\ No newline at end of file
diff --git a/flow/src/wAssign.h b/flow/src/wAssign.h
index 81b8ab0440d96b0627dd8465ff010974000256cd..b6a9c03af240d4f16feb37089a60fb284fa74e57 100644
--- a/flow/src/wAssign.h
+++ b/flow/src/wAssign.h
@@ -36,13 +36,12 @@ class FLOW_API Assign : public tbox::Group
 public:
     std::vector<tbox::Node *> nodes; ///< nodes of the boundary
 #ifndef SWIG
-    F0Ps &f; ///< position-based function to compute the boundary value
+    std::shared_ptr<F0Ps> f; ///< position-based function to compute the boundary value
 #endif
 
-    Assign(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f);
-    Assign(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f);
-
-    ~Assign() { std::cout << "~Assign()\n"; }
+    Assign(std::shared_ptr<tbox::MshData> _msh, int no, std::shared_ptr<F0Ps> _f);
+    Assign(std::shared_ptr<tbox::MshData> _msh, std::string const &name, std::shared_ptr<F0Ps> _f);
+    virtual ~Assign() { std::cout << "~Assign()\n"; }
 
     void apply(std::vector<double> &vec);
 
@@ -57,10 +56,9 @@ private:
 class FLOW_API Initial : public Assign
 {
 public:
-    Initial(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f);
-    Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f);
-
-    ~Initial() { std::cout << "~Initial()\n"; }
+    Initial(std::shared_ptr<tbox::MshData> _msh, int no, std::shared_ptr<F0Ps> _f);
+    Initial(std::shared_ptr<tbox::MshData> _msh, std::string const &name, std::shared_ptr<F0Ps> _f);
+    virtual ~Initial() { std::cout << "~Initial()\n"; }
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
@@ -74,10 +72,9 @@ public:
 class FLOW_API Dirichlet : public Assign
 {
 public:
-    Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, F0Ps &_f);
-    Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, F0Ps &_f);
-
-    ~Dirichlet() { std::cout << "~Dirichlet()\n"; }
+    Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, std::shared_ptr<F0Ps>_f);
+    Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, std::shared_ptr<F0Ps> _f);
+    virtual ~Dirichlet() { std::cout << "~Dirichlet()\n"; }
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
diff --git a/flow/src/wF0Ps.cpp b/flow/src/wF0Ps.cpp
index 64dac02b638b4481c049c6a20c76e2b8f5257220..700102365f32ac5be11d2e2bc6ca0a2510c69195 100644
--- a/flow/src/wF0Ps.cpp
+++ b/flow/src/wF0Ps.cpp
@@ -16,7 +16,6 @@
 
 #include "wF0Ps.h"
 #include <sstream>
-using namespace tbox;
 using namespace flow;
 
 /**
@@ -33,9 +32,9 @@ void F0PsC::write(std::ostream &out) const
 }
 
 F0PsPhiInf::F0PsPhiInf(int _nDim, double _alpha, double _beta) : F0Ps(),
-                                                                     nDim(_nDim),
-                                                                     alpha(_alpha),
-                                                                     beta(_beta)
+                                                                 nDim(_nDim),
+                                                                 alpha(_alpha),
+                                                                 beta(_beta)
 {
     // Sanity check
     if (nDim != 2 && nDim != 3)
diff --git a/flow/src/wF0Ps.h b/flow/src/wF0Ps.h
index b8b176075821f169bc7996dba989b17f010ddad9..99531b8ff0b4f702b3d5607d6a10e30eb8f52233 100644
--- a/flow/src/wF0Ps.h
+++ b/flow/src/wF0Ps.h
@@ -28,11 +28,11 @@ namespace flow
  * @brief Scalar functions depending on spatial position
  * @authors Adrien Crovato
  */
-class FLOW_API F0Ps : public fwk::wObject
+class FLOW_API F0Ps : public fwk::wSharedObject
 {
 public:
 #ifndef SWIG
-    F0Ps() : wObject()
+    F0Ps() : wSharedObject()
     {
     }
 #endif
diff --git a/flow/src/wFreestream.cpp b/flow/src/wFreestream.cpp
index b319b5fe9fb61464db01269425bdde1d9056ef78..b68ecfd4c8093ebce3fb1dca9156cff21229f824 100644
--- a/flow/src/wFreestream.cpp
+++ b/flow/src/wFreestream.cpp
@@ -20,11 +20,11 @@
 using namespace tbox;
 using namespace flow;
 
-Freestream::Freestream(std::shared_ptr<MshData> _msh, int no, Fct1 &_f) : Group(_msh, no), f(_f)
+Freestream::Freestream(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<Fct1> _f) : Group(_msh, no), f(_f)
 {
 }
 
-Freestream::Freestream(std::shared_ptr<MshData> _msh, std::string const &name, Fct1 &_f) : Group(_msh, name), f(_f)
+Freestream::Freestream(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<Fct1> _f) : Group(_msh, name), f(_f)
 {
 }
 
diff --git a/flow/src/wFreestream.h b/flow/src/wFreestream.h
index 1725a2543da0bfc41b19b2c57b323d8bc471cd19..d0d3b7af78295471c90a670b5aa8e1408af1edd1 100644
--- a/flow/src/wFreestream.h
+++ b/flow/src/wFreestream.h
@@ -32,10 +32,10 @@ class FLOW_API Freestream : public tbox::Group
 {
 public:
 #ifndef SWIG
-    tbox::Fct1 &f; ///< vector valued function to compute flux
+    std::shared_ptr<tbox::Fct1> f; ///< vector valued function to compute flux
 #endif
-    Freestream(std::shared_ptr<tbox::MshData> _msh, int no, tbox::Fct1 &_f);
-    Freestream(std::shared_ptr<tbox::MshData> _msh, std::string const &name, tbox::Fct1 &_f);
+    Freestream(std::shared_ptr<tbox::MshData> _msh, int no, std::shared_ptr<tbox::Fct1> _f);
+    Freestream(std::shared_ptr<tbox::MshData> _msh, std::string const &name, std::shared_ptr<tbox::Fct1> _f);
     virtual ~Freestream() { std::cout << "~Freestream()\n"; }
 
 #ifndef SWIG
diff --git a/flow/src/wMedium.cpp b/flow/src/wMedium.cpp
index 9d4fecc2c3dda3016409ee22f5992aaaed0d80c5..beaab96a927f5f696465ce85298269f926ccee2d 100644
--- a/flow/src/wMedium.cpp
+++ b/flow/src/wMedium.cpp
@@ -21,9 +21,12 @@
 using namespace tbox;
 using namespace flow;
 
-Medium::Medium(std::shared_ptr<MshData> _msh, int no, Fct0 &_rho, Fct0 &_mach,
-               Fct0 &_cP, F0PsPhiInf &_phiInf) : Group(_msh, no), rho(_rho),
-                                                 mach(_mach), cP(_cP), phiInf(_phiInf)
+Medium::Medium(std::shared_ptr<MshData> _msh, int no,
+               std::shared_ptr<Fct0> _rho, std::shared_ptr<Fct0> _mach,
+               std::shared_ptr<Fct0> _cP,
+               std::shared_ptr<F0PsPhiInf> _phiInf) : Group(_msh, no),
+                                                      rho(_rho), mach(_mach),
+                                                      cP(_cP), phiInf(_phiInf)
 {
     createMap();
     std::cout << "Fluid is " << *tag << "with " << nCnt << " nodes."
@@ -31,7 +34,11 @@ Medium::Medium(std::shared_ptr<MshData> _msh, int no, Fct0 &_rho, Fct0 &_mach,
               << std::endl;
 }
 Medium::Medium(std::shared_ptr<MshData> _msh, std::string const &name,
-               Fct0 &_rho, Fct0 &_mach, Fct0 &_cP, F0PsPhiInf &_phiInf) : Group(_msh, name), rho(_rho), mach(_mach), cP(_cP), phiInf(_phiInf)
+               std::shared_ptr<Fct0> _rho, std::shared_ptr<Fct0> _mach,
+               std::shared_ptr<Fct0> _cP,
+               std::shared_ptr<F0PsPhiInf> _phiInf) : Group(_msh, name),
+                                                      rho(_rho), mach(_mach),
+                                                      cP(_cP), phiInf(_phiInf)
 {
     createMap();
     std::cout << "Fluid is " << *tag << "with " << nCnt << " nodes."
diff --git a/flow/src/wMedium.h b/flow/src/wMedium.h
index f84b815d56be1c01399978386196adc8f65c558a..77ec1eeb79c0227e83f70001855b093e9c8e6cb2 100644
--- a/flow/src/wMedium.h
+++ b/flow/src/wMedium.h
@@ -36,16 +36,16 @@ public:
     size_t nCnt;                                                ///< number of nodes
 
 #ifndef SWIG
-    tbox::Fct0 &rho;          ///< density
-    tbox::Fct0 &mach;         ///< Mach number
-    tbox::Fct0 &cP;           ///< pressure coefficient
-    flow::F0PsPhiInf &phiInf; ///< freestream potential
+    std::shared_ptr<tbox::Fct0> rho;    ///< density
+    std::shared_ptr<tbox::Fct0> mach;   ///< Mach number
+    std::shared_ptr<tbox::Fct0> cP;     ///< pressure coefficient
+    std::shared_ptr<F0PsPhiInf> phiInf; ///< freestream potential
 #endif
 
     std::vector<std::pair<tbox::Element *, std::vector<tbox::Element *>>> adjMap; ///< map of adjacent elements
 
-    Medium(std::shared_ptr<tbox::MshData> _msh, int no, tbox::Fct0 &_rho, tbox::Fct0 &_mach, tbox::Fct0 &_cP, flow::F0PsPhiInf &_phiInf);
-    Medium(std::shared_ptr<tbox::MshData> _msh, std::string const &name, tbox::Fct0 &_rho, tbox::Fct0 &_mach, tbox::Fct0 &_cP, flow::F0PsPhiInf &_phiInf);
+    Medium(std::shared_ptr<tbox::MshData> _msh, int no, std::shared_ptr<tbox::Fct0> _rho, std::shared_ptr<tbox::Fct0> _mach, std::shared_ptr<tbox::Fct0> _cP, std::shared_ptr<flow::F0PsPhiInf> _phiInf);
+    Medium(std::shared_ptr<tbox::MshData> _msh, std::string const &name, std::shared_ptr<tbox::Fct0> _rho, std::shared_ptr<tbox::Fct0> _mach, std::shared_ptr<tbox::Fct0> _cP, std::shared_ptr<flow::F0PsPhiInf> _phiInf);
     virtual ~Medium() { std::cout << "~Medium()\n"; }
 
 #ifndef SWIG
diff --git a/flow/src/wNewton.cpp b/flow/src/wNewton.cpp
index 61689817c40905b236e91d2751ca3916bd42299e..64a64095748b65dd4c1f19136bd28e5e50005430 100644
--- a/flow/src/wNewton.cpp
+++ b/flow/src/wNewton.cpp
@@ -273,11 +273,11 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
         Element *eU = p.second[0];
 
         // Subsonic contribution (drho*grad_phi*grad_psi + rho*grad_phi*grad_psi)
-        Eigen::MatrixXd Ae1 = e->buildDK(phi, fluid->rho);
-        Eigen::MatrixXd Ae2 = e->buildK(phi, fluid->rho);
+        Eigen::MatrixXd Ae1 = e->buildDK(phi, *fluid->rho);
+        Eigen::MatrixXd Ae2 = e->buildK(phi, *fluid->rho);
 
         // Supersonic contribution
-        double mach = fluid->mach.eval(*e, phi, 0);
+        double mach = fluid->mach->eval(*e, phi, 0);
         if (mach > mCOv)
         {
             // switching function and derivative
@@ -289,13 +289,13 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
             Ae2 *= 1 - mu;
 
             // 3rd term (mu*drhoU*grad_phi*grad_psi)
-            Eigen::MatrixXd Ae3 = mu * fluid->rho.evalD(*eU, phi, 0) * e->buildDs(phi, Fct0C(1.)) * eU->computeGrad(phi, 0).transpose() * eU->getVMem().getJinv(0) * eU->getVCache().getDsf(0);
+            Eigen::MatrixXd Ae3 = mu * fluid->rho->evalD(*eU, phi, 0) * e->buildDs(phi, Fct0C(1.)) * eU->computeGrad(phi, 0).transpose() * eU->getVMem().getJinv(0) * eU->getVCache().getDsf(0);
 
             // 4th term (mu*rhoU*grad_phi*grad_psi)
-            Eigen::MatrixXd Ae4 = mu * fluid->rho.eval(*eU, phi, 0) * e->buildK(phi, Fct0C(1.));
+            Eigen::MatrixXd Ae4 = mu * fluid->rho->eval(*eU, phi, 0) * e->buildK(phi, Fct0C(1.));
 
             // 5th term dmu*(rhoU-rho)*grad_phi*grad_psi
-            Eigen::MatrixXd Ae5 = dmu * (fluid->rho.eval(*eU, phi, 0) - fluid->rho.eval(*e, phi, 0)) * e->buildDK(phi, fluid->mach);
+            Eigen::MatrixXd Ae5 = dmu * (fluid->rho->eval(*eU, phi, 0) - fluid->rho->eval(*e, phi, 0)) * e->buildDK(phi, *fluid->mach);
 
             // Assembly (supersonic)
             tbb::spin_mutex::scoped_lock lock(mutex);
@@ -430,14 +430,14 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
         Element *eU = p.second[0];
 
         // Subsonic contribution
-        Eigen::VectorXd be = e->buildDs(phi, fluid->rho);
+        Eigen::VectorXd be = e->buildDs(phi, *fluid->rho);
 
         // Supersonic contribution
-        double mach = fluid->mach.eval(*e, phi, 0);
+        double mach = fluid->mach->eval(*e, phi, 0);
         if (mach > mCOv)
         {
             double mu = muCv * (1 - (mCOv * mCOv) / (mach * mach));
-            double rhoU = fluid->rho.eval(*eU, phi, 0);
+            double rhoU = fluid->rho->eval(*eU, phi, 0);
             Eigen::VectorXd beU = e->buildDs(phi, Fct0C(1.));
             // scale 1st and 2nd terms
             be *= 1 - mu;
@@ -463,7 +463,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
     {
         tbb::parallel_do(nBC->tag->elems.begin(), nBC->tag->elems.end(), [&](Element *e) {
             // Build elementary flux vector
-            Eigen::VectorXd be = e->builds(phi, nBC->f);
+            Eigen::VectorXd be = e->builds(phi, *nBC->f);
 
             // Assembly
             tbb::spin_mutex::scoped_lock lock(mutex);
diff --git a/flow/src/wPicard.cpp b/flow/src/wPicard.cpp
index cf0b346b1e758cbd87a83467cfc03144e9b6bb63..e2cafb254940b5dc64564dbfea3fac432976a6b4 100644
--- a/flow/src/wPicard.cpp
+++ b/flow/src/wPicard.cpp
@@ -202,7 +202,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
         Element *e = p.first;
 
         // Elementary stiffness
-        Eigen::MatrixXd Ae = e->buildK(phi, fluid->rho);
+        Eigen::MatrixXd Ae = e->buildK(phi, *fluid->rho);
 
         // Assembly
         tbb::spin_mutex::scoped_lock lock(mutex);
@@ -277,7 +277,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
     {
         tbb::parallel_do(nBC->tag->elems.begin(), nBC->tag->elems.end(), [&](Element *e) {
             // Build elementary flux vector
-            Eigen::VectorXd be = e->builds(phi, nBC->f);
+            Eigen::VectorXd be = e->builds(phi, *nBC->f);
 
             // Assembly
             tbb::spin_mutex::scoped_lock lock(mutex);
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index 9e054fcedc2d26ad76223790f137b24e77cd23cd..7b376b1e51491e8ec8c3901adb0781e403fd774d 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -169,7 +169,7 @@ void Solver::computeLoad()
         // Map each element to its value
         std::map<Element *, double> euMap;
         for (auto e : sur->groups[0]->tag->elems)
-            euMap[e] = sur->integrate(e, phi, pbl->medium->cP);
+            euMap[e] = sur->integrate(e, phi, *pbl->medium->cP);
 
         // Compute normalized load (i.e. load / pressure)
         for (size_t i = 0; i < sur->nodes.size(); ++i)
@@ -228,7 +228,7 @@ void Solver::computeFlow()
     tbb::parallel_do(fluid->neMap.begin(), fluid->neMap.end(), [&](std::pair<Node *, std::vector<Element *>> p) {
         Node *n = p.first;
         // Perturbation potential
-        vPhi[n->row] = phi[n->row] - fluid->phiInf.eval(n->pos);
+        vPhi[n->row] = phi[n->row] - fluid->phiInf->eval(n->pos);
         // Velocity (averaged at nodes)
         U[n->row] = Eigen::Vector3d::Zero();
         for (auto e : p.second)
@@ -241,17 +241,17 @@ void Solver::computeFlow()
         // Density (averaged at nodes)
         rho[n->row] = 0.;
         for (auto e : p.second)
-            rho[n->row] += fluid->rho.eval(*e, phi, 0);
+            rho[n->row] += fluid->rho->eval(*e, phi, 0);
         rho[n->row] /= p.second.size();
         // Mach number (averaged at nodes)
         M[n->row] = 0.;
         for (auto e : p.second)
-            M[n->row] += fluid->mach.eval(*e, phi, 0);
+            M[n->row] += fluid->mach->eval(*e, phi, 0);
         M[n->row] /= p.second.size();
         // Pressure coefficient (averaged at nodes)
         Cp[n->row] = 0.;
         for (auto e : p.second)
-            Cp[n->row] += fluid->cP.eval(*e, phi, 0);
+            Cp[n->row] += fluid->cP->eval(*e, phi, 0);
         Cp[n->row] /= p.second.size();
     });
 
diff --git a/heat/_src/heatw.i b/heat/_src/heatw.i
index 33dce9ed4c3dcfac974c00a2d73b81ff7c86cf21..cb15617bab6234c876d37d5ef0c6f38ce62143de 100644
--- a/heat/_src/heatw.i
+++ b/heat/_src/heatw.i
@@ -87,6 +87,10 @@ threads="1"
 %include "wSource.h"
 %include "wBoundary.h"
 %include "wPeriodic.h"
+%shared_ptr(heat::CompiledFct1a);
+%shared_ptr(heat::CompiledFct2a);
+%shared_ptr(heat::CompiledFct1b);
+%shared_ptr(heat::CompiledFct2b);
 %include "wCompiledFct.h"
 %include "wExtractor.h"
 
diff --git a/heat/src/wBoundary.cpp b/heat/src/wBoundary.cpp
index 5d7aeba59d80aef42b028994d6f85d47d2acc017..1963bfccfe22f57f6fe3b314e921953f6d10bd4b 100644
--- a/heat/src/wBoundary.cpp
+++ b/heat/src/wBoundary.cpp
@@ -19,12 +19,12 @@
 #include "wTag.h"
 using namespace heat;
 
-Boundary::Boundary(std::shared_ptr<MshData> _msh, int no, Fct0 &_f) : Group(_msh, no), f(_f)
+Boundary::Boundary(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<Fct0> _f) : Group(_msh, no), f(_f)
 {
     //pbl.bnds.push_back(this);
 }
 
-Boundary::Boundary(std::shared_ptr<MshData> _msh, std::string const &name, Fct0 &_f) : Group(_msh, name), f(_f)
+Boundary::Boundary(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<Fct0> _f) : Group(_msh, name), f(_f)
 {
     //pbl.bnds.push_back(this);
 }
diff --git a/heat/src/wBoundary.h b/heat/src/wBoundary.h
index e0ffaf88c68d8da254945ad3b69f669ce46b7d04..93476b1115dcfc266af89a3737380624361c78e9 100644
--- a/heat/src/wBoundary.h
+++ b/heat/src/wBoundary.h
@@ -34,10 +34,10 @@ class HEAT_API Boundary : public Group
 {
 public:
 #ifndef SWIG
-    Fct0 &f;
+    std::shared_ptr<Fct0> f;
 #endif
-    Boundary(std::shared_ptr<MshData> _msh, int no, Fct0 &_f);
-    Boundary(std::shared_ptr<MshData> _msh, std::string const &name, Fct0 &_f);
+    Boundary(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<Fct0> _f);
+    Boundary(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<Fct0> _f);
     virtual ~Boundary() { std::cout << "~Boundary()\n"; }
 
 #ifndef SWIG
diff --git a/heat/src/wMedium.cpp b/heat/src/wMedium.cpp
index f98bd2b108474d17590f0b5bbf058dda6bfcf618..6e0259251bf23b503dc6693bced97b059380ec14 100644
--- a/heat/src/wMedium.cpp
+++ b/heat/src/wMedium.cpp
@@ -21,13 +21,15 @@
 using namespace heat;
 
 Medium::Medium(std::shared_ptr<MshData> _msh, int no,
-               Fct2 &_k, double _rhoc) : Group(_msh, no), k(_k), rhoc(_rhoc)
+               std::shared_ptr<Fct2> _k, double _rhoc) : Group(_msh, no),
+                                                         k(_k), rhoc(_rhoc)
 {
     //pbl.media.push_back(this);
 }
 
 Medium::Medium(std::shared_ptr<MshData> _msh, std::string const &name,
-               Fct2 &_k, double _rhoc) : Group(_msh, name), k(_k), rhoc(_rhoc)
+               std::shared_ptr<Fct2> _k, double _rhoc) : Group(_msh, name),
+                                                         k(_k), rhoc(_rhoc)
 {
     std::cout << "tag=" << *tag << std::endl;
     //pbl.media.push_back(this);
diff --git a/heat/src/wMedium.h b/heat/src/wMedium.h
index b5ee1751c1a95e560ff83ce179d708d24a8c74b6..30e0f25e46ca4eecc9c0b1ef2b1f9bda0dccea1d 100644
--- a/heat/src/wMedium.h
+++ b/heat/src/wMedium.h
@@ -34,11 +34,11 @@ class HEAT_API Medium : public Group
 {
 public:
 #ifndef SWIG
-    Fct2 &k;
+    std::shared_ptr<Fct2> k;
     double rhoc;
 #endif
-    Medium(std::shared_ptr<MshData> _msh, int no, Fct2 &_k, double _rhoc = 1.0);
-    Medium(std::shared_ptr<MshData> _msh, std::string const &name, Fct2 &_k, double _rhoc = 1.0);
+    Medium(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<Fct2> _k, double _rhoc = 1.0);
+    Medium(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<Fct2> _k, double _rhoc = 1.0);
     virtual ~Medium() { std::cout << "~Medium()\n"; }
 
 #ifndef SWIG
diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp
index fd6c8f78994467ebbc2f37069336a8ee246f7cab..33583c528f0ccddb60e6a86089a86e0044b417cc 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -268,7 +268,7 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
                              //std::cout << "processing element #" << e->no << "\n";
 
                              // (flux moyen . volume) sur l'element
-                             Eigen::VectorXd qV = e->computeqV(T1, mat->k);
+                             Eigen::VectorXd qV = e->computeqV(T1, *mat->k);
 
                              double V = e->computeV();
 
@@ -285,7 +285,7 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
                                  kgradT[i] = Eigen::Vector3d(qV(0) / V, qV(1) / V, 0);
 
                                  // mean k over the element
-                                 Eigen::MatrixXd kmean = e->computeV(T1, mat->k) / V;
+                                 Eigen::MatrixXd kmean = e->computeV(T1, *mat->k) / V;
                                  kmoy11[i] = kmean(0, 0);
                                  kmoy22[i] = kmean(1, 1);
                                  kmoy12[i] = kmean(0, 1);
@@ -383,7 +383,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
                                      return;
                                  //std::cout << "processing element #" << e->no << "\n";
 
-                                 Eigen::MatrixXd Ke = e->buildK(T1, mat->k, false);
+                                 Eigen::MatrixXd Ke = e->buildK(T1, *mat->k, false);
 
                                  // assembly
                                  tbb::spin_mutex::scoped_lock lock(mutex);
@@ -416,7 +416,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
                           [&](Element *e) {
                               if (e->type() != ELTYPE::TRI3)
                                   return;
-                              Eigen::MatrixXd Ke = e->buildK(T1, mat->k, true); // bidon
+                              Eigen::MatrixXd Ke = e->buildK(T1, *mat->k, true); // bidon
                           });
         }
 
@@ -425,7 +425,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
         if (verbose > 1)
             std::cout << "executing job list...\n";
         for (auto mat : pbl->media)
-            mat->k.evalall();
+            mat->k->evalall();
 
         // ***  perform the assembly with the results ("true run")
 
@@ -441,7 +441,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
                           [&](Element *e) {
                               if (e->type() != ELTYPE::TRI3)
                                   return;
-                              Eigen::MatrixXd Ke = e->buildK(T1, mat->k, false);
+                              Eigen::MatrixXd Ke = e->buildK(T1, *mat->k, false);
 
                               // assembly
                               //tbb::spin_mutex::scoped_lock lock(mutex);
@@ -608,7 +608,7 @@ void Solver::buildqint(Eigen::VectorXd &qint)
                                  return;
                              //std::cout << "processing element #" << e->no << "\n";
 
-                             Eigen::VectorXd qe = e->computeqint(T1, mat->k);
+                             Eigen::VectorXd qe = e->computeqint(T1, *mat->k);
 
                              // assembly
                              tbb::spin_mutex::scoped_lock lock(mutex);
@@ -659,7 +659,7 @@ void Solver::builds(Eigen::Map<Eigen::VectorXd> &s)
                              //std::cout << "processing element #" << e->no << "\n";
 
                              // ** se vector => s vector
-                             Eigen::VectorXd se = e->builds(T1, src->f);
+                             Eigen::VectorXd se = e->builds(T1, *src->f);
 
                              // assembly
                              tbb::spin_mutex::scoped_lock lock(mutex);
@@ -699,7 +699,7 @@ void Solver::buildq(Eigen::Map<Eigen::VectorXd> &q)
                              //std::cout << "processing element #" << e->no << "\n";
 
                              // ** qe vector => q vector
-                             Eigen::VectorXd qe = e->builds(T1, bnd->f);
+                             Eigen::VectorXd qe = e->builds(T1, *bnd->f);
 
                              // assembly
                              tbb::spin_mutex::scoped_lock lock(mutex);
diff --git a/heat/src/wSource.cpp b/heat/src/wSource.cpp
index 4bea57003be8c0048da2cf6c07dc60a7ab7b6faa..9e588260c679ba506db92a920311313d66cd9411 100644
--- a/heat/src/wSource.cpp
+++ b/heat/src/wSource.cpp
@@ -19,12 +19,12 @@
 #include "wTag.h"
 using namespace heat;
 
-Source::Source(std::shared_ptr<MshData> _msh, int no, Fct0 &_f) : Group(_msh, no), f(_f)
+Source::Source(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<Fct0> _f) : Group(_msh, no), f(_f)
 {
     //pbl.srcs.push_back(this);
 }
 
-Source::Source(std::shared_ptr<MshData> _msh, std::string const &name, Fct0 &_f) : Group(_msh, name), f(_f)
+Source::Source(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<Fct0> _f) : Group(_msh, name), f(_f)
 {
     //pbl.srcs.push_back(this);
 }
diff --git a/heat/src/wSource.h b/heat/src/wSource.h
index 2dba1b2f9c41be9f616872c2e7eac9a9489c78fd..9de9b77900b3c10bf7bde59039ea901f519c5610 100644
--- a/heat/src/wSource.h
+++ b/heat/src/wSource.h
@@ -36,10 +36,10 @@ class HEAT_API Source : public Group
 {
 public:
 #ifndef SWIG
-    Fct0 &f;
+    std::shared_ptr<Fct0> f;
 #endif
-    Source(std::shared_ptr<MshData> _msh, int no, Fct0 &_f);
-    Source(std::shared_ptr<MshData> _msh, std::string const &name, Fct0 &_f);
+    Source(std::shared_ptr<MshData> _msh, int no, std::shared_ptr<Fct0> _f);
+    Source(std::shared_ptr<MshData> _msh, std::string const &name, std::shared_ptr<Fct0> _f);
     virtual ~Source() { std::cout << "~Source()\n"; }
 
 #ifndef SWIG
diff --git a/tbox/_src/tboxw.i b/tbox/_src/tboxw.i
index db3e8227cd87b721661d0676e95910e1c6875570..0f1f4045450326cd77b0332057b73e35efbe63f4 100644
--- a/tbox/_src/tboxw.i
+++ b/tbox/_src/tboxw.i
@@ -64,34 +64,39 @@ namespace std {
    %template(std_vector_string) std::vector<std::string>;
 }
 
-%pythonappend tbox::PwLf::PwLf "self.thisown=0"   // faire des compteurs de ref au lieu de ça
-
+// %nodefault: for pure virtual (abstract) class
+// %feature("director"): allows C++ to see inherited classe extensions implemented in python
+// %pythonappend ... self.__disown__(): set thisown=0, which directly increments the python reference count by 1 (only for directors)
 %feature("director") tbox::Fct0XYZ;
-%pythonappend tbox::Fct0XYZ::Fct0XYZ "self.__disown__()"    // only for directors
-%feature("director") tbox::Fct0U;
-%pythonappend tbox::Fct0U::Fct0U "self.__disown__()"    // only for directors
-%pythonappend tbox::Fct0C::Fct0C "self.thisown=0"
-%nodefault tbox::Fct0;     // Fct0 is pure virtual but swig doesn't know it.
+%pythonappend tbox::Fct0XYZ::Fct0XYZ "self.__disown__()"
+%shared_ptr(tbox::Fct0);
+%shared_ptr(tbox::Fct0C);
+%shared_ptr(tbox::Fct0XYZ);
+%shared_ptr(tbox::Fct0U);
+%shared_ptr(tbox::PwLf)
+%nodefault tbox::Fct0; // Fct0 is pure virtual but swig doesn't know it.
 %include "wFct0.h"
 
-%pythonappend tbox::Fct1C::Fct1C "self.thisown=0"
-%nodefault tbox::Fct1;     // Fct1 is pure virtual but swig doesn't know it.
+%shared_ptr(tbox::Fct1);
+%shared_ptr(tbox::Fct1C);
+%nodefault tbox::Fct1; // Fct1 is pure virtual but swig doesn't know it.
 %include "wFct1.h"
 
 %feature("director") tbox::Fct2;
 %feature("director") tbox::Fct2XYZ;
-%pythonappend tbox::Fct2XYZ::Fct2XYZ "self.__disown__()"    // only for directors
 %feature("director") tbox::Fct2U;
-%pythonappend tbox::Fct2U::Fct2U "self.__disown__()"    // only for directors
 %feature("director") tbox::Fct2UdU;
-%pythonappend tbox::Fct2UdU::Fct2UdU "self.__disown__()"    // only for directors
-%pythonappend tbox::Fct2C::Fct2C "self.thisown=0"
-%pythonappend tbox::Fct2PwLf::Fct2PwLf "self.thisown=0"
-%nodefault tbox::Fct2;     // Fct2 is pure virtual but swig doesn't know it.
+%shared_ptr(tbox::Fct2);
+%shared_ptr(tbox::Fct2C);
+%shared_ptr(tbox::Fct2XYZ);
+%shared_ptr(tbox::Fct2U);
+%shared_ptr(tbox::Fct2UdU);
+%shared_ptr(tbox::Fct2PwLf);
+%nodefault tbox::Fct2; // Fct2 is pure virtual but swig doesn't know it.
 %include "wFct2.h"
 
 %include "wNode.h"
-%nodefault tbox::Element;     // Element is pure virtual but swig doesn't know it.
+%nodefault tbox::Element; // Element is pure virtual but swig doesn't know it.
 %warnfilter(509); // Shadowed overloaded method MATTYPE. Code compiling and running.
 %include "wElement.h"
 %warnfilter(+509);
diff --git a/tbox/src/wFct0.h b/tbox/src/wFct0.h
index a46d59cb24f504314bccb68273dcc3ade6d166d1..40162ac548fc4d0afdfd5bdaf89705807c3330cc 100644
--- a/tbox/src/wFct0.h
+++ b/tbox/src/wFct0.h
@@ -26,12 +26,12 @@
 namespace tbox
 {
 
-class TBOX_API PwLf : public fwk::wObject
+class TBOX_API PwLf : public fwk::wSharedObject
 {
     std::map<double, double> pts;
 
 public:
-    PwLf() {}
+    PwLf() : wSharedObject() {}
     void add(double x, double y);
     double eval(double x) const;
 #ifndef SWIG
@@ -46,11 +46,11 @@ public:
  * evalD() evaluates the derivative at the kth Gauss point
  * @authors Romain Boman, Adrien Crovato
  */
-class TBOX_API Fct0 : public fwk::wObject
+class TBOX_API Fct0 : public fwk::wSharedObject
 {
 public:
 #ifndef SWIG
-    Fct0() : wObject()
+    Fct0() : wSharedObject()
     {
     }
 #endif
diff --git a/tbox/src/wFct1.h b/tbox/src/wFct1.h
index 0edd85849a14bdd95a99518ef7a38f33b9e38505..feabf1c2e1afb3bd07615bc78b0e4fb31039280d 100644
--- a/tbox/src/wFct1.h
+++ b/tbox/src/wFct1.h
@@ -31,11 +31,11 @@ namespace tbox
  * eval() evaluates the function at kth Gauss point
  * @authors Romain Boman, Adrien Crovato
  */
-class TBOX_API Fct1 : public fwk::wObject
+class TBOX_API Fct1 : public fwk::wSharedObject
 {
 public:
 #ifndef SWIG
-    Fct1() : wObject()
+    Fct1() : wSharedObject()
     {
     }
 #endif
diff --git a/tbox/src/wFct2.cpp b/tbox/src/wFct2.cpp
index b510d11475606713a843e6a86dd203982f9250ec..85844bdadb595de2a1ea09541ff023381db81aa4 100644
--- a/tbox/src/wFct2.cpp
+++ b/tbox/src/wFct2.cpp
@@ -121,9 +121,9 @@ void Fct2UdU::write(std::ostream &out) const
 void Fct2PwLf::eval(double u, Eigen::MatrixXd &v, bool fake) const
 {
     v.resize(2, 2);
-    v(0, 0) = v11.eval(u);
-    v(1, 1) = v22.eval(u);
-    v(1, 0) = v(0, 1) = v12.eval(u);
+    v(0, 0) = v11->eval(u);
+    v(1, 1) = v22->eval(u);
+    v(1, 0) = v(0, 1) = v12->eval(u);
 }
 
 void Fct2PwLf::write(std::ostream &out) const
diff --git a/tbox/src/wFct2.h b/tbox/src/wFct2.h
index b4ac8f9a6a91a49b2e1caa2a6ee60e8c0f7f9f9b..9cffc8c2d6d2f4fe90e81b56fab34e643f3ec4b7 100644
--- a/tbox/src/wFct2.h
+++ b/tbox/src/wFct2.h
@@ -19,6 +19,7 @@
 
 #include "tbox.h"
 #include "wObject.h"
+#include <memory>
 #include <vector>
 #include <Eigen/Dense>
 
@@ -31,11 +32,11 @@ namespace tbox
  * eval() evaluates the function at kth Gauss point
  * @authors Romain Boman
  */
-class TBOX_API Fct2 : public fwk::wObject
+class TBOX_API Fct2 : public fwk::wSharedObject
 {
 public:
 #ifndef SWIG
-    Fct2() : wObject()
+    Fct2() : wSharedObject()
     {
     }
 #endif
@@ -121,12 +122,12 @@ public:
  */
 class TBOX_API Fct2PwLf : public Fct2U
 {
-    PwLf &v11;
-    PwLf &v22;
-    PwLf &v12;
+    std::shared_ptr<PwLf> v11;
+    std::shared_ptr<PwLf> v22;
+    std::shared_ptr<PwLf> v12;
 
 public:
-    Fct2PwLf(PwLf &_v11, PwLf &_v22, PwLf &_v12) : Fct2U(), v11(_v11), v22(_v22), v12(_v12) {}
+    Fct2PwLf(std::shared_ptr<PwLf> _v11, std::shared_ptr<PwLf> _v22, std::shared_ptr<PwLf> _v12) : Fct2U(), v11(_v11), v22(_v22), v12(_v12) {}
 
     virtual void eval(double u, Eigen::MatrixXd &v, bool fake) const override;