diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt index c34eb4aef75aaf305c2100b2201e827f58c7d8b0..694438eda9a5a984a2d31052d0c93a38f64d8db4 100644 --- a/Solver/CMakeLists.txt +++ b/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) diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp index 6e119af54709816b96470916f81004cd8ba8241a..ffef08e52305f9c26abbd03127684a4379ccfe60 100644 --- a/Solver/dgAlgorithm.cpp +++ b/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; } } diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h index fd1441f14e9b1dab4a6f5a3a893803d88e91b382..481d125f23aaed14e338942f17df178b2f642864 100644 --- a/Solver/dgConservationLaw.h +++ b/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 { diff --git a/Solver/dgMainTest.cpp b/Solver/dgMainTest.cpp index cd53ba9c7ac9c1caf47a3d25a6102de78d95c421..3d75e92b0e62eb7bf32d1ed670e5fe02e3ef3667 100644 --- a/Solver/dgMainTest.cpp +++ b/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 diff --git a/Solver/function.cpp b/Solver/function.cpp index 02803ae152a1e8672000c42de0f6b41e556a262f..977ed647982820a215a9f79ab8c2a34fd9c409f4 100644 --- a/Solver/function.cpp +++ b/Solver/function.cpp @@ -1,23 +1,99 @@ #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); - } + diff --git a/Solver/function.h b/Solver/function.h index ab7b6ab56fabdd7113e2252230429fb8cb328820..bc756ea772cd913e3d8c7ae5869d3965451335be 100644 --- a/Solver/function.h +++ b/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