diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 7477ee22ddfcc58f239d01b160a2c51a2a181c62..7f76ede9df3262e6f3f6c53058e3d9069acfe966 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -110,15 +110,9 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
     for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){
       residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
     }
-    if(convectiveFlux)
-      delete convectiveFlux;
-    if(diffusiveFlux)
-      delete diffusiveFlux;
   }
   if(sourceTerm){
     residual.gemm(group.getSourceRedistributionMatrix(),Source);
-    //    FIXME (JF) : for now TEST TEST
-    //    delete sourceTerm;
   }
 }
 
@@ -240,10 +234,6 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
 
   if(riemannSolver || diffusiveFluxLeft)
     residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
-
-  // D ) delete the dataCacheDouble provided by the law
-  if(riemannSolver)
-    delete riemannSolver;
 }
 
 void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
@@ -453,8 +443,6 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw,			// conse
    }
    
    for (int i=0;i<sol._dataSize;i++){
-     //     printf("tempSol[%d] = %g\n",i,(*tempSol._data)(i));
-     //     memcp
      (*sol._data)(i)=(*Unp._data)(i);
    }
    for(int i=0;i<nStages;i++){
@@ -545,9 +533,6 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
   }
   // ----- 3 ---- do the redistribution at face nodes using BLAS3
   residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
