implementation krabbenhoft(2007) for Richards equation

Solver/function.cpp

@@ -206,7 +206,7 @@ void dataCacheDouble::_eval()
 const function * dataCacheMap::_translate(const function *f) const 
 {
   //special case
-  if (f == function::getSolution()) {
+  if (f == function::getSolution() || f == _containerSolution) {
     f = _functionSolution;
     if (f == NULL) {
       dataCacheMap *parent = _parent;
@@ -218,7 +218,7 @@ const function * dataCacheMap::_translate(const function *f) const
       if (f == NULL) 
         Msg::Error ("solution function has not been set");
     }
-  } else if (f == function::getSolutionGradient()) {
+  } else if (f == function::getSolutionGradient() || f == _containerSolutionGradient) {
     f = _functionSolutionGradient;
     if (f == NULL) {
       dataCacheMap *parent = _parent;

Solver/function.h

@@ -201,7 +201,7 @@ class dataCacheDouble {
 };
 
 class dataCacheMap {
-  const function *_functionSolution, *_functionSolutionGradient, *_functionCoordinates;
+  const function *_functionSolution, *_functionSolutionGradient, *_functionCoordinates, *_containerSolution, *_containerSolutionGradient;
   //handle function solution and funciton solution gradient
   //we should get rid of them
   const function * _translate (const function *) const;
@@ -215,9 +215,9 @@ class dataCacheMap {
   std::vector<dataCacheDouble*> _toInvalidateOnElement;
   MElement *_element;
   dataCacheMap() {
-    _functionSolution = _functionSolutionGradient = _functionCoordinates = NULL;
+    _functionSolution = _functionSolutionGradient = _functionCoordinates = _containerSolution = _containerSolutionGradient = NULL;
     _nbEvaluationPoints = 0;
-    _parent=NULL;
+    _parent = NULL;
   }
   ~dataCacheMap();
   void addDataCacheDouble(dataCacheDouble *data, bool invalidatedOnElement)
@@ -269,6 +269,13 @@ class dataCacheMap {
     _functionSolution = functionSolution;
     _functionSolutionGradient = functionSolutionGradient;
   }
+  inline void setReferenceSolutionFunction(const function *functionSolution, const function *functionSolutionGradient) {
+    _containerSolution = functionSolution;
+    _containerSolutionGradient = functionSolutionGradient;
+    for(std::list<dataCacheMap*>::iterator it = _children.begin(); it != _children.end(); it++) {
+      (*it)->setReferenceSolutionFunction(functionSolution, functionSolutionGradient);
+    }
+  }
   void setNbEvaluationPoints(int nbEvaluationPoints);
   inline int getNbEvaluationPoints() { return _nbEvaluationPoints; }
 };