diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index be697b95cd7be8f09b2e70f4a4004ca5132fcfc8..d663c64837f26dcf7da8612aad288a357155ed67 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -12,66 +12,6 @@
 #include "dgTransformNodalValue.h"
 #include "meshPartition.h"
 
-// works for any number of groups 
-
-/*
-void dgAlgorithm::residualForSomeGroups( const dgConservationLaw &claw,
-          std::vector<dgGroupOfElements *>groups,
-          dgGroupCollection &groupCollection,
-          dgDofContainer &solution,
-          dgDofContainer &residu)
-{
-  solution.scatter();
-  int nbFields=claw.getNbFields();
-  //volume term
-  std::set<dgGroupOfFaces *>groupOfInternalFacesToUpdate;
-  std::vector<dgGroupOfFaces *>groupOfBoundaryFacesToUpdate;
-  for(size_t i=0;i<groups.size() ; i++) {
-    int groupId=groups[i]->getId();
-    residu.getGroupProxy(groupId).scale(0);
-    residualVolume(claw,*groups[i],solution.getGroupProxy(groupId),residu.getGroupProxy(groupId));
-    groupOfInternalFacesToUpdate.insert(groups[i]->getGroupsOfNeighboringFaces()->begin(),groups[i]->getGroupsOfNeighboringFaces()->end());
-    groupOfBoundaryFacesToUpdate.insert(groupOfBoundaryFacesToUpdate.end(),groups[i]->getGroupsOfBoundaryFaces()->begin(),groups[i]->getGroupsOfBoundaryFaces()->end());
-    groupOfInternalFacesToUpdate.insert(groups[i]->getGroupOfInsideFaces());
-  }
-  for(std::set<dgGroupOfFaces *>::iterator it=groupOfInternalFacesToUpdate.begin();it!=groupOfInternalFacesToUpdate.end();it++) {
-    dgGroupOfFaces *faces = *it;
-    int iGroupLeft = -1, iGroupRight = -1;
-    for(size_t j=0;j<groupCollection.getNbElementGroups(); j++) {
-      dgGroupOfElements *groupj = groupCollection.getElementGroup(j);
-      if (groupj == &faces->getGroupLeft()) iGroupLeft = j;
-      if (groupj == &faces->getGroupRight()) iGroupRight= j;
-    }
-    for(size_t j=0;j<groupCollection.getNbGhostGroups(); j++) {
-      dgGroupOfElements *groupj = groupCollection.getGhostGroup(j);
-      if (groupj == &faces->getGroupLeft()) iGroupLeft = j;
-      if (groupj == &faces->getGroupLeft())iGroupLeft = j + groupCollection.getNbElementGroups();
-      if (groupj == &faces->getGroupRight())iGroupRight= j + groupCollection.getNbElementGroups();
-    }
-    fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*nbFields);
-    fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*nbFields);
-    faces->mapToInterface(nbFields, solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight), solInterface);
-    residualInterface(claw,*faces,solInterface,solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight),residuInterface);
-    faces->mapFromInterface(nbFields, residuInterface,residu.getGroupProxy(iGroupLeft), residu.getGroupProxy(iGroupRight));
-  }
-  //boundaries
-  for(size_t i=0;i<groupOfBoundaryFacesToUpdate.size() ; i++) {
-    dgGroupOfFaces *faces = groupOfBoundaryFacesToUpdate[i];
-    int iGroupLeft = -1, iGroupRight = -1;
-    for(size_t j=0;j<groupCollection.getNbElementGroups(); j++) {
-      dgGroupOfElements *groupj = groupCollection.getElementGroup(j);
-      if (groupj == &faces->getGroupLeft())iGroupLeft = j;
-      if (groupj == &faces->getGroupRight())iGroupRight= j;
-    }
-    fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*nbFields);
-    fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*nbFields);
-    faces->mapToInterface(nbFields, solution.getGroupProxy(iGroupLeft), solution.getGroupProxy(iGroupRight), solInterface);
-    residualBoundary(claw,*faces,solInterface,solution.getGroupProxy(iGroupLeft),residuInterface);
-    faces->mapFromInterface(nbFields, residuInterface, residu.getGroupProxy(iGroupLeft), residu.getGroupProxy(iGroupRight));
-  }
-}
-*/
-
 void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
                                               const dgConservationLaw &claw,   // the conservation law
                                               const dgGroupOfElements & group, 
@@ -82,8 +22,8 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
   const int nbFields = claw.getNbFields();
   dataCacheMap cacheMap;
   cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
-  dataCacheDouble &sol = cacheMap.provideSolution(nbFields);
-  dataCacheDouble &UVW = cacheMap.provideParametricCoordinates();
+  dataCacheDouble &sol = cacheMap.get(function::getSolution(), NULL);
+  dataCacheDouble &UVW = cacheMap.get(function::getParametricCoordinates(), NULL);
   UVW.set()=group.getIntegrationPointsMatrix();
   // provided dataCache
   dataCacheDouble *maxConvectiveSpeed = claw.newMaxConvectiveSpeed(cacheMap);
@@ -98,8 +38,8 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
   double l_red_sq = l_red*l_red;
   DT.resize(group.getNbElements());
   for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
-    sol.set().setAsProxy(solution, iElement*nbFields, nbFields);
     cacheMap.setElement(group.getElement(iElement));
+    sol.set().setAsProxy(solution, iElement*nbFields, nbFields);
     const double L = group.getInnerRadius(iElement);
     double spectralRadius = 0.0;
     if (maximumDiffusivity){
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index 6d0837ed314f74a3a32dfb3f36832af9278e7ca7..511c721654c262a26a28dbb9131fcbcae18f085d 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/Solver/dgConservationLaw.cpp
@@ -11,7 +11,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
     public:
     term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const function *outsideValueFunction):
       dataCacheDouble(cacheMapLeft, 1, claw->getNbFields()),