-  delete boundaryTerm;
-  if (dirichlet) delete dirichlet;
-  if (neumann) delete neumann;
 }
 
 /*
@@ -753,11 +738,8 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
       double c = (*maxConvectiveSpeed)(0,0);
       for (int k=1;k<group.getNbNodes();k++) c = std::max((*maxConvectiveSpeed)(k,0), c);
       spectralRadius += 4.0*c*l_red/L; 
-      //      printf("coucou %g %g %g %g\n",c,l_red,L,spectralRadius);
     }
     DT[iElement] = 1./spectralRadius;
   }
-  delete maximumDiffusivity;
-  delete maxConvectiveSpeed;
 }
 
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index 853b6756bbff9af99428e0e3c16a1f8599298010..c9a71dc9fc97fbc95c7dd4b831e500835cbc0f81 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/Solver/dgConservationLaw.cpp
@@ -10,7 +10,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
     dgConservationLaw *_claw;
     public:
     term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
-      dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()),
+      dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()),
       cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
       solutionRight(cacheMapRight.provideData("Solution")),
       outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
@@ -44,7 +44,7 @@ class dgBoundarySymmetry : public dgBoundaryCondition {
     dgConservationLaw *_claw;
   public:
     term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft):
-      dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()), _claw(claw)
+      dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()), _claw(claw)
     {
       riemannSolver=_claw->newRiemannSolver(cacheMapLeft,cacheMapLeft);
       riemannSolver->addMeAsDependencyOf(this);
@@ -69,7 +69,7 @@ class dgBoundaryCondition0Flux : public dgBoundaryCondition {
   class term : public dataCacheDouble {
   public:
     term(dgConservationLaw *claw,dataCacheMap &cacheMapLeft):
-      dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()) {}
+      dataCacheDouble(cacheMapLeft,cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()) {}
     void _eval() {
     }
   };
@@ -95,6 +95,7 @@ class dgBoundaryCondition::dirichlet_ : public dataCacheDouble {
   dataCacheDouble &sol;
 public:
   dirichlet_(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this)){}
   void _eval () { 
     int nQP = sol().size1();
@@ -112,6 +113,7 @@ class dgBoundaryCondition::neumann_ : public dataCacheDouble {
   dataCacheDouble *diffusiveFlux;
 public:
   neumann_(dataCacheMap &cacheMap, dgConservationLaw *claw):
+    dataCacheDouble(cacheMap),
     _claw (claw),
     sol(cacheMap.get("Solution",this)),
     normals(cacheMap.get("Normals",this)){
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index 76b61eb5df8a539aa6d13b221a0e3edb38a82b52..0217e9412854cf90f065158c6137ac3df984bf39 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -10,7 +10,7 @@ class dgConservationLawAdvection::advection : public dataCacheDouble {
   dataCacheDouble &sol, &v;
   public:
   advection(std::string vFunctionName, dataCacheMap &cacheMap):
-    dataCacheDouble(cacheMap.getNbEvaluationPoints(),3),
+    dataCacheDouble(cacheMap,cacheMap.getNbEvaluationPoints(),3),
     sol(cacheMap.get("Solution",this)),
     v(cacheMap.get(vFunctionName,this))
   {};
@@ -26,7 +26,7 @@ class dgConservationLawAdvection::riemann : public dataCacheDouble {
   dataCacheDouble &normals, &solLeft, &solRight,&v;
   public:
   riemann(std::string vFunctionName, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
-    dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),2),
+    dataCacheDouble(cacheMapLeft,cacheMapLeft.getNbEvaluationPoints(),2),
     normals(cacheMapLeft.get("Normals", this)),
     solLeft(cacheMapLeft.get("Solution", this)),
     solRight(cacheMapRight.get("Solution", this)),
@@ -49,6 +49,7 @@ class dgConservationLawAdvection::diffusion : public dataCacheDouble {
   dataCacheDouble &solgrad, &nu;
   public:
   diffusion(std::string nuFunctionName, dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     solgrad(cacheMap.get("SolutionGradient",this)),
     nu(cacheMap.get(nuFunctionName,this))
   {};
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index 22e974fde1585e71426d24c757df928122fe275c..fc9a1331f58719b36e9f9f7f0c911a08772a70be 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -328,6 +328,7 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
   dataCacheDouble &sol;
   public:
   advection(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this))
   {};
   void _eval () { 
@@ -365,6 +366,7 @@ class dgPerfectGasLaw2d::diffusion : public dataCacheDouble {
   diffusion(dataCacheMap &cacheMap, 
 	    const std::string &muFunctionName,
 	    const std::string &kappaFunctionName ):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this)),
     grad(cacheMap.get("SolutionGradient",this)),
     mu(cacheMap.get(muFunctionName,this)),
@@ -425,6 +427,7 @@ class dgPerfectGasLaw2d::source : public dataCacheDouble {
   dataCacheDouble &sol,&s;
 public:
   source(dataCacheMap &cacheMap, const std::string &sourceFunctionName):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this)),
     s(cacheMap.get(sourceFunctionName,this))
   {};
@@ -446,6 +449,7 @@ class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
   double _presMin, _rhoMin;
 public:
   clipToPhysics(dataCacheMap &cacheMap, double presMin, double rhoMin):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this))
   {
     _presMin=presMin;
@@ -479,6 +483,7 @@ class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
   dataCacheDouble &normals, &solL, &solR;
   public:
   riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+    dataCacheDouble(cacheMapLeft),
     normals(cacheMapLeft.get("Normals", this)),
     solL(cacheMapLeft.get("Solution", this)),
     solR(cacheMapRight.get("Solution", this))
@@ -512,6 +517,7 @@ class dgPerfectGasLaw2d::riemannGodunov : public dataCacheDouble {
   dataCacheDouble &normals, &solL, &solR;
   public:
   riemannGodunov(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+    dataCacheDouble(cacheMapLeft),
     normals(cacheMapLeft.get("Normals", this)),
     solL(cacheMapLeft.get("Solution", this)),
     solR(cacheMapRight.get("Solution", this))
@@ -578,6 +584,7 @@ class dgPerfectGasLaw2d::maxConvectiveSpeed : public dataCacheDouble {
   dataCacheDouble &sol;
   public:
   maxConvectiveSpeed(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this))
   {
   };
@@ -599,6 +606,7 @@ class dgPerfectGasLaw2d::maxDiffusivity : public dataCacheDouble {
   dataCacheDouble &sol,&mu,&kappa;
   public:
   maxDiffusivity(dataCacheMap &cacheMap, const std::string &muFunctionName, const std::string &kappaFunctionName ):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this)),
     mu(cacheMap.get(muFunctionName,this)),
     kappa(cacheMap.get(kappaFunctionName,this))
@@ -662,6 +670,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
 // NON VISCOUS BOUNDARY FLUX
 //-------------------------------------------------------------------------------
     term(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
       sol(cacheMap.get("Solution",this)),
       normals(cacheMap.get("Normals",this)){}
     void _eval () { 
@@ -709,6 +718,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
     dataCacheDouble &sol;
     public:
     dirichletNonSlip(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this)){}
     void _eval () { 
       int nQP = sol().size1();
@@ -736,6 +746,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
     dataCacheDouble *diffusiveFlux;
   public:
     neumannNonSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
+      dataCacheDouble(cacheMap),
       _claw (claw),
       sol(cacheMap.get("Solution",this)),
       normals(cacheMap.get("Normals",this)){
@@ -771,6 +782,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
     dataCacheDouble &sol;
     public:
     dirichletSlip(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this)){}
     void _eval () { 
       int nQP = sol().size1();
@@ -795,6 +807,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
   public:
     neumannSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
       _claw (claw),
+      dataCacheDouble(cacheMap),
       sol(cacheMap.get("Solution",this)),
       normals(cacheMap.get("Normals",this)){
     }
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index b0fe3019c2208778b8f8d661ccc62765067d136c..7fef3c89e6f82dbe184a9f0156b46bf68701d2a9 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -7,6 +7,7 @@ class dgConservationLawShallowWater2d::advection: public dataCacheDouble {
   dataCacheDouble &sol;
   public:
   advection(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this))
   {};
   void _eval () { 
@@ -37,6 +38,7 @@ class dgConservationLawShallowWater2d::source: public dataCacheDouble {
   dataCacheDouble &xyz, &solution;
   public :
   source(dataCacheMap &cacheMap) : 
+    dataCacheDouble(cacheMap),
     xyz(cacheMap.get("XYZ",this)),
     solution(cacheMap.get("Solution",this))
   {}
@@ -61,6 +63,7 @@ class dgConservationLawShallowWater2d::riemann:public dataCacheDouble {
   dataCacheDouble &normals, &solL, &solR;
   public:
   riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+    dataCacheDouble(cacheMapLeft),
     normals(cacheMapLeft.get("Normals", this)),
     solL(cacheMapLeft.get("Solution", this)),
     solR(cacheMapRight.get("Solution", this))
@@ -99,6 +102,7 @@ class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition
     dataCacheDouble &sol,&normals;
     public:
     term(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this)),
     normals(cacheMap.get("Normals",this)){}
     void _eval () { 
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index 30d34e3a2ee835b27194c438534ff6042e4730b2..61f7fe076e6115df7629870e087e41f8f211447b 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -10,6 +10,7 @@ class dgConservationLawWaveEquation::hyperbolicFlux : public dataCacheDouble {
   const int _DIM,_nbf;    
   public:
   hyperbolicFlux(dataCacheMap &cacheMap,int DIM):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this)),_DIM(DIM),_nbf(DIM+1)
   {};
   void _eval () {
@@ -31,6 +32,7 @@ class dgConservationLawWaveEquation::maxConvectiveSpeed : public dataCacheDouble
   dataCacheDouble &sol;
   public:
   maxConvectiveSpeed(dataCacheMap &cacheMap):
+    dataCacheDouble(cacheMap),
     sol(cacheMap.get("Solution",this))
   {
   };
@@ -47,6 +49,7 @@ class dgConservationLawWaveEquation::riemann : public dataCacheDouble {
   const int _DIM,_nbf;
   public:
   riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, int DIM):
+    dataCacheDouble(cacheMapLeft),
     normals(cacheMapLeft.get("Normals", this)),
     solL(cacheMapLeft.get("Solution", this)),
     solR(cacheMapRight.get("Solution", this)),
@@ -107,6 +110,7 @@ class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition {
     dataCacheDouble &sol,&normals;    
     public:
     term(dataCacheMap &cacheMap, int DIM):
+      dataCacheDouble(cacheMap),
       sol(cacheMap.get("Solution",this)),
       normals(cacheMap.get("Normals",this)),
       _DIM(DIM){}
diff --git a/Solver/function.cpp b/Solver/function.cpp
index c1506b9bc256f63bc3bb3696ccd103275dbacbf6..00a19a8b9e155c50e51875e249d3e643f3b67b23 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -4,6 +4,9 @@
 #include <sstream>
 
 // dataCache members
+dataCache::dataCache(dataCacheMap *cacheMap) : _valid(false) {
+  cacheMap->addDataCache(this); //this dataCache can be deleted when the dataCacheMap is deleted
+}
 
 void dataCache::addMeAsDependencyOf (dataCache *newDep)
 {
@@ -45,8 +48,8 @@ function *function::get(std::string functionName, bool acceptNull)
 dataCacheElement &dataCacheMap::getElement(dataCache *caller) 
 {
   if(caller)
-    _cacheElement.addMeAsDependencyOf(caller);
-  return _cacheElement;
+    _cacheElement->addMeAsDependencyOf(caller);
+  return *_cacheElement;
 }
 
 dataCacheDouble &dataCacheMap::get(const std::string &functionName, dataCache *caller) 
@@ -64,15 +67,15 @@ dataCacheDouble &dataCacheMap::provideData(std::string name)
   dataCacheDouble *&r= _cacheDoubleMap[name];
   if(r!=NULL)
     throw;
-  r = new providedDataDouble;
+  r = new providedDataDouble(*this);
   return *r;
 }
 
 dataCacheMap::~dataCacheMap()
 {
-  for (std::map<std::string, dataCacheDouble*>::iterator it = _cacheDoubleMap.begin();
-      it!=_cacheDoubleMap.end(); it++) {
-    delete it->second;
+  for (std::set<dataCache*>::iterator it = _toDelete.begin();
+      it!=_toDelete.end(); it++) {
+    delete *it;
   }
 }
 
@@ -88,7 +91,7 @@ class functionXYZ : public function {
   int count;
    public:
     data(dataCacheMap *m) : 
-      dataCacheDouble(m->getNbEvaluationPoints(),3),
+      dataCacheDouble(*m, m->getNbEvaluationPoints(),3),
       _element(m->getElement(this)), _uvw(m->get("UVW", this))
     {
     }
@@ -118,7 +121,7 @@ class functionConstant::data : public dataCacheDouble {
  const functionConstant *_function;
  public:
  data(const functionConstant * function,dataCacheMap *m):
-   dataCacheDouble(m->getNbEvaluationPoints(),function->_source.size1()){
+   dataCacheDouble(*m,m->getNbEvaluationPoints(),function->_source.size1()){
      _function = function;
    }
  void _eval() {
diff --git a/Solver/function.h b/Solver/function.h
index 268bccd3734f226c5ce94058d17c70e7b8c82a8e..adeabe48737e6f3fe565b2c3cb6d8cd44c5ec6ea 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -43,7 +43,7 @@ protected :
       it!=_dependOnMe.end(); it++)
         (*it)->_valid=false;
   }
-  dataCache() : _valid(false) {}
+  dataCache(dataCacheMap *cacheMap);
   virtual ~dataCache(){};
 public :
   // dataCacheMap is the only one supposed to call this
@@ -105,8 +105,8 @@ class dataCacheDouble : public dataCache {
     }
     return _value;
   }
-  dataCacheDouble(){};
-  dataCacheDouble(int size1, int size2):_value(size1,size2){};
+  dataCacheDouble(dataCacheMap &map):dataCache(&map){};
+  dataCacheDouble(dataCacheMap &map,int size1, int size2):dataCache(&map),_value(size1,size2){};
   virtual ~dataCacheDouble(){};
 };
 
@@ -141,15 +141,17 @@ class dataCacheElement : public dataCache {
     _element=ele;
   };
   inline MElement *operator () () { return _element; }
+  dataCacheElement(dataCacheMap *map):dataCache(map){}
 };
 
 // more explanation at the head of this file
 class dataCacheMap {
+  friend class dataCache;
  private:
   int _nbEvaluationPoints;
   // keep track of the current element and all the dataCaches that
   // depend on it
-  dataCacheElement _cacheElement;
+  dataCacheElement *_cacheElement;
   std::map<std::string, dataCacheDouble*> _cacheDoubleMap;
   class providedDataDouble : public dataCacheDouble
   // for data provided by the algorithm and that does not have an _eval function
@@ -157,15 +159,23 @@ class dataCacheMap {
   {
     void _eval() {throw;};
     public:
-    providedDataDouble() {
+    providedDataDouble(dataCacheMap &map):dataCacheDouble(map) {
       _valid=true;
     }
   };
+  std::set<dataCache*> _toDelete;
+ protected:
+  void addDataCache(dataCache *data){
+    _toDelete.insert(data);
+  }
  public:
   dataCacheDouble &get(const std::string &functionName, dataCache *caller=0);
   dataCacheElement &getElement(dataCache *caller=0);
   dataCacheDouble &provideData(std::string name);
-  dataCacheMap(int nbEvaluationPoints):_nbEvaluationPoints(nbEvaluationPoints){}
+  dataCacheMap(int nbEvaluationPoints):_nbEvaluationPoints(nbEvaluationPoints){
+    _cacheElement= new dataCacheElement(this);
+}
+
   inline int getNbEvaluationPoints(){return _nbEvaluationPoints;}
   ~dataCacheMap();
 };
diff --git a/Solver/luaFunction.cpp b/Solver/luaFunction.cpp
index 2cf123532e98f3c3717e372f4a52e6cb1470b1fc..c0518cadef4233ad691305a6da001666e25240a3 100644
--- a/Solver/luaFunction.cpp
+++ b/Solver/luaFunction.cpp
@@ -14,7 +14,7 @@ class functionLua::data : public dataCacheDouble{
   std::vector<dataCacheDouble *> _dependencies;
   public:
   data(const functionLua *f, dataCacheMap *m):
-    dataCacheDouble(m->getNbEvaluationPoints(),f->_nbCol),
+    dataCacheDouble(*m,m->getNbEvaluationPoints(),f->_nbCol),
     _function(f)    
   {
     _dependencies.resize ( _function->_dependenciesName.size());