From 40171e0cbdbdaed3b7bf95c7fef413418bdee884 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Tue, 21 Apr 2020 16:26:14 +0200
Subject: [PATCH] Add share_ptr to tbox, heat and flow functions

---
 CMakeLists.txt           |  2 +-
 flow/_src/floww.i        | 37 +++++++++++++------------------------
 flow/src/wAdjoint.cpp    | 16 ++++++++--------
 flow/src/wAssign.cpp     | 18 +++++++++---------
 flow/src/wAssign.h       | 23 ++++++++++-------------
 flow/src/wF0Ps.cpp       |  7 +++----
 flow/src/wF0Ps.h         |  4 ++--
 flow/src/wFreestream.cpp |  4 ++--
 flow/src/wFreestream.h   |  6 +++---
 flow/src/wMedium.cpp     | 15 +++++++++++----
 flow/src/wMedium.h       | 12 ++++++------
 flow/src/wNewton.cpp     | 20 ++++++++++----------
 flow/src/wPicard.cpp     |  4 ++--
 flow/src/wSolver.cpp     | 10 +++++-----
 heat/_src/heatw.i        |  4 ++++
 heat/src/wBoundary.cpp   |  4 ++--
 heat/src/wBoundary.h     |  6 +++---
 heat/src/wMedium.cpp     |  6 ++++--
 heat/src/wMedium.h       |  6 +++---
 heat/src/wSolver.cpp     | 18 +++++++++---------
 heat/src/wSource.cpp     |  4 ++--
 heat/src/wSource.h       |  6 +++---
 tbox/_src/tboxw.i        | 37 +++++++++++++++++++++----------------
 tbox/src/wFct0.h         |  8 ++++----
 tbox/src/wFct1.h         |  4 ++--
 tbox/src/wFct2.cpp       |  6 +++---
 tbox/src/wFct2.h         | 13 +++++++------
 27 files changed, 152 insertions(+), 148 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3306b7a6..0b553686 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 9a6b59d5..9bfe0570 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 cb04c1c5..a31bed41 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 919faf26..f97a4ebe 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 81b8ab04..b6a9c03a 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 64dac02b..70010236 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 b8b17607..99531b8f 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 b319b5fe..b68ecfd4 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 1725a254..d0d3b7af 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 9d4fecc2..beaab96a 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 f84b815d..77ec1eeb 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 61689817..64a64095 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 cf0b346b..e2cafb25 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 9e054fce..7b376b1e 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 33dce9ed..cb15617b 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 5d7aeba5..1963bfcc 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 e0ffaf88..93476b11 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 f98bd2b1..6e025925 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 b5ee1751..30e0f25e 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 fd6c8f78..33583c52 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 4bea5700..9e588260 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 2dba1b2f..9de9b779 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 db3e8227..0f1f4045 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 a46d59cb..40162ac5 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 0edd8584..feabf1c2 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 b510d114..85844bda 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 b4ac8f9a..9cffc8c2 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;
 
-- 
GitLab