diff --git a/Solver/function.cpp b/Solver/function.cpp
index 9b0ca8199ddb0eb6b472053e8e9272080fea7c7f..c135be15a815f291234402d5cb276da703cac1e5 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -10,16 +10,9 @@
 #endif
 #include "Bindings.h"
 
-void function::call (dataCacheMap *m, fullMatrix<double> &res, std::vector<const fullMatrix<double>*> &depM) {
-  switch (arguments.size()) {
-    case 0 : call(m, res); break;
-    case 1 : call(m, *depM[0], res); break;
-    case 2 : call(m, *depM[0], *depM[1], res); break;
-    case 3 : call(m, *depM[0], *depM[1], *depM[2], res); break;
-    case 4 : call(m, *depM[0], *depM[1], *depM[2], *depM[3], res); break;
-    case 5 : call(m, *depM[0], *depM[1], *depM[2], *depM[3], *depM[4], res); break;
-    case 6 : call(m, *depM[0], *depM[1], *depM[2], *depM[3], *depM[4], *depM[5], res); break;
-    default : Msg::Error("function are not implemented for %i arguments\n", arguments.size());
+function::~function() {
+  for (int i=0; i<arguments.size(); i++) {
+    delete arguments[i];
   }
 }
 function::function(int nbCol, bool invalidatedOnElement):_nbCol(nbCol), _invalidatedOnElement(invalidatedOnElement){};
@@ -50,15 +43,23 @@ dataCacheDouble::dataCacheDouble(dataCacheMap *m, function *f):
   m->addDataCacheDouble(this, f->isInvalitedOnElement());
   _function = f;
   _dependencies.resize ( _function->arguments.size());
-  _depM.resize (_function->arguments.size());
-  for (unsigned int i=0;i<_function->arguments.size();i++)
-    _dependencies[i] = &m[_function->arguments[i].first].get(_function->arguments[i].second,this);
+  for (unsigned int i=0;i<_function->arguments.size();i++) {
+    int iCache = _function->arguments[i]->iMap;
+    const function *f = _function->arguments[i]->f;
+    _dependencies[i] = &m->getSecondaryCache(iCache)->get(f,this);
+  }
 }
 
 void dataCacheDouble::resize() {
   _value = fullMatrix<double>(_nRowByPoint==0?1:_nRowByPoint*_cacheMap.getNbEvaluationPoints(),_value.size2());
 }
 
+void dataCacheDouble::_eval() {
+  for(unsigned int i=0;i<_dependencies.size(); i++){
+    _function->arguments[i]->val.setAsProxy((*_dependencies[i])());
+  }
+  _function->call(&_cacheMap, _value);
+}
 
 //dataCacheMap members
 dataCacheDouble &dataCacheMap::get(const function *f, dataCacheDouble *caller) 
@@ -106,8 +107,9 @@ class functionConstant : public function {
   fullMatrix<double> _source;
   void call(dataCacheMap *m, fullMatrix<double> &val) {
     for(int i=0;i<val.size1();i++)
-      for(int j=0;j<_source.size1();j++)
+      for(int j=0;j<_source.size1();j++){
         val(i,j)=_source(j,0);
+        }
   }
   functionConstant(std::vector<double> source):function(source.size()){
     _source = fullMatrix<double>(source.size(),1);
@@ -130,7 +132,8 @@ function *functionConstantNew(const std::vector<double> &v) {
 // get XYZ coordinates
 class functionCoordinates : public function {
   static functionCoordinates *_instance;
-  void call (dataCacheMap *m, const fullMatrix<double> &uvw, fullMatrix<double> &xyz){
+  const fullMatrix<double> &uvw;
+  void call (dataCacheMap *m, 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);
@@ -139,8 +142,8 @@ class functionCoordinates : public function {
       xyz(i, 2) = p.z();
     }
   }
-  functionCoordinates():function(3){
-    addArgument(function::getParametricCoordinates());
+  functionCoordinates():function(3),
+    uvw(addArgument(function::getParametricCoordinates())){
   };// constructor is private only 1 instance can exists, call get to access the instance
  public:
   static function *get() {
@@ -228,6 +231,7 @@ function *function::getNormals() {
 }
 
 class functionStructuredGridFile : public function {
+  const fullMatrix<double> &coord;
   public:
   int n[3];
   double d[3],o[3];
@@ -235,7 +239,7 @@ class functionStructuredGridFile : public function {
     return v[(i*n[1]+j)*n[2]+k];
   }
   double *v;
-  void call(dataCacheMap *m, const fullMatrix<double> &coord, fullMatrix<double> &val){
+  void call(dataCacheMap *m, fullMatrix<double> &val){
     for(int pt=0;pt<val.size1();pt++){
       double xi[3];
       int id[3];
@@ -257,9 +261,8 @@ class functionStructuredGridFile : public function {
         +get(id[0]+1 ,id[1]+1 ,id[2]+1 )*(  xi[0])*(  xi[1])*(  xi[2]);
     }
   }
-  functionStructuredGridFile(const std::string filename, const function *coordFunction): function(1){
+  functionStructuredGridFile(const std::string filename, const function *coordFunction): function(1), coord(addArgument(coordFunction)){
     std::ifstream input(filename.c_str());
-    addArgument(coordFunction);
     if(!input)
       Msg::Error("cannot open file : %s",filename.c_str());
     if(filename.substr(filename.size()-4,4)!=".bin") {
@@ -288,12 +291,12 @@ class functionLua : public function {
   lua_State *_L;
   std::string _luaFunctionName;
  public:
-  void call (dataCacheMap *m, fullMatrix<double> &res, std::vector<const fullMatrix<double>*> &depM) {
+  void call (dataCacheMap *m, fullMatrix<double> &res) {
     lua_getfield(_L, LUA_GLOBALSINDEX, _luaFunctionName.c_str());
-    for (int i=0;i< depM.size(); i++)
-      luaStack<const fullMatrix<double>*>::push(_L, depM[i]);
+    for (int i=0;i< arguments.size(); i++)
+      luaStack<const fullMatrix<double>*>::push(_L, &arguments[i]->val);
     luaStack<const fullMatrix<double>*>::push(_L, &res);
-    lua_call(_L, depM.size()+1, 0);
+    lua_call(_L, arguments.size()+1, 0);
   }
   functionLua (int nbCol, std::string luaFunctionName, std::vector<const function*> dependencies, lua_State *L)
     : function(nbCol), _luaFunctionName(luaFunctionName), _L(L)
@@ -322,40 +325,40 @@ void dataCacheMap::setNbEvaluationPoints(int nbEvaluationPoints) {
 class functionC : public function {
   void (*callback)(void);
   public:
-  void call (dataCacheMap *m, fullMatrix<double> &val, std::vector<const fullMatrix<double>*> &depM) {
-    switch (depM.size()) {
+  void call (dataCacheMap *m, fullMatrix<double> &val) {
+    switch (arguments.size()) {
       case 0 : 
         ((void (*)(fullMatrix<double> &))(callback))(val);
         break;
       case 1 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&))
-         (callback)) (val, *depM[0]);
+         (callback)) (val, arguments[0]->val);
         break;
       case 2 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double> &))
-         (callback)) (val, *depM[0], *depM[1]);
+         (callback)) (val, arguments[0]->val, arguments[1]->val);
         break;
       case 3 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&))
-         (callback)) (val, *depM[0], *depM[1], *depM[2]);
+         (callback)) (val, arguments[0]->val, arguments[1]->val, arguments[2]->val);
         break;
       case 4 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&,
                    const fullMatrix<double>&))
-         (callback)) (val, *depM[0], *depM[1], *depM[2], *depM[3]);
+         (callback)) (val, arguments[0]->val, arguments[1]->val, arguments[2]->val, arguments[3]->val);
         break;
       case 5 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&,
                    const fullMatrix<double>&, const fullMatrix<double>&))
-         (callback)) (val, *depM[0], *depM[1], *depM[2], *depM[3], *depM[4]);
+         (callback)) (val, arguments[0]->val, arguments[1]->val, arguments[2]->val, arguments[3]->val, arguments[4]->val);
         break;
       case 6 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&,
                    const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&))
