dg: add the number of evaluating points on dataCacheMap

Solver/dgAlgorithm.cpp

@@ -40,7 +40,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
   fullMatrix<double> Fuvw[3] = {fullMatrix<double> ( group.getNbIntegrationPoints(), group.getNbElements() * nbFields),
 				fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields),
 				fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields)};
-  dataCacheMap cacheMap;
+  dataCacheMap cacheMap(group.getNbIntegrationPoints());
   dataCacheElement &cacheElement = cacheMap.getElement();
   // provided dataCache
   cacheMap.provideData("UVW").set(group.getIntegrationPointsMatrix());
@@ -146,7 +146,8 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
   fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), nbFaces*nbFields*2);
 
   // create one dataCache for each side
-  dataCacheMap cacheMapLeft, cacheMapRight;
+  dataCacheMap cacheMapLeft(group.getNbIntegrationPoints());
+  dataCacheMap cacheMapRight(group.getNbIntegrationPoints());
 
   // data this algorithm provide to the cache map (we can maybe move each data to a separate function but
   // I think It's easier like this)
@@ -310,7 +311,7 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
   // ----- 2 ----  compute normal fluxes  at integration points
   fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields);
 
-  dataCacheMap cacheMapLeft;
+  dataCacheMap cacheMapLeft(group.getNbIntegrationPoints());
   // provided dataCache
   cacheMapLeft.provideData("UVW").set(group.getIntegrationPointsMatrix());
   dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");

Solver/dgConservationLaw.cpp

@@ -4,10 +4,10 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
   dgConservationLaw &_claw;
   std::string _outsideValueFunctionName;
   class term : public dataCacheDouble {
-    dataCacheMap cacheMapRight; // new cacheMap to  pass to the Riemann solver
     dataCacheDouble &solutionRight;
     dataCacheDouble &solutionLeft;
     dataCacheDouble &outsideValue;
+    dataCacheMap cacheMapRight; // new cacheMap to  pass to the Riemann solver
     dataCacheDouble *riemannSolver;
     dgConservationLaw &_claw;
     public:
@@ -15,6 +15,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
       solutionRight(cacheMapRight.provideData("Solution")),
       solutionLeft(cacheMapLeft.get("Solution",this)),
       outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
+      cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
       _claw(claw)
     {
       riemannSolver=_claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
@@ -26,10 +27,12 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
         _value = fullMatrix<double>(solutionLeft().size1(),_claw.nbFields());
       }
       solutionRight.set(outsideValue());
+      if(riemannSolver){
         for(int i=0;i<_value.size1(); i++)
           for(int j=0;j<_value.size2(); j++)
             _value(i,j) = (*riemannSolver)(i,j*2);
       }
+    }
   };
   public:
   dgBoundaryConditionOutsideValue(dgConservationLaw &claw,const std::string outsideValueFunctionName): _claw(claw),

Solver/function.cpp

@@ -86,10 +86,9 @@ class functionXYZ : public function {
     dataCacheDouble &_uvw;
    public:
     data(dataCacheMap *m) : 
+      dataCacheDouble(m->getNbEvaluationPoints(),3),
       _element(m->getElement(this)), _uvw(m->get("UVW", this))
-    {
-      _value = fullMatrix<double> (_uvw().size1(), 3);
-    }
+    {}
     void _eval()
     {
       for(int i = 0; i < _uvw().size1(); i++){
@@ -113,15 +112,12 @@ class functionConstant : public function {
  private :
   class data : public dataCacheDouble {
     const functionConstant *_function;
-    dataCacheDouble &_uvw;
     public:
-    data(const functionConstant * function,dataCacheMap *m) :_uvw(m->get("UVW",this)){
+    data(const functionConstant * function,dataCacheMap *m):
+      dataCacheDouble(m->getNbEvaluationPoints(),function->_source.size1()){
       _function = function;
     }
     void _eval() {
-      if(_value.size1()!=_uvw().size1()){
-        _value=fullMatrix<double>(_uvw().size1(),_function->_source.size1());
-      }
       for(int i=0;i<_value.size1();i++)
         for(int j=0;j<_function->_source.size1();j++)
           _value(i,j)=_function->_source(j,0);

Solver/function.h

@@ -97,6 +97,8 @@ class dataCacheDouble : public dataCache {
     }
     return _value;
   }
+  dataCacheDouble(){};
+  dataCacheDouble(int size1, int size2):_value(size1,size2){};
   virtual ~dataCacheDouble(){};
 };
 
@@ -133,6 +135,7 @@ class dataCacheElement : public dataCache {
 // more explanation at the head of this file
 class dataCacheMap {
  private:
+  int _nbEvaluationPoints;
   // keep track of the current element and all the dataCaches that
   // depend on it
   dataCacheElement _cacheElement;
@@ -140,7 +143,7 @@ class dataCacheMap {
   class providedDataDouble : public dataCacheDouble
   // for data provided by the algorithm and that does not have an _eval function
   // (typically "UVW") this class is not stricly necessary, we could write
-  // a function for each case but I think it's more practical like this
+  // a function for each case 
   {
     void _eval() {throw;};
     public:
@@ -152,6 +155,8 @@ class dataCacheMap {
   dataCacheDouble &get(const std::string &functionName, dataCache *caller=0);
   dataCacheElement &getElement(dataCache *caller=0);
   dataCacheDouble &provideData(std::string name);
+  dataCacheMap(int nbEvaluationPoints):_nbEvaluationPoints(nbEvaluationPoints){}
+  inline int getNbEvaluationPoints(){return _nbEvaluationPoints;}
   ~dataCacheMap();
 };
 #endif

Solver/luaFunction.cpp

@@ -13,12 +13,12 @@ class functionLua : public function {
  private:
   class data : public dataCacheDouble{
     private:    
-    dataCacheDouble &_uvw;
     const functionLua *_function;
     std::vector<dataCacheDouble *> _dependencies;
     public:
-    data(const functionLua *f, dataCacheMap *m)
-      : _function(f),_uvw(m->get("UVW",this))
+    data(const functionLua *f, dataCacheMap *m):
+      dataCacheDouble(m->getNbEvaluationPoints(),f->_nbCol),
+      _function(f)    
     {
       _dependencies.resize ( _function->_dependenciesName.size());
       for (int i=0;i<_function->_dependenciesName.size();i++)
@@ -26,9 +26,6 @@ class functionLua : public function {
     }
     void _eval()
     {
-      if(_value.size1()!=_uvw().size1() || _value.size2() != _function->_nbCol){
-        _value=fullMatrix<double>(_uvw().size1(),_function->_nbCol);
-      }
       lua_getfield(_function->_L, LUA_GLOBALSINDEX, _function->_luaFunctionName.c_str());
       for (int i=0;i< _dependencies.size();i++){
         const fullMatrix<double> *data = &(*_dependencies[i])();