gmshdg : add functions and replace the sourceTerm by a normal function

other terms will follows....

Solver/CMakeLists.txt

@@ -14,6 +14,7 @@ set(SRC
   dgGroupOfElements.cpp
   dgAlgorithm.cpp
   dgConservationLawAdvection.cpp
+  function.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) 

Solver/dgAlgorithm.cpp

@@ -6,6 +6,7 @@
 #include "MElement.h"
 #include "GModel.h"
 #include "MEdge.h"
+#include "function.h"
 /*
   compute 
     \int \vec{f} \cdot \grad \phi dv   
@@ -44,6 +45,11 @@ 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;
+  dataCacheElement &cacheElement = cacheMap.getElement();
+  cacheMap.provideData("UVW").set(group.getIntegrationPointsMatrix());
+  dataCacheDouble *sourceTerm = claw.newSourceTerm(cacheMap);
+
   // ----- 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
@@ -52,6 +58,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
     if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradientSolutionQP, 3*iElement*nbFields,3*nbFields );      
     else gradSolutionQPe = new fullMatrix<double>;
     dgElement DGE( group.getElement(iElement), solutionQPe, *gradSolutionQPe, group.getIntegrationPointsMatrix());
+    cacheElement.set(group.getElement(iElement));
     if(claw.convectiveFlux() || claw.diffusiveFlux()) {
       // ----- 2.3.2 --- compute fluxes in XYZ coordinates
       if (claw.convectiveFlux()) (*claw.convectiveFlux())(DGE,fConv);
@@ -75,13 +82,13 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
         } 
       }
     }
-    if (claw.sourceTerm()){
+    if (sourceTerm){
       fullMatrix<double> source(Source, iElement*claw.nbFields(),claw.nbFields());
-      (*claw.sourceTerm())(DGE,&source); 
       for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
         const double detJ = group.getDetJ (iElement, iPt);
-        for (int k=0;k<nbFields;k++)
-          source(iPt,k) *= detJ;
+        for (int k=0;k<nbFields;k++){
+          source(iPt,k) = (*sourceTerm)(iPt,k)*detJ;
+        }
       }
     }
     delete gradSolutionQPe;
@@ -92,8 +99,9 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
       residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
     }
   }
-  if(claw.sourceTerm()){
+  if(sourceTerm){
     residual.gemm(group.getSourceRedistributionMatrix(),Source);
+    delete sourceTerm;
   }
 }
 

Solver/dgConservationLaw.h

@@ -9,6 +9,7 @@
 #include "fullMatrix.h"
 class dgElement; 
 class dgFace;
+class dataCacheDouble;
 
 class dgTerm{
  public:
@@ -30,6 +31,7 @@ public:
   virtual boundaryType  type() const =0;
 };
 
+class dataCacheMap;
 class dgConservationLaw {
   protected :
   int _nbf;
@@ -43,9 +45,10 @@ public:
 
   int nbFields() const {return _nbf;}
 
+  virtual dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const { return NULL; }
+
   inline const dgTerm     * convectiveFlux () const {return _convective;}
   inline const dgTerm     * diffusiveFlux  () const {return _diffusive;}
-  inline const dgTerm     * sourceTerm     () const {return _source;}
   inline const dgFaceTerm * riemannSolver  () const {return _riemannSolver;}
   inline const dgTerm     * maxConvectiveSpeed () const {return _maxConvectiveSpeed;}
   inline const dgBoundaryCondition *boundaryCondition(const std::string tag) const {

Solver/dgMainTest.cpp

@@ -5,31 +5,33 @@
 #include "dgAlgorithm.h"
 #include "dgConservationLaw.h"
 #include "Gmsh.h"
+#include "function.h"
 
 
 #include "MElement.h"
 void print (const char *filename,const dgGroupOfElements &els, double *v);
 
-class testSourceTerm : public dgTerm {
-  void operator () (const dgElement &el, fullMatrix<double> fcx[]) const{
-    const fullMatrix<double> &sol = el.solution();
-    const fullMatrix<double> &qp = el.integration();
-    SPoint3 p;
-    for(int i=0; i< sol.size1(); i++) {
-      el.element()->pnt(qp(i,0),qp(i,1),qp(i,2),p);
-      fcx[0](i,0)=exp(-(pow(p.x()-0.2,2)+pow(p.y()-0.3,2))*100);
-    }
-  }
-};
 
 class dgConservationLawInitialCondition : public dgConservationLaw {
+  class gaussian : public dataCacheDouble {
+    dataCacheDouble &xyz;
+    double _xc,_yc;
+    public:
+    gaussian(dataCacheMap &cacheMap,double xc, double yc):xyz(cacheMap.get("XYZ",this)),_xc(xc),_yc(yc){};
+    void _eval () { 
+      if(_value.size1() != xyz().size1())
+        _value=fullMatrix<double>(xyz().size1(),1);
+      for(int i=0; i< _value.size1(); i++) {
+        _value(i,0)=exp(-(pow(xyz(i,0)-_xc,2)+pow(xyz(i,1)-_yc,2))*100);
+      }
+    }
+  };
   public:
   dgConservationLawInitialCondition() {
     _nbf = 1;
-    _source = new testSourceTerm;
   }
-  ~dgConservationLawInitialCondition() {
-    delete _source;
+  dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
+    return new gaussian(cacheMap,0.2,0.3);
   }
 };
 
@@ -42,6 +44,7 @@ int main(int argc, char **argv){
   int order=1;
   int dimension=2;
   dgAlgorithm algo;
+  function::registerAllFunctions();
   algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
 
   //for now, we suppose there is only one group of elements

Solver/function.cpp

 #include "function.h"
+#include "SPoint3.h"
+#include "MElement.h"
 
-class functionXYZ : public function{
+// dataCache members
+
+void dataCache::addMeAsDependencyOf (dataCache *newDep)
+{
+  _dependOnMe.insert(&newDep->_valid);
+  newDep->_iDependOn.insert(this);
+  for(std::set<dataCache*>::iterator it = _iDependOn.begin();
+      it != _iDependOn.end(); it++) {
+    (*it)->_dependOnMe.insert(&newDep->_valid);
+    newDep->_iDependOn.insert(*it);
+  }
+}
+
+// function static members
+
+std::map<std::string,function*> function::_allFunctions;
+
+bool function::add(const std::string functionName, function *f)
+{
+  if(_allFunctions.find(functionName)!=_allFunctions.end())
+    return false;
+  _allFunctions[functionName]=f;
+  return true;
+}
+
+function *function::get(std::string functionName, bool acceptNull)
+{
+  std::map<std::string, function*>::iterator it=_allFunctions.find(functionName);
+  if (it==_allFunctions.end()) {
+    if (acceptNull)
+      return NULL;
+    else
+      throw;
+  }
+  return it->second;
+}
+
+//dataCacheMap members
+
+dataCacheElement &dataCacheMap::getElement(dataCache *caller) 
+{
+  if(caller)
+    _cacheElement.addMeAsDependencyOf(caller);
+  return _cacheElement;
+}
+
+dataCacheDouble &dataCacheMap::get(const std::string &functionName, dataCache *caller) 
+{
+  dataCacheDouble *&r= _cacheDoubleMap[functionName];
+  if(r==NULL)
+    r = function::get(functionName)->newDataCache(this);
+  if(caller)
+    r->addMeAsDependencyOf(caller);
+  return *r;
+}
+
+dataCacheDouble &dataCacheMap::provideData(std::string name)
+{
+  dataCacheDouble *&r= _cacheDoubleMap[name];
+  if(r!=NULL)
+    throw;
+  r = new providedDataDouble;
+  return *r;
+}
+
+dataCacheMap::~dataCacheMap()
+{
+  for (std::map<std::string, dataCacheDouble*>::iterator it = _cacheDoubleMap.begin();
+      it!=_cacheDoubleMap.end(); it++) {
+    delete it->second;
+  }
+}
+
+// now some example of functions
+
+class functionXYZ : public function {
  private:
   class data : public dataCacheDouble{
    private:
     dataCacheElement &_element;
     dataCacheDouble &_uvw;
    public:
-    data(dataCacheMap *m) 
-      : dataCache(m), _element(m->getElement(this)), _uvw(m->get("UVW", this))
+    data(dataCacheMap *m) :
+      _element(m->getElement(this)), _uvw(m->get("UVW", this))
     {
-      _value.resize(_uvw.size1(), 3);
+      _value = fullMatrix<double> (_uvw().size1(), 3);
     }
-    void eval()
+    void _eval()
     {
-      MElement *ele = _element();
-      for(int i = 0; i < _uvw.size1(); i++){
+      for(int i = 0; i < _uvw().size1(); i++){
         SPoint3 p;
-        ele->pnt(_uvw(i, 0), _uvw(i, 1), _uvw(i, 2), p);
+        _element()->pnt(_uvw(i, 0), _uvw(i, 1), _uvw(i, 2), p);
         _value(i, 0) = p.x();
         _value(i, 1) = p.y();
         _value(i, 2) = p.z();
@@ -25,15 +101,15 @@ class functionXYZ : public function{
     }
   };
  public:
-  dataCache &getDataCache(dataCacheMap *m)
+  dataCacheDouble *newDataCache(dataCacheMap *m)
   {
     return new data(m);
   }
 };
 
 
-void registerFunctions()
+void function::registerAllFunctions()
 {
   function::add("XYZ", new functionXYZ);
-
 }
+

Solver/function.h

@@ -2,43 +2,79 @@
 #define _FUNCTION_H_
 
 #include <map>
+#include <set>
 #include <string>
+#include <fullMatrix.h>
+class dataCacheMap;
+class MElement;
 
-// An abstract interface to functions
-class function{
- private:
-  static std::map<std::string, function*> _allFunctions;
-  std::string _name;
- public:
-  virtual dataCache &getDataCache(dataCacheMap *m) = 0;
-  static bool add(std::string functionName, function *f);
-};
-
-// A class to store function data and cache it intelligently (by
-// computing dependencies)
-class dataCache{
- protected:
+// A class to store function data and cache it (by computing dependencies)
+class dataCache {
+  friend class dataCacheMap;
+  // pointers to the "_valid" flag of all dataCache depending on me
+  std::set<bool*> _dependOnMe;
+  std::set<dataCache*> _iDependOn;
+protected :
   bool _valid;
-  std::set<dataCache*> _iDependOnThem, _theyDependOnMe;
-  // validates all the dependencies
-  void _validate();
-  void _invalidateDependencies();
+  // invalidates all the cached data that depends on me
+  inline void _invalidateDependencies()
+  {
+    // if this is too slow we can keep a C array cache of the _dependOnMe set
+    for(std::set<bool*>::iterator it = _dependOnMe.begin();
+      it!=_dependOnMe.end(); it++)
+        **it=false;
+  }
+  // dataCacheMap is the only one supposed to call this
+  void addMeAsDependencyOf (dataCache *newDep);
+  dataCache() : _valid(false) {}
+  virtual ~dataCache(){};
 };
 
 // A node in the dependency tree for which all the leafs depend on the
 // given double value
 class dataCacheDouble : public dataCache {
- private:
+ protected:
   fullMatrix<double> _value;
+  // do the actual computation and put the result into _value
+  virtual void _eval()=0;
  public:
-  void set(const fullMatrix<double> &mat);
-  void set(int i, int j, double val);
+  //set the value (without computing it by _eval) and invalidate the dependencies
+  inline void set(const fullMatrix<double> &mat) {
+    _invalidateDependencies();
+    _value=mat;
+    _valid=true;
+  }
+  inline const double &operator () (int i, int j)
+  {
+    if(!_valid) {
+      _eval();
+      _valid=true;
+    }
+    return _value(i,j);
+  }
+  //access _value and compute it if necessary
   inline const fullMatrix<double> &operator () ()
   {
-    if(!_valid) _validate();
+    if(!_valid) {
+      _eval();
+      _valid=true;
+    }
     return _value;
   }
-  virtual void eval() = 0;
+  virtual ~dataCacheDouble(){};
+};
+
+// An abstract interface to functions
+class function {
+ private:
+  static std::map<std::string, function*> _allFunctions;
+  std::string _name;
+ public:
+  static void registerAllFunctions();
+  static bool add(const std::string functionName, function *f);
+  static function *get(std::string functionName, bool acceptNull=false);
+
+  virtual dataCacheDouble *newDataCache(dataCacheMap *m) =0;
 };
 
 // A special node in the dependency tree for which all the leafs
@@ -47,20 +83,30 @@ class dataCacheElement : public dataCache {
  private:
   MElement *_element;
  public:
-  void set(const MElement *ele);
-  inline const MElement *operator () () { return _element; }
+  void set(MElement *ele){
+    _invalidateDependencies();
+    _element=ele;
+  };
+  inline MElement *operator () () { return _element; }
 };
 
-class dataCacheMap{
+class dataCacheMap {
  private:
   // keep track of the current element and all the dataCaches that
   // depend on it
-  dataCacheElement *_dependOnElement;
+  dataCacheElement _cacheElement;
+  std::map<std::string, dataCacheDouble*> _cacheDoubleMap;
+  class providedDataDouble : public dataCacheDouble
+  {
+    void _eval() {throw;};
+  };
  public:
-  dataCacheDouble &get(const std::string &functionName,
-                       dataCache *caller=0);
+  dataCacheDouble &get(const std::string &functionName, dataCache *caller=0);
   dataCacheElement &getElement(dataCache *caller=0);
-
+  // for data provided by the user 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
+  dataCacheDouble &provideData(std::string name);
+  ~dataCacheMap();
 };
-
 #endif