-         (callback)) (val, *depM[0], *depM[1], *depM[2], *depM[3], *depM[4], *depM[5]);
+         (callback)) (val, arguments[0]->val, arguments[1]->val, arguments[2]->val, arguments[3]->val, arguments[4]->val, arguments[5]->val);
         break;
       default :
-        Msg::Error("C callback not implemented for %i argurments", depM.size());
+        Msg::Error("C callback not implemented for %i argurments", arguments.size());
     }
   }
   functionC (std::string file, std::string symbol, int nbCol, std::vector<const function *> dependencies):
diff --git a/Solver/function.h b/Solver/function.h
index 77e32267240e7f5653a670528d2ccb93dcb36ab7..193eddeee47097e6a022a99a7128c45b3e5f774d 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -34,25 +34,29 @@ class dataCacheDouble;
 // An abstract interface to functions 
 // more explanation at the head of this file
 class function {
+  class argument {
+    //iMap is the id of the dataCacheMap, e.g. on interfaces
+    public:
+    int iMap;
+    const function *f;
+    fullMatrix<double> val;
+    argument(int iMap_, const function *f_) {
+      iMap = iMap_;
+      f = f_;
+    }
+  };
   int _nbCol;
   bool _invalidatedOnElement;
   protected :
-  virtual void call (dataCacheMap *m, 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;};
-  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, const fullMatrix<double> &arg5, fullMatrix<double> &res) {throw;};
   public :
