reduce interface of dataCacheDouble

Solver/dgAlgorithm.cpp

@@ -84,7 +84,7 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
   cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
   dataCacheDouble &sol = cacheMap.provideSolution(nbFields);
   dataCacheDouble &UVW = cacheMap.provideParametricCoordinates();
-  UVW.set(group.getIntegrationPointsMatrix());
+  UVW.set()=group.getIntegrationPointsMatrix();
   // provided dataCache
   dataCacheDouble *maxConvectiveSpeed = claw.newMaxConvectiveSpeed(cacheMap);
   dataCacheDouble *maximumDiffusivity = claw.newMaximumDiffusivity(cacheMap);
@@ -98,7 +98,7 @@ 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.setAsProxy(solution, iElement*nbFields, nbFields);
+    sol.set().setAsProxy(solution, iElement*nbFields, nbFields);
     cacheMap.setElement(group.getElement(iElement));
     const double L = group.getInnerRadius(iElement);
     double spectralRadius = 0.0;

Solver/dgConservationLaw.cpp

@@ -22,7 +22,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
     }
 
     void _eval() {
-      solutionRight.set(outsideValue());
+      solutionRight.set()=outsideValue();
       if(riemannSolver){
         for(int i=0;i<_value.size1(); i++)
           for(int j=0;j<_value.size2(); j++)
@@ -62,7 +62,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
     }
     void _eval() {
       if(maxDif){
-        solutionRight.set(outsideValue());
+        solutionRight.set()=outsideValue();
         for(int i=0;i<_value.size1(); i++)
           for(int j=0;j<_value.size2(); j++)
             _value(i,j) = (*maxDif)(i,j);

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.provideParametricCoordinates().set()=group.getIntegrationPointsMatrix();
     dataCacheDouble &sourceTerm = cacheMap.get(f);
     fullMatrix<double> source;
     for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {

Solver/dgFunctionIntegrator.cpp

@@ -20,13 +20,13 @@ void dgFunctionIntegrator::compute(dgDofContainer *sol,fullMatrix<double> &resul
     dgGroupOfElements &group=*sol->getGroups()->getElementGroup(iGroup);
     fullMatrix<double> &solProxy=sol->getGroupProxy(&group);
     cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
-    UVW.set(group.getIntegrationPointsMatrix());
+    UVW.set()=group.getIntegrationPointsMatrix();
     fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
     group.getCollocationMatrix().mult(solProxy  , solutionQP); 
     fullMatrix<double> IPMatrix = group.getIntegrationPointsMatrix();
     for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
       cacheMap.setElement(group.getElement(iElement));
-      solutionQPe.setAsProxy(solutionQP, iElement*nbFields, nbFields );
+      solutionQPe.set().setAsProxy(solutionQP, iElement*nbFields, nbFields );
       for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
         const double detJ = group.getDetJ (iElement, iPt);
         for (int k=0;k<nbRowResult;k++){

Solver/dgLimiter.cpp

@@ -118,8 +118,8 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution)
     dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
     if (solutionEClipped){
       for (int iElement=0 ; iElement<egroup->getNbElements() ;++iElement) {
-        solutionE.setAsProxy(solGroup, iElement*nbFields, nbFields );
-        solutionE.set((*solutionEClipped)());    
+        solutionE.set().setAsProxy(solGroup, iElement*nbFields, nbFields );
+        solutionE.set()=(*solutionEClipped)();
       }
     }
   }  

Solver/dgResidual.cpp

@@ -28,7 +28,7 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
 {
   residual.scale(0);
   _cacheMap->setNbEvaluationPoints(group.getNbIntegrationPoints());
-  _UVW.set(group.getIntegrationPointsMatrix());
+  _UVW.set()=group.getIntegrationPointsMatrix();
   // ----- 1 ----  get the solution at quadrature points
   // ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints) 
   fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
@@ -53,7 +53,7 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
   // ----- 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
-    _solutionQPe.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
+    _solutionQPe.set().setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
 
     if(_gradientSolutionQPe.somethingDependOnMe()){
       dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
@@ -112,7 +112,7 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
 {
   residual.scale(0);
   _cacheMap->setNbEvaluationPoints(group.getNbIntegrationPoints());
-  _UVW.set(group.getIntegrationPointsMatrix());
+  _UVW.set()=group.getIntegrationPointsMatrix();
   // ----- 1 ----  get the solution at quadrature points
   // ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints) 
   fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
@@ -145,7 +145,7 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
   // ----- 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
-    _solutionQPe.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
+    _solutionQPe.set().setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
     a.setAsProxy(A, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
     b.setAsProxy(B, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
 
@@ -321,9 +321,9 @@ 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).setAsProxy(connections[i]->getIntegrationPointsOnElement(iFace));
-      caches[i].getSolution(NULL).setAsProxy(solutionQP, (iFace*nbConnections+i)*_nbFields, _nbFields);
-      caches[i].getNormals(NULL).setAsProxy(connections[i]->getNormals(), iFace*group.getNbIntegrationPoints(), group.getNbIntegrationPoints());
+      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());
       // B2 ) compute the gradient of the solution
       if(caches[i].getSolutionGradient(NULL).somethingDependOnMe()) {
         dofs.setAsProxy(*solutionOnElements[i], _nbFields*connections[i]->getElementId(iFace), _nbFields);

Solver/function.h

@@ -76,13 +76,6 @@ public :
   inline bool somethingDependOnMe() {
     return !_dependOnMe.empty();
   }
-  inline int howManyDependOnMe() {
-    return _dependOnMe.size();
-  }
-  inline int howManyDoIDependOn() {
-    return _iDependOn.size();
-  }
-
   std::vector<dataCacheDouble*> _dependencies;
   std::vector<const fullMatrix<double>*> _depM;
 
@@ -92,6 +85,7 @@ public :
  protected:
   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(int i=0;i<_dependencies.size(); i++)
@@ -100,31 +94,11 @@ public :
   }
  public:
   //set the value (without computing it by _eval) and invalidate the dependencies
-  inline void set(const fullMatrix<double> &mat) {
-    if(_valid)
-      _invalidateDependencies();
-    _value=mat;
-    _valid=true;
-  }
-  // take care if you use this you must ensure that the value pointed to are not modified
-  // without further call to setAsProxy because the dependencies won't be invalidate
-  inline void setAsProxy(const fullMatrix<double> &mat, int cShift, int c) {
-    if(_valid)
-      _invalidateDependencies();
-    _value.setAsProxy(mat,cShift,c);
-    _valid=true;
-  }
-  // take care if you use this you must ensure that the value pointed to are not modified
-  // without further call to setAsProxy because the dependencies won't be invalidate
-  inline void setAsProxy(const fullMatrix<double> &mat) {
-    if(_valid)
-      _invalidateDependencies();
-    _value.setAsProxy(mat,0,mat.size2());
-    _valid=true;
-  }
   // this function is needed to be able to pass the _value to functions like gemm or mult
   // but you cannot keep the reference to the _value, you should always use the set function 
   // to modify the _value
+  // take care if you use this to set a proxy you must ensure that the value pointed to are not modified
+  // without further call to set because the dependencies won't be invalidate
   inline fullMatrix<double> &set() {
     if(_valid)
       _invalidateDependencies();

Solver/functionDerivator.cpp

@@ -9,10 +9,10 @@ void functionDerivator::compute() {
     _xDx = _xRef;
     for (int i=0;i<_fRef.size1();i++)
       _xDx(i,j) += _epsilon;
-    _x.set(_xDx);
+    _x.set()=_xDx;
     for (int i=0; i<_fRef.size1(); i++) 
       for (int k=0; k<_fRef.size2(); k++) 
         _dfdx(i, k*_xRef.size2()+j) = (_f(i,k)-_fRef(i,k))/_epsilon;
   }
-  _x.set(_xRef);
+  _x.set()=_xRef;
 };