diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6beafc1dfe06defe2646d43a2004ae8711aea19d..55d4ead63149b74f4e7a5fe815fa4fd5c3b5c843 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -32,7 +32,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 3e969f4799c26c62549499d49f29ebdcf5fe7012..b3f6bd80dfec11ff02486b60bb666f53442ea0bd 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-2021
 (c) ULg - A&M",
-directors="1",
+directors="0",
 threads="1"
 ) floww
 %{
 
-#include <string>
-#include <sstream>
-#include <typeinfo>
 #include "flow.h"
 
 #include "fwkw.h"
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index 991721ff36e3fb6d317562733df113c2997b054c..29c598dd3b7e4a47b7fed6693f99508dacd516bf 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -57,9 +57,9 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<LinearSolver> _lin
     std::cout << "**                          | __| | | /  \\ \\ \\_/ /                           **" << std::endl;
     std::cout << "**                          |_|   |_| \\__/  \\/ \\/                            **" << std::endl;
     std::cout << "*******************************************************************************" << std::endl;
-    std::cout << "** Hi! My name is Flow v1.7-2011                                             **" << std::endl;
+    std::cout << "** Hi! My name is Flow v1.8-21.06                                            **" << std::endl;
     std::cout << "** Adrien Crovato & Romain Boman                                             **" << std::endl;
-    std::cout << "** ULiege 2018-2020                                                          **" << std::endl;
+    std::cout << "** ULiege 2018-2021                                                          **" << std::endl;
     std::cout << "*******************************************************************************" << std::endl
               << std::endl;
 
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 7bca36520e0f46825597482e9c8043974721e470..4a308e88ae73602c2e7dcbb78cddacc8ba955455 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -263,7 +263,7 @@ void Solver::start(MshExport *mshWriter)
                              //std::cout << "processing element #" << e->no << "\n";
 
                              // (flux moyen . volume) sur l'element
-                             Eigen::VectorXd qV = HeatTerm::computeFlux(*e, T1, mat->k);
+                             Eigen::VectorXd qV = HeatTerm::computeFlux(*e, T1, *mat->k);
 
                              double V = e->getVol();
 
@@ -280,7 +280,7 @@ void Solver::start(MshExport *mshWriter)
                                  kgradT[i] = Eigen::Vector3d(qV(0) / V, qV(1) / V, 0);
 
                                  // mean k over the element
-                                 Eigen::MatrixXd kmean = HeatTerm::computeMatrix(*e, T1, mat->k) / V;
+                                 Eigen::MatrixXd kmean = HeatTerm::computeMatrix(*e, T1, *mat->k) / V;
                                  kmoy11[i] = kmean(0, 0);
                                  kmoy22[i] = kmean(1, 1);
                                  kmoy12[i] = kmean(0, 1);
@@ -378,7 +378,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
                                      return;
                                  //std::cout << "processing element #" << e->no << "\n";
 
-                                 Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, mat->k, false);
+                                 Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, *mat->k, false);
 
                                  // assembly
                                  tbb::spin_mutex::scoped_lock lock(mutex);
@@ -411,7 +411,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
                           [&](Element *e) {
                               if (e->type() != ELTYPE::TRI3)
                                   return;
-                              Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, mat->k, true); // bidon
+                              Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, *mat->k, true); // bidon
                           });
         }
 
@@ -420,7 +420,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")
 
@@ -436,7 +436,7 @@ void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map
                           [&](Element *e) {
                               if (e->type() != ELTYPE::TRI3)
                                   return;
-                              Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, mat->k, false);
+                              Eigen::MatrixXd Ke = HeatTerm::build(*e, T1, *mat->k, false);
 
                               // assembly
                               //tbb::spin_mutex::scoped_lock lock(mutex);
@@ -602,7 +602,7 @@ void Solver::buildqint(Eigen::VectorXd &qint)
                                  return;
                              //std::cout << "processing element #" << e->no << "\n";
 
-                             Eigen::VectorXd qe = HeatTerm::build2(*e, T1, mat->k);
+                             Eigen::VectorXd qe = HeatTerm::build2(*e, T1, *mat->k);
 
                              // assembly
                              tbb::spin_mutex::scoped_lock lock(mutex);
@@ -653,7 +653,7 @@ void Solver::builds(Eigen::Map<Eigen::VectorXd> &s)
                              //std::cout << "processing element #" << e->no << "\n";
 
                              // ** se vector => s vector
-                             Eigen::VectorXd se = HeatTerm::build(*e, T1, src->f);
+                             Eigen::VectorXd se = HeatTerm::build(*e, T1, *src->f);
 
                              // assembly
                              tbb::spin_mutex::scoped_lock lock(mutex);
@@ -693,7 +693,7 @@ void Solver::buildq(Eigen::Map<Eigen::VectorXd> &q)
                              //std::cout << "processing element #" << e->no << "\n";
 
                              // ** qe vector => q vector
-                             Eigen::VectorXd qe = HeatTerm::build(*e, T1, bnd->f);
+                             Eigen::VectorXd qe = HeatTerm::build(*e, 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 5203eb403e367357cde94f6745f98d0bb090730c..58f05923f227df10ac295e802c7b94b670190fa9 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 eb6d41661c5309808cb262d89b42cdea41a2406d..591b6d515c7dfe288d97dea2456d6a57b9303294 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
@@ -45,11 +45,11 @@ public:
  * eval() evaluates the function at kth Gauss point
  * @authors Romain Boman
  */
-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 0a3c546dc6a54630f9a17075a1d39480a92a1ccc..f82320d1fa59a4ef2d6ba1844fefe543b178e6b4 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 2bf6ac7adb8454b9a666494f23b02a375b73fa31..9abe21bcab28f8a8edb6f3ea6c55f562b3814421 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 4b902552555a0b2e521b3ddbae3f8ccd44afa526..ed200d74d3bae6a569e369132c21a3503a671ca7 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;