-  std::vector<std::pair<int, const function*> > arguments;
-  void addArgument(const function *f, int iMap = 0) {
-    //iMap is the id of the dataCacheMap, e.g. on interfaces
-    arguments.push_back(std::pair<int, const function*>(iMap, f));
+  virtual void call (dataCacheMap *m, fullMatrix<double> &res)=0;
+  std::vector<argument*> arguments;
+  const fullMatrix<double> &addArgument(const function *f, int iMap = 0) {
+    arguments.push_back(new argument(iMap, f));
+    return arguments.back()->val;
   }
-  virtual ~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, bool invalidatedOnElement = true);
   inline int getNbCol()const {return _nbCol;}
   inline bool isInvalitedOnElement() { return _invalidatedOnElement;}
@@ -91,21 +95,15 @@ public :
     return (_iDependOn.find(&other)!=_iDependOn.end());
   }
   std::vector<dataCacheDouble*> _dependencies;
-  std::vector<const fullMatrix<double>*> _depM;
 
   int _nRowByPoint;
-  dataCacheMap &_cacheMap;
   function *_function;
  protected:
+  dataCacheMap &_cacheMap;
   fullMatrix<double> _value;
   // do the actual computation and put the result into _value
   // still virtual because it is overrided by conservation law terms, as soon as conservation law terms will be regular functions, we will remove this
-  virtual void _eval()
-  {
-    for(unsigned int i=0;i<_dependencies.size(); i++)
-      _depM[i] = &(*_dependencies[i])();
-    _function->call(&_cacheMap, _value, _depM);
-  }
+  virtual void _eval();
  public:
   //set the value (without computing it by _eval) and invalidate the dependencies
   // this function is needed to be able to pass the _value to functions like gemm or mult
@@ -143,11 +141,13 @@ public :
 };
 
 
+class dgDataCacheMap;
 // more explanation at the head of this file
 class dataCacheMap {
   friend class dataCacheDouble;
   dataCacheMap  *_parent;
   std::list<dataCacheMap*> _children;
+  std::vector<dataCacheMap*> _secondaryCaches;
   int _nbEvaluationPoints;
   std::map<const function*, dataCacheDouble*> _cacheDoubleMap;
   std::set<dataCacheDouble*> _allDataCaches;
@@ -162,9 +162,21 @@ class dataCacheMap {
       _toInvalidateOnElement.insert(data);
   }
  public:
+  virtual dgDataCacheMap *asDgDataCacheMap() {
+    Msg::Error("I'm not a dgDataCacheMap\n");
+    return NULL;
+  }
+  dataCacheMap *getSecondaryCache(int i) {
+    if (i==0)
+      return this;
+    return _secondaryCaches[i-1];
+  }
+  void addSecondaryCache(dataCacheMap *s) {
+    _secondaryCaches.push_back(s);
+  }
   dataCacheDouble &get(const function *f, dataCacheDouble *caller=0);
   dataCacheDouble &substitute(const function *f);
-  inline void setElement(MElement *element) {
+  virtual void setElement(MElement *element) {
     _element=element;
     for(std::set<dataCacheDouble*>::iterator it = _toInvalidateOnElement.begin(); it!= _toInvalidateOnElement.end(); it++) {
       (*it)->_valid=false;