-      solutionRight(cacheMapRight.provideSolution(claw->getNbFields())),
+      solutionRight(cacheMapRight.get(function::getSolution(),NULL)),
       outsideValue(cacheMapLeft.get(outsideValueFunction,this)),
       _claw(claw)
     {
@@ -51,7 +51,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
     public:
     maximumDiffusivity(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const function *outsideValueFunction):
       dataCacheDouble(cacheMapLeft, 1, claw->getNbFields()),
-      solutionRight(cacheMapRight.provideSolution(claw->getNbFields())),
+      solutionRight(cacheMapRight.get(function::getSolution(), NULL)),
       outsideValue(cacheMapLeft.get(outsideValueFunction,this)),
       _claw(claw)
     {
@@ -168,7 +168,7 @@ class dgBoundaryCondition::neumann_ : public dataCacheDouble {
 public:
   neumann_(dataCacheMap &cacheMap, dgConservationLaw *claw):
     dataCacheDouble(cacheMap,1,claw->getNbFields()),
-    normals(cacheMap.getNormals(this))
+    normals(cacheMap.get(function::getNormals(),this))
   {
     diffusiveFlux=claw->newDiffusiveFlux(cacheMap);
     if (diffusiveFlux)diffusiveFlux->addMeAsDependencyOf(this);
@@ -189,7 +189,7 @@ public:
 };
 
 dataCacheDouble *dgBoundaryCondition::newDiffusiveDirichletBC(dataCacheMap &cacheMapLeft) const {
-  return &cacheMapLeft.getSolution(NULL);
+  return &cacheMapLeft.get(function::getSolution(),NULL);
 }
 dataCacheDouble *dgBoundaryCondition::newDiffusiveNeumannBC(dataCacheMap &cacheMapLeft) const {
   return new dgBoundaryCondition::neumann_(cacheMapLeft, _claw);
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index f01104f9fcfe8ecb9502c12b697df94660c68369..7348a47f80912b679e0d2184576c6c5fa961b569 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -11,7 +11,7 @@ class dgConservationLawAdvectionDiffusion::advection : public dataCacheDouble {
   public:
   advection(const function *vFunction, dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,3),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(),this)),
     v(cacheMap.get(vFunction,this))
   {};
   void _eval () { 
@@ -44,9 +44,9 @@ class dgConservationLawAdvectionDiffusion::riemann : public dataCacheDouble {
   public:
   riemann(const function *vFunction, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
     dataCacheDouble(cacheMapLeft,1,2),
-    normals(cacheMapLeft.getNormals(this)),
-    solLeft(cacheMapLeft.getSolution(this)),
-    solRight(cacheMapRight.getSolution(this)),
+    normals(cacheMapLeft.get(function::getNormals(), this)),
+    solLeft(cacheMapLeft.get(function::getSolution(), this)),
+    solRight(cacheMapRight.get(function::getSolution(), this)),
     v(cacheMapLeft.get(vFunction,this))
   {};
   void _eval () { 
@@ -62,12 +62,13 @@ class dgConservationLawAdvectionDiffusion::riemann : public dataCacheDouble {
     }
   }
 };
+
 class dgConservationLawAdvectionDiffusion::diffusion : public dataCacheDouble {
   dataCacheDouble &solgrad, &nu;
   public:
   diffusion(const function * nuFunction, dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,3),
-    solgrad(cacheMap.getSolutionGradient(this)),
+    solgrad(cacheMap.get(function::getSolutionGradient(), this)),
     nu(cacheMap.get(nuFunction,this))
   {};
   void _eval () { 
@@ -78,6 +79,7 @@ class dgConservationLawAdvectionDiffusion::diffusion : public dataCacheDouble {
     }
   }
 };
+
 dataCacheDouble *dgConservationLawAdvectionDiffusion::newConvectiveFlux( dataCacheMap &cacheMap) const {
   return _vFunction ? new advection(_vFunction,cacheMap) : NULL;
 }
diff --git a/Solver/dgConservationLawMaxwell.cpp b/Solver/dgConservationLawMaxwell.cpp
index 49743b6da85efb5799895c9a0841676573c512c7..96da2f543397fad6c0ffadc8ce9738c21f67554a 100644
--- a/Solver/dgConservationLawMaxwell.cpp
+++ b/Solver/dgConservationLawMaxwell.cpp
@@ -8,7 +8,7 @@ class dgConservationLawMaxwell::advection : public dataCacheDouble {
   public:
   advection(const function *mu_epsFunction,dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,18),
-    sol(cacheMap.getSolution(this)), 
+    sol(cacheMap.get(function::getSolution(), this)), 
     mu_eps(cacheMap.get(mu_epsFunction,this))
   {};
   void _eval () { 
@@ -53,7 +53,7 @@ class dgConservationLawMaxwell::source: public dataCacheDouble {
   public :
   source(dataCacheMap &cacheMap) : 
     dataCacheDouble(cacheMap,1,6),
-    solution(cacheMap.getSolution(this))
+    solution(cacheMap.get(function::getSolution(), this))
   {}
   void _eval () {
     int nQP = _value.size1();
@@ -75,9 +75,9 @@ class dgConservationLawMaxwell::riemann:public dataCacheDouble {
   public:
   riemann(const function *mu_epsFunction,dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
     dataCacheDouble(cacheMapLeft,1,12),
-    normals(cacheMapLeft.getNormals( this)),
-    solL(cacheMapLeft.getSolution( this)),
-    solR(cacheMapRight.getSolution( this)),
+    normals(cacheMapLeft.get(function::getNormals(), this)),
+    solL(cacheMapLeft.get(function::getSolution(), this)),
+    solR(cacheMapRight.get(function::getSolution(), this)),
     mu_eps(cacheMapLeft.get(mu_epsFunction,this))
   {};
   void _eval () { 
@@ -179,8 +179,8 @@ class dgBoundaryConditionMaxwellWall : public dgBoundaryCondition {
     public:
     term(dataCacheMap &cacheMap):
       dataCacheDouble(cacheMap,1,6),
-      sol(cacheMap.getSolution(this)),
-      normals(cacheMap.getNormals(this))
+      sol(cacheMap.get(function::getSolution(), this)),
+      normals(cacheMap.get(function::getNormals(), this))
       {}
     void _eval () { 
       int nQP = sol().size1();
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index b6fc735e69acb626f8b66ee898b4725b4f803c63..2ec3831ca8eea0b6643aea11142399c912fe39a5 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -365,7 +365,7 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
   public:
   advection(dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,8),
-    sol(cacheMap.getSolution(this))
+    sol(cacheMap.get(function::getSolution(), this))
   {};
   void _eval () { 
     const int nQP = _value.size1();      
@@ -401,8 +401,8 @@ class dgPerfectGasLaw2d::diffusion : public dataCacheDouble {
             const function * muFunction,
             const function * kappaFunction ):
     dataCacheDouble(cacheMap,1,12),
-    sol(cacheMap.getSolution(this)),
-    grad(cacheMap.getSolutionGradient(this)),
+    sol(cacheMap.get(function::getSolution(), this)),
+    grad(cacheMap.get(function::getSolutionGradient(), this)),
     mu(cacheMap.get(muFunction,this)),
     kappa(cacheMap.get(kappaFunction,this))
   {};
@@ -459,7 +459,7 @@ class dgPerfectGasLaw2d::source : public dataCacheDouble {
 public:
   source(dataCacheMap &cacheMap, const function *sourceFunction):
     dataCacheDouble(cacheMap,1,4),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(), this)),
     s(cacheMap.get(sourceFunction,this))
   {};
   void _eval () { 
@@ -479,7 +479,7 @@ class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
 public:
   clipToPhysics(dataCacheMap &cacheMap, double presMin, double rhoMin):
     dataCacheDouble(cacheMap,1,4),
-    sol(cacheMap.getSolution(this))
+    sol(cacheMap.get(function::getSolution(), this))
   {
     _presMin=presMin;
     _rhoMin=rhoMin;
@@ -512,9 +512,9 @@ class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
   public:
   riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
     dataCacheDouble(cacheMapLeft,1,12),
-    normals(cacheMapLeft.getNormals( this)),
-    solL(cacheMapLeft.getSolution( this)),
-    solR(cacheMapRight.getSolution( this))
+    normals(cacheMapLeft.get(function::getNormals(), this)),
+    solL(cacheMapLeft.get(function::getSolution(), this)),
+    solR(cacheMapRight.get(function::getSolution(), this))
   {};
   void _eval () { 
     int nQP = solL().size1();
@@ -543,9 +543,9 @@ class dgPerfectGasLaw2d::riemannGodunov : public dataCacheDouble {
   public:
   riemannGodunov(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
     dataCacheDouble(cacheMapLeft,1,12),
-    normals(cacheMapLeft.getNormals( this)),
-    solL(cacheMapLeft.getSolution( this)),
-    solR(cacheMapRight.getSolution( this))
+    normals(cacheMapLeft.get(function::getNormals(), this)),
+    solL(cacheMapLeft.get(function::getSolution(), this)),
+    solR(cacheMapRight.get(function::getSolution(), this))
   {};
   void _eval () { 
     int nQP = solL().size1();
@@ -607,7 +607,7 @@ class dgPerfectGasLaw2d::maxConvectiveSpeed : public dataCacheDouble {
   public:
   maxConvectiveSpeed(dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,1),
-    sol(cacheMap.getSolution(this))
+    sol(cacheMap.get(function::getSolution(), this))
   {
   };
   void _eval () {
@@ -627,7 +627,7 @@ class dgPerfectGasLaw2d::maxDiffusivity : public dataCacheDouble {
   public:
   maxDiffusivity(dataCacheMap &cacheMap, const function *muFunction, const function *kappaFunction ):
     dataCacheDouble(cacheMap,1,1),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(), this)),
     mu(cacheMap.get(muFunction,this)),
     kappa(cacheMap.get(kappaFunction,this))
   {
@@ -687,8 +687,8 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
 //-------------------------------------------------------------------------------
     term(dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,4),
-      sol(cacheMap.getSolution(this)),
-      normals(cacheMap.getNormals(this)){}
+      sol(cacheMap.get(function::getSolution(), this)),
+      normals(cacheMap.get(function::getNormals(), this)){}
     void _eval () { 
       int nQP = sol().size1();
       for(int i=0; i< nQP; i++) {
@@ -732,7 +732,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
     public:
     dirichletNonSlip(dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,4),
-    sol(cacheMap.getSolution(this)){}
+    sol(cacheMap.get(function::getSolution(), this)){}
     void _eval () { 
       int nQP = sol().size1();
       for(int i=0; i< nQP; i++) {
@@ -758,8 +758,8 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
     neumannNonSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
       dataCacheDouble(cacheMap,1,4),
       _claw (claw),
-      sol(cacheMap.getSolution(this)),
-      normals(cacheMap.getNormals(this)){
+      sol(cacheMap.get(function::getSolution(), this)),
+      normals(cacheMap.get(function::getNormals(), this)){
       diffusiveFlux=_claw->newDiffusiveFlux(cacheMap);
       if (diffusiveFlux)diffusiveFlux->addMeAsDependencyOf(this);
     }
@@ -789,8 +789,8 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
     dataCacheDouble &sol;
     public:
     dirichletSlip(dataCacheMap &cacheMap):
-    dataCacheDouble(cacheMap,1,cacheMap.getSolution(NULL)().size2()),
-    sol(cacheMap.getSolution(this)){}
+    dataCacheDouble(cacheMap,1,cacheMap.get(function::getSolution(), NULL)().size2()),
+    sol(cacheMap.get(function::getSolution(), this)){}
     void _eval () { 
       int nQP = sol().size1();
       for(int i=0; i< nQP; i++) {
@@ -812,8 +812,8 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
     neumannSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
       _claw (claw),
       dataCacheDouble(cacheMap,1,4),
-      sol(cacheMap.getSolution(this)),
-      normals(cacheMap.getNormals(this)){
+      sol(cacheMap.get(function::getSolution(), this)),
+      normals(cacheMap.get(function::getNormals(), this)){
     }
     void _eval () { 
       int nQP = sol().size1();
diff --git a/Solver/dgConservationLawShallowWater1d.cpp b/Solver/dgConservationLawShallowWater1d.cpp
index f9d987f4de5b2fa54d56489531df4e796489dfda..a998cda842bfcaa27a114251e45f933ace0e5aba 100644
--- a/Solver/dgConservationLawShallowWater1d.cpp
+++ b/Solver/dgConservationLawShallowWater1d.cpp
@@ -39,7 +39,7 @@ class dgConservationLawShallowWater1d::clipToPhysics : public dataCacheDouble {
 public:
   clipToPhysics(dataCacheMap &cacheMap, const function *bathymetry, double hMin):
     dataCacheDouble(cacheMap,1,3),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(),this)),
     _bathymetry(cacheMap.get(bathymetry,this))
   {
     _hMin=hMin;
@@ -63,7 +63,7 @@ class dgConservationLawShallowWater1d::maxConvectiveSpeed : public dataCacheDoub
   public:
   maxConvectiveSpeed(dataCacheMap &cacheMap, const function *celerity):
     dataCacheDouble(cacheMap,1,1),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(),this)),
     _celerity(cacheMap.get(celerity,this))
   {};
   void _eval () {
@@ -79,7 +79,7 @@ class dgConservationLawShallowWater1d::advection: public dataCacheDouble {
   public:
   advection(dataCacheMap &cacheMap, const function *bathymetry, const function *pressure):
     dataCacheDouble(cacheMap,1,6),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(),this)),
     _bathymetry(cacheMap.get(bathymetry,this)), 
     _pressure(cacheMap.get(pressure,this))
   {};
@@ -110,7 +110,7 @@ class dgConservationLawShallowWater1d::source: public dataCacheDouble {
   source(dataCacheMap &cacheMap,  const function *linearDissipation, const function *bathymetry) : 
     dataCacheDouble(cacheMap,1,2),
     _linearDissipation(cacheMap.get(linearDissipation,this)),
-    solution(cacheMap.getSolution(this)),
+    solution(cacheMap.get(function::getSolution(),this)),
     _bathymetry(cacheMap.get(bathymetry,this))
   {}
   void _eval () {
@@ -131,9 +131,9 @@ class dgConservationLawShallowWater1d::riemann:public dataCacheDouble {
   riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, const function *bathymetry, 
           const function *pressure, const function *celerity):
     dataCacheDouble(cacheMapLeft,1,4),
-    normals(cacheMapLeft.getNormals(this)),
-    solL(cacheMapLeft.getSolution(this)),
-    solR(cacheMapRight.getSolution(this)),
+    normals(cacheMapLeft.get(function::getNormals(), this)),
+    solL(cacheMapLeft.get(function::getSolution(),this)),
+    solR(cacheMapRight.get(function::getSolution(),this)),
     _bathymetryL(cacheMapLeft.get(bathymetry, this)), 
     _bathymetryR(cacheMapRight.get(bathymetry, this)), 
     _pressureL(cacheMapLeft.get(pressure, this)), 
@@ -222,9 +222,9 @@ class dgConservationLawShallowWater1d::boundaryWall : public dgBoundaryCondition
     public:
     term(dataCacheMap &cacheMap, const function *pressure):
     dataCacheDouble(cacheMap,1,2),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(),this)),
     _pressure(cacheMap.get(pressure,this)),
-    normals(cacheMap.getNormals(this))
+    normals(cacheMap.get(function::getNormals(), this))
     {}
     void _eval () { 
       int nQP = sol().size1();
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index 476cc50bd3ba440286067095eb252537751b12d8..46ac472ed4598410c267dcc457fb340e84049c38 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -10,6 +10,7 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
   class boundaryWall;
   class maxConvectiveSpeed;
   const function *_linearDissipation, *_quadraticDissipation, *_source, *_coriolisFactor, *_bathymetry;
+  mutable function *_advection, *_riemann;
   public:
   dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const;
   dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
@@ -22,17 +23,38 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
   inline void setQuadraticDissipation(const function *quadraticDissipation){ _quadraticDissipation = quadraticDissipation;}
   inline void setSource(const function *source){ _source = source;}
   inline void setBathymetry(const function *bathymetry) { _bathymetry = bathymetry;}
-  dgConservationLawShallowWater2d() 
+  dgConservationLawShallowWater2d();
+  dgBoundaryCondition *newBoundaryWall();
+};
+
+class dgConservationLawShallowWater2d::advection: public function {
+  public:
+  advection(const function *bathymetry, const function *solution):function(9)
   {
-    _nbf = 3; // eta u v
-    function *fzero = functionConstantNew(0.);
-    _bathymetry = fzero;
-    _linearDissipation = fzero;
-    _coriolisFactor = fzero;
-    _quadraticDissipation = fzero;
-    _source = fzero;
+    addArgument(bathymetry);
+    addArgument(solution);
+  }
+  void call (dataCacheMap *m, const fullMatrix<double> &bathymetry, const fullMatrix<double> &sol, fullMatrix<double> &value) { 
+    int nQP = value.size1();
+    for(int i=0; i< nQP; i++) {
+      double h = bathymetry(i,0);
+      double eta = sol(i,0);
+      double u = sol(i,1);
+      double v = sol(i,2);
+      // flux_x
+      value(i,0) = (h+eta)*u;
+      value(i,1) = g*eta+u*u;
+      value(i,2) = u*v;
+      // flux_y
+      value(i,3) = (h+eta)*v;
+      value(i,4) = v*u;
+      value(i,5) = g*eta+v*v;
+      // flux_z
+      value(i,6) = 0;
+      value(i,7) = 0;
+      value(i,8) = 0;
+    }
   }
-  dgBoundaryCondition *newBoundaryWall();
 };
 
 class dgConservationLawShallowWater2d::clipToPhysics : public dataCacheDouble {
@@ -41,7 +63,7 @@ class dgConservationLawShallowWater2d::clipToPhysics : public dataCacheDouble {
 public:
   clipToPhysics(dataCacheMap &cacheMap, const function *bathymetry, double hMin):
     dataCacheDouble(cacheMap,1,3),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(), this)),
     _bathymetry(cacheMap.get(bathymetry,this))
   {
     _hMin=hMin;
@@ -66,7 +88,7 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public dataCacheDoubl
   public:
   maxConvectiveSpeed(dataCacheMap &cacheMap, const function *bathymetry):
     dataCacheDouble(cacheMap,1,1),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(), this)),
     _bathymetry(cacheMap.get(bathymetry,this))
   {};
   void _eval () { 
@@ -81,12 +103,13 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public dataCacheDoubl
   }
 };
 
-class dgConservationLawShallowWater2d::advection: public dataCacheDouble {
+
+/*class dgConservationLawShallowWater2d::advection: public dataCacheDouble {
   dataCacheDouble &sol, &_bathymetry;
   public:
   advection(dataCacheMap &cacheMap, const function *bathymetry):
     dataCacheDouble(cacheMap,1,9),
-    sol(cacheMap.getSolution(this)),
+    sol(cacheMap.get(function::getSolution(), this)),
     _bathymetry(cacheMap.get(bathymetry,this))
   {};
   void _eval () { 
@@ -110,7 +133,7 @@ class dgConservationLawShallowWater2d::advection: public dataCacheDouble {
       _value(i,8) = 0;
     }
   }
-};
+};*/
 
 class dgConservationLawShallowWater2d::source: public dataCacheDouble {
   dataCacheDouble &solution,&solutionGradient;
@@ -122,8 +145,8 @@ class dgConservationLawShallowWater2d::source: public dataCacheDouble {
     _source(cacheMap.get(source,this)),
     _linearDissipation(cacheMap.get(linearDissipation,this)),
     _quadraticDissipation(cacheMap.get(quadraticDissipation,this)),
-    solution(cacheMap.getSolution(this)),
-    solutionGradient(cacheMap.getSolutionGradient(this))
+    solution(cacheMap.get(function::getSolution(), this)),
+    solutionGradient(cacheMap.get(function::getSolutionGradient(), this))
   {
   }
   void _eval () {
@@ -143,14 +166,87 @@ class dgConservationLawShallowWater2d::source: public dataCacheDouble {
     }
   }
 };
+
+#if 0
+class dgConservationLawShallowWater2d::riemann:public function {
+  public:
+  riemann (const function *bathymetry): function(6)
+  {
+    addArgument (function::getNormals(), 0);
+    addArgument (function::getSolution(), 0);
+    addArgument (function::getSolution(), 1);
+    addArgument (bathymetry, 0);
+  };
+  void call  (dataCacheMap *m, const fullMatrix<double> &normals, const fullMatrix<double> &solL, const fullMatrix<double> &solR, const fullMatrix<double> &_bathymetry, fullMatrix<double> &value) { 
+    int nQP = solL.size1();
+    for(int i=0; i< nQP; i++) {
+      double h = _bathymetry(i,0);
+      double nx = normals(0,i);
+      double ny = normals(1,i);
+      double HR = solR(i,0) + h, HL = solL(i,0) + h;
+      double uL = nx*solL(i,1) + ny*solL(i,2), uR = nx*solR(i,1) + ny*solR(i,2);
+      double vL = -ny*solL(i,1) + nx*solL(i,2), vR = - ny*solR(i,1) + nx*solR(i,2);
+      double HuL = uL*HL, HuR = uR*HR;
+      double HvL = vL*HL, HvR = vR*HR;
+      double HM = (HL+HR)/2, HJ = (HL-HR)/2;
+      double HuM = (HuL+HuR)/2, HuJ = (HuL-HuR)/2;
+      double HvM = (HvL+HvR)/2, HvJ = (HvL-HvR)/2;
+      double sqHL = sqrt(HL), sqHR = sqrt(HR);
+      double u_roe = (sqHL*uL + sqHR*uR) / (sqHL + sqHR);
+      double v_roe = (sqHL*vL + sqHR*vR) / (sqHL + sqHR);
+      double c_roe = sqrt(g*HM);
+      double H = HM + (HuJ - u_roe *HJ) / c_roe;
+      double Hu = HuM + (c_roe - u_roe*u_roe/c_roe) *HJ + u_roe/c_roe *HuJ;
+      double Hv = -v_roe*u_roe/c_roe*HJ + v_roe/c_roe*HuJ + (u_roe>0 ? -v_roe*HJ+HvL : v_roe*HJ+HvR);
+      double u = Hu / H;
+      double v = Hv / H;
+      double eta  = H-h;
+
+      /*double u = (uL+uR)/2; // 
+      double v = (vL+vR)/2;
+      double eta = (HL+HR)/2-h;
+      double c = hypot(u,v)+sqrt(g*(h+eta));
+      double Fun = -g*eta - u*u - (uL-uR)/2 * c;
+      double Fut = -u*( u>0 ? vL:vR);
+      double Feta = -(h+eta)*u - (HL-HR)/2 * c;*/
+
+      double Fun = -g*eta - u*u;
+      double Fut = -u*v;
+      double Feta = -(h+eta)*u;
+
+      /* //linear equations 
+      double unL = nx*solL(i,1) + ny*solL(i,2);
+      double unR = nx*solR(i,1) + ny*solR(i,2);
+      double etaR = solR(i,0);
+      double etaL = solL(i,0);
+      double sq_g_h = sqrt(g/h);
+      double etaRiemann = (etaL+etaR + (unL-unR)/sq_g_h)/2;
+      double unRiemann = (unL+unR + (etaL-etaR)*sq_g_h)/2;
+      double Fun = -g*etaRiemann;
+      double Fut = 0;
+      double Feta = -h*unRiemann;
+      */
+
+      value(i,0) = Feta;
+      value(i,1) = Fun*nx - Fut*ny;
+      value(i,2) = Fun*ny + Fut*nx;
+      value(i,3) = -value(i,0);
+      value(i,4) = -value(i,1);
+      value(i,5) = -value(i,2);
+    }
+  }
+};
+#endif
+
+#if 1
 class dgConservationLawShallowWater2d::riemann:public dataCacheDouble {
   dataCacheDouble &normals, &solL, &solR, &_bathymetry;
   public:
   riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, const function *bathymetry):
     dataCacheDouble(cacheMapLeft,1,6),
-    normals(cacheMapLeft.getNormals( this)),
-    solL(cacheMapLeft.getSolution( this)),
-    solR(cacheMapRight.getSolution( this)),
+    normals(cacheMapLeft.get(function::getNormals(), this)),
+    solL(cacheMapLeft.get(function::getSolution(), this)),
+    solR(cacheMapRight.get(function::getSolution(), this)),
     _bathymetry(cacheMapLeft.get(bathymetry, this))
   {};
   void _eval () { 
@@ -212,6 +308,7 @@ class dgConservationLawShallowWater2d::riemann:public dataCacheDouble {
     }
   }
 };
+#endif
 
 class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition {
   class term : public dataCacheDouble {
@@ -219,8 +316,8 @@ class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition
     public:
     term(dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,3),
-    sol(cacheMap.getSolution(this)),
-    normals(cacheMap.getNormals(this)){}
+    sol(cacheMap.get(function::getSolution(), this)),
+    normals(cacheMap.get(function::getNormals(), this)){}
     void _eval () { 
       int nQP = sol().size1();
       for(int i=0; i< nQP; i++) {
@@ -242,9 +339,15 @@ class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition
 };
 
 dataCacheDouble *dgConservationLawShallowWater2d::newConvectiveFlux( dataCacheMap &cacheMap) const {
-  return new advection(cacheMap, _bathymetry);
+  if(_advection==NULL)
+    _advection = new advection(_bathymetry, function::getSolution());
+  return &cacheMap.get(_advection, NULL);
 }
+
 dataCacheDouble *dgConservationLawShallowWater2d::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
+  //if(_riemann==NULL)
+    //_riemann = new riemann(_bathymetry);
+  //return &cacheMapLeft.get(_riemann, NULL);
   return new riemann(cacheMapLeft, cacheMapRight, _bathymetry);
 }
 dataCacheDouble *dgConservationLawShallowWater2d::newDiffusiveFlux( dataCacheMap &cacheMap) const {
@@ -264,6 +367,17 @@ dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(){
 dataCacheDouble *dgConservationLawShallowWater2d::newClipToPhysics( dataCacheMap &cacheMap) const {
   return new clipToPhysics(cacheMap, _bathymetry, 1e-5);
 }
+dgConservationLawShallowWater2d::dgConservationLawShallowWater2d()
+{
+  _nbf = 3; // eta u v
+  function *fzero = functionConstantNew(0.);
+  _bathymetry = fzero;
+  _linearDissipation = fzero;
+  _coriolisFactor = fzero;
+  _quadraticDissipation = fzero;
+  _source = fzero;
+  _advection = NULL;
+}
 
 #include "Bindings.h"
 void dgConservationLawShallowWater2dRegisterBindings (binding *b){
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index dfaf0553d75ba8347f4aad13c1e9125d32faeaa7..419d0ef3135053c13ea878720d6aaf7130609389 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -12,7 +12,7 @@ class dgConservationLawWaveEquation::hyperbolicFlux : public dataCacheDouble {
   public:
   hyperbolicFlux(dataCacheMap &cacheMap,int DIM):
     dataCacheDouble(cacheMap,1,3*(DIM+1)),
-    sol(cacheMap.getSolution(this)),_DIM(DIM),_nbf(DIM+1)
+    sol(cacheMap.get(function::getSolution(), this)),_DIM(DIM),_nbf(DIM+1)
   {};
   void _eval () {
     int nQP = sol().size1();
@@ -32,7 +32,7 @@ class dgConservationLawWaveEquation::maxConvectiveSpeed : public dataCacheDouble
   public:
   maxConvectiveSpeed(dataCacheMap &cacheMap):
     dataCacheDouble(cacheMap,1,1),
-    sol(cacheMap.getSolution(this))
+    sol(cacheMap.get(function::getSolution(), this))
   {
   };
   void _eval () {
@@ -46,9 +46,9 @@ class dgConservationLawWaveEquation::riemann : public dataCacheDouble {
   public:
   riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, int DIM):
     dataCacheDouble(cacheMapLeft,1,2*(DIM+1)),
-    normals(cacheMapLeft.getNormals( this)),
-    solL(cacheMapLeft.getSolution( this)),
-    solR(cacheMapRight.getSolution( this)),
+    normals(cacheMapLeft.get(function::getNormals(), this)),
+    solL(cacheMapLeft.get(function::getSolution(), this)),
+    solR(cacheMapRight.get(function::getSolution(), this)),
     _DIM(DIM),_nbf(DIM+1)
   {};
   void _eval () { 
@@ -105,8 +105,8 @@ class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition {
     public:
     term(dataCacheMap &cacheMap, int DIM):
       dataCacheDouble(cacheMap,1,DIM+1),
-      sol(cacheMap.getSolution(this)),
-      normals(cacheMap.getNormals(this)),
+      sol(cacheMap.get(function::getSolution(), this)),
+      normals(cacheMap.get(function::getNormals(), this)),
       _DIM(DIM){}
     void _eval () { 
       int nQP = sol().size1();
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 61b0d671ba23558bc879fa0dcb061fdb24261900..cef75b4a78450d9c88e6260a66c15a4118b6b6ff 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -229,7 +229,7 @@ void dgDofContainer::L2Projection(const function *f){
     fullMatrix<double> Source = fullMatrix<double> (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
     dataCacheMap cacheMap;
     cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
-    cacheMap.provideParametricCoordinates().set()=group.getIntegrationPointsMatrix();
+    cacheMap.get(function::getParametricCoordinates(), NULL).set()=group.getIntegrationPointsMatrix();
     dataCacheDouble &sourceTerm = cacheMap.get(f);
     fullMatrix<double> source;
     for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
diff --git a/Solver/dgFunctionIntegrator.cpp b/Solver/dgFunctionIntegrator.cpp
index 71b7559b91c1a66f27860a78ad5aa42122ad5d23..d72984bae92a536afb1d25f8276df033b1028dd2 100644
--- a/Solver/dgFunctionIntegrator.cpp
+++ b/Solver/dgFunctionIntegrator.cpp
@@ -11,8 +11,8 @@ dgFunctionIntegrator::dgFunctionIntegrator(const function *f):_f(f){}
 void dgFunctionIntegrator::compute(dgDofContainer *sol,fullMatrix<double> &result){
   int nbFields=sol->getNbFields();
   dataCacheMap cacheMap;
-  dataCacheDouble &UVW=cacheMap.provideParametricCoordinates();
-  dataCacheDouble &solutionQPe=cacheMap.provideSolution(nbFields);
+  dataCacheDouble &UVW=cacheMap.get(function::getParametricCoordinates(),NULL);
+  dataCacheDouble &solutionQPe=cacheMap.get(function::getSolution(),NULL);
   dataCacheDouble &F=cacheMap.get(_f);
   int nbRowResult=result.size1();
   result.scale(0.0);
diff --git a/Solver/dgLimiter.cpp b/Solver/dgLimiter.cpp
index c8cb8461fe0b69948458f204dd9a0dfff787ced7..8d4bba90ed13231df12667eade34953278bc7950 100644
--- a/Solver/dgLimiter.cpp
+++ b/Solver/dgLimiter.cpp
@@ -114,7 +114,7 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution)
 
     dataCacheMap cacheMap;
     cacheMap.setNbEvaluationPoints(egroup->getNbNodes());//nbdofs for each element
-    dataCacheDouble &solutionE = cacheMap.provideSolution(nbFields);
+    dataCacheDouble &solutionE = cacheMap.get(function::getSolution(), NULL);
     dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
     if (solutionEClipped){
       for (int iElement=0 ; iElement<egroup->getNbElements() ;++iElement) {
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index 71d15e8ae228329ace8d7f9f9303fb12a9b76d7a..36bc083ecadab9981e073234f89d6950c36db382 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -10,9 +10,9 @@ dgResidualVolume::dgResidualVolume(const dgConservationLaw &claw):
   _cacheMap(new dataCacheMap),
   _claw(claw),
   _nbFields(_claw.getNbFields()),
-  _UVW(_cacheMap->provideParametricCoordinates()),
-  _solutionQPe(_cacheMap->provideSolution(_nbFields)),
-  _gradientSolutionQPe(_cacheMap->provideSolutionGradient(_nbFields)),
+  _UVW(_cacheMap->get(function::getParametricCoordinates(), NULL)),
+  _solutionQPe(_cacheMap->get(function::getSolution(), NULL)),
+  _gradientSolutionQPe(_cacheMap->get(function::getSolutionGradient(),NULL)),
   _sourceTerm(_claw.newSourceTerm(*_cacheMap)),
   _convectiveFlux(_claw.newConvectiveFlux(*_cacheMap)),
   _diffusiveFlux(_claw.newDiffusiveFlux(*_cacheMap))
@@ -50,18 +50,18 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
   fullMatrix<double> fuvwe;
   fullMatrix<double> source;
   fullMatrix<double> dPsiDx,dofs, gradSolProxy; 
+  _gradientSolutionQPe.set().resize(group.getNbIntegrationPoints(),3*_nbFields);
   // ----- 2.3 --- iterate on elements
   for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
       // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
+    _cacheMap->setElement(group.getElement(iElement));
     _solutionQPe.set().setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
-
     if(_gradientSolutionQPe.somethingDependOnMe()){
       dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
       dofs.setAsProxy(solution, _nbFields*iElement, _nbFields);
       gradSolProxy.setAsShapeProxy(_gradientSolutionQPe.set(),group.getNbIntegrationPoints()*3, _nbFields);
       dPsiDx.mult(dofs, gradSolProxy);
     }
-    _cacheMap->setElement(group.getElement(iElement));
     if(_convectiveFlux || _diffusiveFlux) {
       // ----- 2.3.3 --- compute fluxes in UVW coordinates
       for (int iUVW=0;iUVW<group.getDimUVW();iUVW++) {
@@ -281,16 +281,18 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
   std::vector<dataCacheMap> caches(nbConnections);
   std::vector<const dgGroupOfConnections*> connections(nbConnections);
   std::vector<const dgGroupOfElements*> elementGroups(nbConnections);
-  std::vector<dataCacheDouble*> diffusiveFluxes(nbConnections), maximumDiffusivities(nbConnections);
+  std::vector<dataCacheDouble*> diffusiveFluxes(nbConnections), maximumDiffusivities(nbConnections), solutionGradients(nbConnections),
+    normals(nbConnections), solutions(nbConnections), parametricCoordinates(nbConnections);
   for (int i=0; i<nbConnections; i++) {
     connections[i] = &group.getGroupOfConnections(i);
     elementGroups[i] = &connections[i]->getGroupOfElements();
-    caches[i].provideParametricCoordinates();
-    caches[i].provideSolution(_nbFields);
-    caches[i].provideSolutionGradient(_nbFields);
-    caches[i].provideNormals();
     diffusiveFluxes[i] =_claw.newDiffusiveFlux(caches[i]);
     maximumDiffusivities[i] = _claw.newMaximumDiffusivity(caches[i]); 
+    parametricCoordinates[i] = &caches[i].get(function::getParametricCoordinates(),NULL);
+    normals[i] = &caches[i].get(function::getNormals(), NULL);
+    solutions[i] = &caches[i].get(function::getSolution(), NULL);
+    solutionGradients[i] = &caches[i].get(function::getSolutionGradient(), NULL);
+    solutionGradients[i]->set().resize(group.getNbIntegrationPoints(), 3*_nbFields);
     caches[i].setNbEvaluationPoints(group.getNbIntegrationPoints());
   }
   dataCacheDouble *riemannSolver = NULL, *neumann = NULL, *dirichlet = NULL, *maximumDiffusivityOut = NULL;
@@ -321,14 +323,14 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
     for (int i=0; i<nbConnections; i++) {
       // B1 )  adjust the proxies for this element
       caches[i].setElement(connections[i]->getElement(iFace));
-      caches[i].getParametricCoordinates(NULL).set() = connections[i]->getIntegrationPointsOnElement(iFace);
-      caches[i].getSolution(NULL).set().setAsProxy(solutionQP, (iFace*nbConnections+i)*_nbFields, _nbFields);
-      caches[i].getNormals(NULL).set().setAsProxy(connections[i]->getNormals(), iFace*group.getNbIntegrationPoints(), group.getNbIntegrationPoints());
+      parametricCoordinates[i]->set() = connections[i]->getIntegrationPointsOnElement(iFace);
+      solutions[i]->set().setAsProxy(solutionQP, (iFace*nbConnections+i)*_nbFields, _nbFields);
+      normals[i]->set().setAsProxy(connections[i]->getNormals(), iFace*group.getNbIntegrationPoints(), group.getNbIntegrationPoints());
       // B2 ) compute the gradient of the solution
-      if(caches[i].getSolutionGradient(NULL).somethingDependOnMe()) {
+      if(solutionGradients[i]->somethingDependOnMe()) {
         dofs.setAsProxy(*solutionOnElements[i], _nbFields*connections[i]->getElementId(iFace), _nbFields);
         dPsiDx.setAsProxy(connections[i]->getDPsiDx(),iFace*elementGroups[i]->getNbNodes(),elementGroups[i]->getNbNodes());
-        gradSolProxy.setAsShapeProxy(caches[i].getSolutionGradient(NULL).set(),group.getNbIntegrationPoints()*3, _nbFields);
+        gradSolProxy.setAsShapeProxy(solutionGradients[i]->set(),group.getNbIntegrationPoints()*3, _nbFields);
         dPsiDx.mult(dofs, gradSolProxy);
       }
     }
@@ -340,7 +342,7 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
         //just for the lisibility :
         const fullMatrix<double> &dfl = (*diffusiveFluxes[0])();
         const fullMatrix<double> &dfr = (*diffusiveFluxes[1])();
-        const fullMatrix<double> &n = (caches[0].getNormals(NULL))();
+        const fullMatrix<double> &n = (*normals[0])();
         for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
           const double detJ = group.getDetJ (iFace, iPt);
           for (int k=0;k<_nbFields;k++) { 
@@ -353,7 +355,7 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
                 )/group.getInterfaceSurface(iFace);
             double nu = std::max((*maximumDiffusivities[1])(iPt,0),(*maximumDiffusivities[0])(iPt,0));
             double mu = (p+1)*(p+dim)/dim*nu/minl;
-            double solutionJumpPenalty = (caches[0].getSolution(NULL)(iPt,k)-caches[1].getSolution(NULL)(iPt,k))/2*mu;
+            double solutionJumpPenalty = ((*solutions[0])(iPt,k)-(*solutions[1])(iPt,k))/2*mu;
             normalFluxQP(iPt,k) -= (meanNormalFlux+solutionJumpPenalty)*detJ;
             normalFluxQP(iPt,k+_nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ;
           }
@@ -365,7 +367,7 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
             double minl = elementGroups[0]->getElementVolume(connections[0]->getElementId(iFace))/group.getInterfaceSurface(iFace);
             double nu = std::max((*maximumDiffusivities[0])(iPt,0),(*maximumDiffusivityOut)(iPt,0));
             double mu = (p+1)*(p+dim)/dim*nu/minl;
-            double solutionJumpPenalty = (caches[0].getSolution(NULL)(iPt,k)-(*dirichlet)(iPt,k))/2*mu;
+            double solutionJumpPenalty = ((*solutions[0])(iPt,k)-(*dirichlet)(iPt,k))/2*mu;
             normalFluxQP(iPt,k) -= ((*neumann)(iPt,k)+solutionJumpPenalty)*detJ;
           }
         }
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 5d0ec00787bedf63e0fc96ab9a49305b49d407fd..fbb9a6df45ec61b96e22ecc6c87a4bbc790865c1 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -24,7 +24,7 @@ void function::call (dataCacheMap *m, fullMatrix<double> &res, std::vector<const
     default : Msg::Error("function are not implemented for %i arguments\n", arguments.size());
   }
 }
-function::function(int nbCol):_nbCol(nbCol){};
+function::function(int nbCol, bool invalidatedOnElement):_nbCol(nbCol), _invalidatedOnElement(invalidatedOnElement){};
 
 void dataCacheDouble::addMeAsDependencyOf (dataCacheDouble *newDep)
 {
@@ -37,18 +37,19 @@ void dataCacheDouble::addMeAsDependencyOf (dataCacheDouble *newDep)
   }
 }
 
-dataCacheDouble::dataCacheDouble(dataCacheMap &map,int nRowByPoint, int nCol):
-  _cacheMap(map),_value(nRowByPoint==0?1:nRowByPoint*map.getNbEvaluationPoints(),nCol)
+dataCacheDouble::dataCacheDouble(dataCacheMap &m, int nRowByPoint, int nbCol):
+  _cacheMap(m),_value(m.getNbEvaluationPoints()*nRowByPoint,nbCol)
 {
-    _nRowByPoint=nRowByPoint;
-    map.addDataCacheDouble(this);
-};
+  _nRowByPoint=nRowByPoint;
+  m.addDataCacheDouble(this,true);
+  _function = NULL;
+}
 
 dataCacheDouble::dataCacheDouble(dataCacheMap *m, function *f):
   _cacheMap(*m),_value(m->getNbEvaluationPoints(),f->getNbCol())
 {
   _nRowByPoint=1;
-  m->addDataCacheDouble(this);
+  m->addDataCacheDouble(this, f->isInvalitedOnElement());
   _function = f;
   _dependencies.resize ( _function->arguments.size());
   _depM.resize (_function->arguments.size());
@@ -62,57 +63,14 @@ void dataCacheDouble::resize() {
 
 
 //dataCacheMap members
-
-static dataCacheDouble &returnDataCacheDouble(dataCacheDouble *data, dataCacheDouble *caller)
-{
-  if(data==NULL) throw;
-  if(caller)
-    data->addMeAsDependencyOf(caller);
-  return *data;
-}
 dataCacheDouble &dataCacheMap::get(const function *f, dataCacheDouble *caller) 
 {
   dataCacheDouble *&r= _cacheDoubleMap[f];
   if(r==NULL)
     r = new dataCacheDouble(this, const_cast<function*>(f));
-  return returnDataCacheDouble(r,caller);
-}
-dataCacheDouble &dataCacheMap::getSolution(dataCacheDouble *caller) 
-{
-  return returnDataCacheDouble(_solution,caller);
-}
-dataCacheDouble &dataCacheMap::getSolutionGradient(dataCacheDouble *caller)
-{
-  return returnDataCacheDouble(_solutionGradient,caller);
-}
-dataCacheDouble &dataCacheMap::getParametricCoordinates(dataCacheDouble *caller) 
-{
-  return returnDataCacheDouble(_parametricCoordinates,caller);
-}
-dataCacheDouble &dataCacheMap::getNormals(dataCacheDouble *caller)
-{
-  returnDataCacheDouble(_normals,caller);
-  return *_normals;
-}
-
-dataCacheDouble &dataCacheMap::provideSolution(int nbFields)
-{
-  _solution = new providedDataDouble(*this,1, nbFields);
-  return *_solution;
-}
-dataCacheDouble &dataCacheMap::provideSolutionGradient(int nbFields){
-  _solutionGradient =  new providedDataDouble(*this,1, 3*nbFields);
-  return *_solutionGradient;
-}
-dataCacheDouble &dataCacheMap::provideParametricCoordinates()
-{
-  _parametricCoordinates = new providedDataDouble(*this,1, 3);
-  return *_parametricCoordinates;
-}
-dataCacheDouble &dataCacheMap::provideNormals()
-{
-  _normals = new providedDataDouble(*this,1, 3);
-  return *_normals;
+  if(caller)
+    r->addMeAsDependencyOf(caller);
+  return *r;
 }
 
 dataCacheMap::~dataCacheMap()
@@ -155,8 +113,7 @@ function *functionConstantNew(const std::vector<double> &v) {
 // get XYZ coordinates
 class functionCoordinates : public function {
   static functionCoordinates *_instance;
-  void call (dataCacheMap *m, fullMatrix<double> &xyz){
-    const fullMatrix<double> &uvw = m->getParametricCoordinates(NULL)();
+  void call (dataCacheMap *m, const fullMatrix<double> &uvw, fullMatrix<double> &xyz){
     for(int i = 0; i < uvw.size1(); i++){
       SPoint3 p;
       m->getElement()->pnt(uvw(i, 0), uvw(i, 1), uvw(i, 2), p);
@@ -165,7 +122,9 @@ class functionCoordinates : public function {
       xyz(i, 2) = p.z();
     }
   }
-  functionCoordinates():function(3){};// constructor is private only 1 instance can exists, call get to access the instance
+  functionCoordinates():function(3){
+    addArgument(function::getParametricCoordinates());
+  };// constructor is private only 1 instance can exists, call get to access the instance
  public:
   static function *get() {
     if(!_instance)
@@ -180,7 +139,8 @@ class functionSolution : public function {
   functionSolution():function(0){};// constructor is private only 1 instance can exists, call get to access the instance
  public:
   void call(dataCacheMap *m, fullMatrix<double> &sol) {
-    sol.setAsProxy(m->getSolution(NULL)());
+    Msg::Error("a function requires the solution but this algorithm does not provide the solution");
+    throw;
   }
   static function *get() {
     if(!_instance)
@@ -189,13 +149,17 @@ class functionSolution : public function {
   }
 };
 functionSolution *functionSolution::_instance = NULL;
+function *function::getSolution() {
+  return functionSolution::get();
+}
 
 class functionSolutionGradient : public function {
   static functionSolutionGradient *_instance;
   functionSolutionGradient():function(0){};// constructor is private only 1 instance can exists, call get to access the instance
  public:
   void call(dataCacheMap *m, fullMatrix<double> &sol) {
-    sol.setAsProxy(m->getSolutionGradient(NULL)());
+    Msg::Error("a function requires the gradient of the solution but this algorithm does not provide the gradient of the solution");
+    throw;
   }
   static function *get() {
     if(!_instance)
@@ -204,6 +168,47 @@ class functionSolutionGradient : public function {
   }
 };
 functionSolutionGradient *functionSolutionGradient::_instance = NULL;
+function *function::getSolutionGradient() {
+  return functionSolutionGradient::get();
+}
+
+class functionParametricCoordinates : public function {
+  static functionParametricCoordinates *_instance;
+  functionParametricCoordinates():function(0, false){};// constructor is private only 1 instance can exists, call get to access the instance
+ public:
+  void call(dataCacheMap *m, fullMatrix<double> &sol) {
+    Msg::Error("a function requires the parametric coordinates but this algorithm does not provide the parametric coordinates");
+    throw;
+  }
+  static function *get() {
+    if(!_instance)
+      _instance = new functionParametricCoordinates();
+    return _instance;
+  }
+};
+functionParametricCoordinates *functionParametricCoordinates::_instance = NULL;
+function *function::getParametricCoordinates() {
+  return functionParametricCoordinates::get();
+}
+
+class functionNormals : public function {
+  static functionNormals *_instance;
+  functionNormals():function(0){};// constructor is private only 1 instance can exists, call get to access the instance
+ public:
+  void call(dataCacheMap *m, fullMatrix<double> &sol) {
+    Msg::Error("a function requires the normals coordinates but this algorithm does not provide the normals");
+    throw;
+  }
+  static function *get() {
+    if(!_instance)
+      _instance = new functionNormals();
+    return _instance;
+  }
+};
+functionNormals *functionNormals::_instance = NULL;
+function *function::getNormals() {
+  return functionNormals::get();
+}
 
 class functionStructuredGridFile : public function {
   public:
diff --git a/Solver/function.h b/Solver/function.h
index 4cfdce9b2259e39da49c2b55c2a5523de15ceb6d..fd7b39a79c01e6485aa28b86d7bd22c77214f212 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -35,10 +35,11 @@ class dataCacheDouble;
 // more explanation at the head of this file
 class function {
   int _nbCol;
+  bool _invalidatedOnElement;
   protected :
   virtual void call (dataCacheMap *m, fullMatrix<double> &res) {throw;}
-  virtual void call (dataCacheMap *m, const fullMatrix<double> &arg0, const fullMatrix<double> &res) {throw;};
-  virtual void call (dataCacheMap *m, const fullMatrix<double> &arg0, const fullMatrix<double> &arg1, const fullMatrix<double> &res) {throw;};
+  virtual void call (dataCacheMap *m, const fullMatrix<double> &arg0, fullMatrix<double> &res) {throw;};
+  virtual void call (dataCacheMap *m, const fullMatrix<double> &arg0, const fullMatrix<double> &arg1, fullMatrix<double> &res) {throw;};
   virtual void call (dataCacheMap *m, const fullMatrix<double> &arg0, const fullMatrix<double> &arg1, const fullMatrix<double> &arg2, fullMatrix<double> &res) {throw;};
   virtual void call (dataCacheMap *m, const fullMatrix<double> &arg0, const fullMatrix<double> &arg1, const fullMatrix<double> &arg2, const fullMatrix<double> &arg3, fullMatrix<double> &res) {throw;};
   virtual void call (dataCacheMap *m, const fullMatrix<double> &arg0, const fullMatrix<double> &arg1, const fullMatrix<double> &arg2, const fullMatrix<double> &arg3, const fullMatrix<double> &arg4, fullMatrix<double> &res) {throw;};
@@ -52,8 +53,14 @@ class function {
   virtual ~function(){};
   static void registerBindings(binding *b);
   virtual void call (dataCacheMap *m, fullMatrix<double> &res, std::vector<const fullMatrix<double>*> &depM);
-  function(int nbCol);
+  function(int nbCol, bool invalidatedOnElement = true);
   inline int getNbCol()const {return _nbCol;}
+  inline bool isInvalitedOnElement() { return _invalidatedOnElement;}
+  
+  static function *getSolution();
+  static function *getSolutionGradient();
+  static function *getParametricCoordinates();
+  static function *getNormals();
 };
 
 // dataCache when the value is a  matrix of double 
@@ -65,7 +72,6 @@ class dataCacheDouble {
   std::set<dataCacheDouble*> _dependOnMe;
   std::set<dataCacheDouble*> _iDependOn;
 protected :
-  bool _valid;
   // invalidates all the cached data that depends on me
   inline void _invalidateDependencies()
   {
@@ -75,6 +81,7 @@ protected :
         (*it)->_valid=false;
   }
 public :
+  bool _valid;
   // dataCacheMap is the only one supposed to call this
   void addMeAsDependencyOf (dataCacheDouble *newDep);
   inline bool somethingDependOnMe() {
@@ -127,8 +134,8 @@ public :
     return _value;
   }
   void resize();
-  dataCacheDouble(dataCacheMap &map,int nRowByPoints, int nCol);
   dataCacheDouble(dataCacheMap *,function *f);
+  dataCacheDouble(dataCacheMap &m, int nRowByPoint, int nbCol);
   virtual ~dataCacheDouble(){};
 };
 
@@ -141,17 +148,6 @@ class dataCacheMap {
   // keep track of the current element and all the dataCaches that
   // depend on it
   std::map<const function*, dataCacheDouble*> _cacheDoubleMap;
-  class providedDataDouble : public dataCacheDouble
-  // for data provided by the algorithm and that does not have an _eval function
-  // (typically "UVW")
-  {
-    void _eval() {throw;};
-    public:
-    providedDataDouble(dataCacheMap &map, int nRowByPoints, int ncol):dataCacheDouble(map,nRowByPoints,ncol) {
-      _valid=true;
-      map._toInvalidateOnElement.erase(this);
-    }
-  };
   std::set<dataCacheDouble*> _toDelete;
   std::set<dataCacheDouble*> _toResize;
   std::set<dataCacheDouble*> _toInvalidateOnElement;
@@ -162,21 +158,12 @@ class dataCacheMap {
   void addDataCache(dataCacheDouble *data){
     _toDelete.insert(data);
   }
-  void addDataCacheDouble(dataCacheDouble *data){
+  void addDataCacheDouble(dataCacheDouble *data, bool invalidatedOnElement){
     _toResize.insert(data);
-    _toInvalidateOnElement.insert(data);
+    if(invalidatedOnElement)
+      _toInvalidateOnElement.insert(data);
   }
  public:
-  dataCacheDouble *_solution, *_solutionGradient, *_parametricCoordinates, *_normals;
-  dataCacheDouble &getSolution(dataCacheDouble *caller);
-  dataCacheDouble &getSolutionGradient(dataCacheDouble *caller);
-  dataCacheDouble &getParametricCoordinates(dataCacheDouble *caller);
-  dataCacheDouble &getNormals(dataCacheDouble *caller);
-  dataCacheDouble &provideSolution(int nbFields);
-  dataCacheDouble &provideSolutionGradient(int nbFields);
-  dataCacheDouble &provideParametricCoordinates();
-  dataCacheDouble &provideNormals();
-
   dataCacheDouble &get(const function *f, dataCacheDouble *caller=0);
   inline void setElement(MElement *element) {
     _element=element;
@@ -186,7 +173,6 @@ class dataCacheMap {
   }
   inline MElement *getElement() {return _element;}
   dataCacheMap(){
-    _normals = _solution = _solutionGradient = _parametricCoordinates = 0;
     _nbEvaluationPoints = 0;
   }
   void setNbEvaluationPoints(int nbEvaluationPoints);