dg : functor to replace function (simpler and OpenMP compatible)

git-svn-id: https://geuz.org/svn/gmsh/trunk@19168 59c8cabf-10d6-4a28-b43f-3936b85956c4

CMakeLists.txt

@@ -161,7 +161,6 @@ add_subdirectory(modules/shallowWater)
 add_subdirectory(modules/slim3d)
 add_subdirectory(modules/slimFunction)
 add_subdirectory(modules/slimIce)
-add_subdirectory(modules/subsurface)
 add_subdirectory(modules/supra)
 add_subdirectory(modules/wave)
 

benchmarks/Supra/DiffusionSupra3Dvectorial.py

@@ -204,7 +204,7 @@ solve.getNewtonZ().setVerb(10)
 
 
 startcpu = time.clock()
-for i in range (0, 10001):
+for i in range (0, 101):
  
   #print'*** solve ***'
   #sys = linearSystemPETScBlockDouble ()

conservationLaw/dgConservationLaw.h

@@ -14,7 +14,7 @@ template<class t> class dgFullMatrix;
 #include <vector>
 
 //@TODO : function and dgExtrusion should not be needed
-class function;
+class functor;
 class dgExtrusion;
 
 class dgConservationLaw {
@@ -59,8 +59,8 @@ class dgConservationLaw {
   /** @todo everything bellow still relies on functions (but is optional) => should be adapted:
    *for dgLimiter
    */
-  virtual const function *getOutsideValue(const dgGroupOfFaces*, const function *in) const {return NULL;}
-  virtual const function *getClipToPhysics(const std::string tag) const {return NULL;}
+  virtual const functor *getOutsideValue(const dgGroupOfFaces*, const functor *in) const {return NULL;}
+  virtual const functor *getClipToPhysics(const std::string tag) const {return NULL;}
 };
 
 #endif

conservationLaw/dgConservationLawFunction.h

@@ -3,30 +3,63 @@
 
 #include "dgConservationLaw.h"
 #include "GmshMessage.h"
+#include "fullMatrix.h"
 #include <map>
 
+class functor;
 class dataCacheMap;
-class function;
-class dgDofContainer;
-class dgDofManager;
+class dgConservationLawFunction;
 class dgTermAssembler;
 class dgFaceTermAssembler;
-class dgExtrusion;
-class dgBoundaryCondition;
-class dgStrongBoundaryCondition;
-class dgGroupCollection;
-template <class t> class fullMatrix;
+/*==============================================================================
+ * Class for boundary conditions
+ *============================================================================*/
+
+/**
+ * A weak boundary condition. 
+ * Boundary conditions should be associated with tag using dgConservationLawFunction::addBoundaryCondition.
+ */
+class dgBoundaryCondition {
+ protected:
+  dgConservationLawFunction *_claw;
+  const functor *_term0, *_term1N;
+  std::string _type;
+  std::map<const functor*,const functor*> _outsideValue;
+ public: 
+  
+  dgBoundaryCondition () {
+    _term0 = NULL;
+    _term1N = NULL;
+    _type = "undefined";
+  }
+  virtual ~dgBoundaryCondition () {}
+  
+  // Get
+  virtual const functor * getTerm0 () { return _term0; }
+  virtual const functor * getTerm1N () { return _term1N; }
+  virtual const functor * getOutsideValue (const functor* f) { return _outsideValue[f]; } 
+  inline std::string getType() const { return _type; }
+};
+
 
 /**
- * A conservation law is defined a convective flux (f), a diffusive flux(g), a source term(r)
- * and a set of boundary conditions.\n
- * \f[\partial_t L(u) = \nabla \cdot (\mathbf{f}(u,\nabla u,forcings)) + r(u,forcings))\f]
+ * A strong boundary condition. 
  */
+class dgStrongBoundaryCondition {
+ public:
+  int dim;
+  std::string physicalTag;
+  const functor *f;
+  dgStrongBoundaryCondition (int dim_, const std::string physicalTag_, const functor *f_) :
+      dim(dim_), physicalTag(physicalTag_), f(f_) {};
+};
+
+
 class dgConservationLawFunction : public dgConservationLaw {
   friend class dgConservationLawFunctionGather;
   
  public:
-  typedef std::map<const std::string, const function*> termMap;
+  typedef std::map<const std::string, const functor*> termMap;
   
  private:
   std::vector<dgStrongBoundaryCondition> _strongBoundaryConditions;
@@ -54,16 +87,16 @@ class dgConservationLawFunction : public dgConservationLaw {
   
   std::vector<std::string> _undeletableBoundaryTerm0;
   
-  std::map<const std::string, std::pair<const function*, massFactorType> > _massFactor;
+  std::map<const std::string, std::pair<const functor*, massFactorType> > _massFactor;
   
-  static inline const function * functionForTag(const termMap &map, const std::string tag)
+  static inline const functor * functionForTag(const termMap &map, const std::string tag)
   {
     termMap::const_iterator it = map.find(tag);
     if (it == map.end()) it = map.find("");
     if (it == map.end()) return NULL;
     return it->second;
   }
-  static inline const function * functionForTagEnforced(const termMap &map, const std::string tag)
+  static inline const functor * functionForTagEnforced(const termMap &map, const std::string tag)
   {
     termMap::const_iterator it = map.find(tag);
     if (it == map.end()) return NULL;
@@ -98,16 +131,16 @@ class dgConservationLawFunction : public dgConservationLaw {
   
   virtual dgTermAssembler *getTerm(const dgGroupOfElements &group) const;
   virtual dgFaceTermAssembler *getFaceTerm(const dgGroupOfFaces &group) const;
-  const function *getClipToPhysics(const std::string tag) const;
-  const function *getInterfaceTerm0(const std::string tag) const;
-  const function *getdiffusivity(const std::string tag) const;
-  const function *getMaximumConvectiveSpeed(const std::string tag) const;
+  const functor *getClipToPhysics(const std::string tag) const;
+  const functor *getInterfaceTerm0(const std::string tag) const;
+  const functor *getdiffusivity(const std::string tag) const;
+  const functor *getMaximumConvectiveSpeed(const std::string tag) const;
   
   /** Associate a mesh tag with a boundary condition.*/
   void addBoundaryCondition(std::string tag, dgBoundaryCondition *condition);
   void addBoundaryCondition(const std::vector<std::string> tags, dgBoundaryCondition *condition);
   dgBoundaryCondition *getBoundaryCondition(const std::string tag, bool allowNull=false) const;
-  virtual const function *getOutsideValue(const dgGroupOfFaces*, const function *in) const;
+  virtual const functor *getOutsideValue(const dgGroupOfFaces*, const functor *in) const;
   bool haveFaceTerm(const dgGroupOfFaces &faces) const;
   
   inline void checkBoundaryConditionUndefined(std::string TAG) const {
@@ -119,17 +152,17 @@ class dgConservationLawFunction : public dgConservationLaw {
    * Create a new boundary condition based on external data.
    * The interface term with the provided tag computes the fluxes using the provided values as external solution");
    */
-  dgBoundaryCondition *newOutsideValueBoundary(const std::string tag, const function *outsideValueFunction);
-  dgBoundaryCondition *newOutsideValueBoundary(const function *outsideValueFunction)
+  dgBoundaryCondition *newOutsideValueBoundary(const std::string tag, const functor *outsideValueFunction);
+  dgBoundaryCondition *newOutsideValueBoundary(const functor *outsideValueFunction)
     {return newOutsideValueBoundary("", outsideValueFunction);}
   
   /**
    * Create a new boundary condition based on external data. Generic version
    * The interface term with the provided tag computes the fluxes using the provided values as external solution");
    */
-  dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::string tag, const std::vector<const function*> toReplace,
-                                                                             const std::vector<const function*> replaceBy);
-  dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::vector<const function*> toReplace, const std::vector<const function*> replaceBy)
+  dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::string tag, const std::vector<const functor*> toReplace,
+                                                                             const std::vector<const functor*> replaceBy);
+  dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::vector<const functor*> toReplace, const std::vector<const functor*> replaceBy)
     {return newOutsideValueBoundaryGeneric("", toReplace, replaceBy);}
   
   /**
@@ -137,18 +170,18 @@ class dgConservationLawFunction : public dgConservationLaw {
    * The interface term with the provided tag computes the fluxes using the provided values as external solution");
    */
   dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::string tag,
-                                                      const std::vector<const function*> toReplace1,
-                                                      const std::vector<const function*> replaceBy1,
-                                                      const std::vector<const function*> toReplace2,
-                                                      const std::vector<const function*> replaceBy2);
+                                                      const std::vector<const functor*> toReplace1,
+                                                      const std::vector<const functor*> replaceBy1,
+                                                      const std::vector<const functor*> toReplace2,
+                                                      const std::vector<const functor*> replaceBy2);
   dgBoundaryCondition *newOutsideValueBoundaryGeneric(const std::string tag,
-                                                      const std::map<const function*, const function*> *map,
-                                                      const std::map<const function*, const function*> *map2=NULL);
+                                                      const std::map<const functor*, const functor*> *map,
+                                                      const std::map<const functor*, const functor*> *map2=NULL);
   /**
    * Create a new boundary condition with a given flux (no other fluxes will be computed)
    */
-  dgBoundaryCondition *newNeumannBoundary(const function *fluxFunction);
-  void setBoundaryNeumann(std::string, const function *);
+  dgBoundaryCondition *newNeumannBoundary(const functor *fluxFunction);
+  void setBoundaryNeumann(std::string, const functor *);
   
   /**
    * Create a new boundary condition which set to 0 all the fluxes through the interfaces.
@@ -168,7 +201,7 @@ class dgConservationLawFunction : public dgConservationLaw {
    * Fields with flags >0 are fixed
    * Strong boundary condition are not implemented in Discontinuous Galerkin
    */
-  void addStrongBoundaryCondition (int dim, const std::string tag, const function *f);
+  void addStrongBoundaryCondition (int dim, const std::string tag, const functor *f);
   
   /** Return a fullMatrix list of the group's elements' innerRadii. */
   static fullMatrix<double> innerRadius(dataCacheMap *m);
@@ -196,10 +229,10 @@ class dgConservationLawFunction : public dgConservationLaw {
       of a 2d mesh and horizontal/vertical phenomena have different timescales. */
   double getMinOfTimeSteps(dgDofContainer *solution, dgExtrusion *extrusion = NULL) const;
   void fixStrongBC(dgDofManager &dof, dgDofContainer *oldSolution = NULL, double alpha = 1.) const;
-  const function * getMassFactor(const std::string tag, massFactorType &type) const;
+  const functor * getMassFactor(const std::string tag, massFactorType &type) const;
   massFactorType haveMassFactor(const dgGroupOfElements &group) const;
-  void setVolumeTerm(int dorder, const function *f, const std::string tag="");
-  void setInterfaceTerm(int dorder, const function *f, const std::string tag="");
+  void setVolumeTerm(int dorder, const functor *f, const std::string tag="");
+  void setInterfaceTerm(int dorder, const functor *f, const std::string tag="");
 };
 
 class dgConservationLawPython : public dgConservationLawFunction {
@@ -212,59 +245,15 @@ class dgConservationLawPython : public dgConservationLawFunction {
 };
 
 
-/*==============================================================================
- * Class for boundary conditions
- *============================================================================*/
-
-/**
- * A weak boundary condition. 
- * Boundary conditions should be associated with tag using dgConservationLawFunction::addBoundaryCondition.
- */
-class dgBoundaryCondition {
- protected:
-  dgConservationLawFunction *_claw;
-  const function *_term0, *_term1N;
-  std::string _type;
-  std::map<const function*,const function*> _outsideValue;
- public: 
-  
-  dgBoundaryCondition () {
-    _term0 = NULL;
-    _term1N = NULL;
-    _type = "undefined";
-  }
-  virtual ~dgBoundaryCondition () {}
-  
-  // Get
-  virtual const function * getTerm0 () { return _term0; }
-  virtual const function * getTerm1N () { return _term1N; }
-  virtual const function * getOutsideValue (const function* f) { return _outsideValue[f]; } 
-  inline std::string getType() const { return _type; }
-};
-
-
-/**
- * A strong boundary condition. 
- */
-class dgStrongBoundaryCondition {
- public:
-  int dim;
-  std::string physicalTag;
-  const function *f;
-  dgStrongBoundaryCondition (int dim_, const std::string physicalTag_, const function *f_) :
-      dim(dim_), physicalTag(physicalTag_), f(f_) {};
-};
-
-
 /*==============================================================================
  * IP Term
  *============================================================================*/
 
-function *dgNewIpTerm (int nbFields, const function *diffusiveFlux, const function *diffusivity, bool _useSip=true);
-function *dgNewIpTermAnisotropic(int nbFields, const function *diffusiveFlux, const function *ipPenalty);
-function *dgNewIpTermErnRef(int nbFields, const function *diffusiveFlux,const function *ipPenalty);
-function *dgNewSymmetricIpTerm (int nbFields, const function *diffusivity);
-function *dgNewSymIpTerm3d(int nbFields, const function *diffusivity);
+functor *dgNewIpTerm (int nbFields, const functor *diffusiveFlux, const functor *diffusivity, bool _useSip=true);
+functor *dgNewIpTermAnisotropic(int nbFields, const functor *diffusiveFlux, const functor *ipPenalty);
+functor *dgNewIpTermErnRef(int nbFields, const functor *diffusiveFlux,const functor *ipPenalty);
+functor *dgNewSymmetricIpTerm (int nbFields, const functor *diffusivity);
+functor *dgNewSymIpTerm3d(int nbFields, const functor *diffusivity);
 
 #endif
 

conservationLaw/dgConservationLawFunctionGather.cpp

@@ -22,7 +22,7 @@ class dgConservationLawFunctionGather::term: public function {
   std::vector<fullMatrix<double> > _childTerms;
   int _dim;
   public:
-  term (lawList &laws, std::vector<const function*> &childs, int size, int dim):function(size*dim),_dim(dim){
+  term (lawList &laws, std::vector<const functor*> &childs, int size, int dim):function(size*dim),_dim(dim){
     fflush(stdout);
     _childTerms.resize(childs.size());
     for (size_t iLaw=0; iLaw < childs.size(); iLaw++){

conservationLaw/dgConservationLawFunctionGather.h

 #ifndef _DG_CONSERVATION_LAW_FUNCTION_GATHER_H_
 #define _DG_CONSERVATION_LAW_FUNCTION_GATHER_H_
 #include "dgConservationLawFunction.h"
-#include "function.h"
 #include "dgDofContainer.h"
 
 /** Gathers several conservation laws together in a single one */
 class dgConservationLawFunctionGather : public dgConservationLawFunction {
  public:
   typedef std::vector<const dgConservationLawFunction*> lawList;
-  typedef std::map<std::string, std::vector<const function*> > functionMat;
+  typedef std::map<std::string, std::vector<const functor*> > functionMat;
  protected:
   std::vector<const dgConservationLawFunction*> iFieldToLaw;
   std::vector<int> gToLField;

conservationLaw/dgCoupling.cpp

@@ -79,7 +79,7 @@ function *dgCoupling_SW1D_SW2D::computeOutsideSolution2New(const function *integ
   return new computeOutsideSolution2(integrals);
 }
 
-void dgDofOverAnotherCollectionFunction::call(dataCacheMap *m, dataCacheDouble *d, fullMatrix<double> &val) {
+void dgDofOverAnotherCollectionFunction::call(dataCacheMap *m, fullMatrix<double> &val) {
   int ig, ie;
   int nbFields = _dofC->nField(*m->getGroupOfElements());
 
@@ -156,7 +156,7 @@ void dgDofOverAnotherCollectionFunction::call(dataCacheMap *m, dataCacheDouble *
   }
 }
 
-dgCoupling::dgCoupling(dgFunctionCoupled * f, dataCacheMap::mode m) : _f(f) {
+dgCoupling::dgCoupling(dgFunctionCoupled * f, dataCacheMap::Mode m) : _f(f) {
   GModel * model = NULL;
   // check if the same GModel is used
   for(int i = 0; i < _f->nFunction(); ++i) {

conservationLaw/dgCoupling.h

@@ -134,7 +134,7 @@ class dgDofOverAnotherCollectionFunction : public function {
   const dgDofContainer * _dofC;
  public:
   dgDofOverAnotherCollectionFunction(const dgDofContainer * dofC) : function(dofC->nFieldMax()), _dofC(dofC) {}
-  void call(dataCacheMap *m, dataCacheDouble *d, fullMatrix<double> &val);
+  void call(dataCacheMap *m, fullMatrix<double> &val);
 };
 
 class dgFunctionCoupled {
@@ -183,7 +183,7 @@ class dgCoupling {
   std::map<const dgGroupCollection *, function *> _vFunc;
   int _dimFace;
  public:
-  dgCoupling(dgFunctionCoupled * f, dataCacheMap::mode m = dataCacheMap::INTEGRATION_GROUP_MODE);
+  dgCoupling(dgFunctionCoupled * f, dataCacheMap::Mode m = dataCacheMap::INTEGRATION_GROUP_MODE);
   void compute(std::vector<std::string> physicalTags=std::vector<std::string>());
 
   /* get dgFunctionCoupled in the required group collection */

conservationLaw/dgTerm.cpp

@@ -41,9 +41,9 @@ dgInterfaceTerm::~dgInterfaceTerm() {}
 
 class dgPsiTerm : public dgTerm {
   bool _isLinear;
-  const function *_f;
+  const functor *_f;
 public:
-  dgPsiTerm(const function *f, bool isLinear) {
+  dgPsiTerm(const functor *f, bool isLinear) {
     _f = f;
     _isLinear = isLinear;
   }
@@ -52,12 +52,13 @@ public:
     dataCacheDouble &term = *cache.get(_f);
     const dgIntegrationMatrices &integ = cache.getIntegrationMatrices();
     const dgMeshJacobian &meshj = cache.getJacobians();
+    //term.get().print("element psi");
     residual.gemm(integ.psiW(), meshj.multJElement(cache.getGroupId(), term.get()));
-    if (jac && term.doIDependOn(function::getSolution())) {
+    if (jac && cache.dependOn(*_f, *function::getSolution())) {
       functionDerivator derivator(term, *cache.get(function::getSolution()), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE);
       gemmAndAllocate(*jac, integ.psiPsiW(), meshj.multJElement(cache.getGroupId(), derivator.compute()));
     }
-    if (jac && term.doIDependOn(function::getSolutionGradient())) {
+    if (jac && cache.dependOn(*_f, *function::getSolutionGradient())) {
       functionDerivator derivatorGradSol(term, *cache.get(function::getSolutionGradient()), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE);
       gemmAndAllocate(*jac, integ.dPsiPsiW(), meshj.multDXiDXJRightElement(cache.getGroupId(), derivatorGradSol.compute()));
     }
@@ -66,9 +67,9 @@ public:
 
 class dgGradPsiTerm : public dgTerm {
   bool _isLinear;
-  const function *_f;
+  const functor *_f;
   public:
-  dgGradPsiTerm(const function *f, bool isLinear) {
+  dgGradPsiTerm(const functor *f, bool isLinear) {
     _f = f;
     _isLinear = isLinear;
   }
@@ -77,12 +78,13 @@ class dgGradPsiTerm : public dgTerm {
     dataCacheDouble &term = *cache.get(_f);
     const dgIntegrationMatrices &integ = cache.getIntegrationMatrices();
     const dgMeshJacobian &meshj = cache.getJacobians();
+    //term.get().print("element grad psi");
     residual.gemm(integ.dPsiW(), meshj.multDXiDXJElement(cache.getGroupId(), term.get()));
-    if (jac && term.doIDependOn(function::getSolution())) {
+    if (jac && cache.dependOn(*_f, *function::getSolution())) {
       functionDerivator derivator(term, *cache.get(function::getSolution()), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE);
       gemmAndAllocate(*jac, integ.psiDPsiW(),meshj.multDXiDXJElement(cache.getGroupId(), derivator.compute()));
     }
-    if (jac && term.doIDependOn(function::getSolutionGradient())) {
+    if (jac && cache.dependOn(*_f, *function::getSolutionGradient())) {
       functionDerivator derivatorGradSol(term, *cache.get(function::getSolutionGradient()), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE);
       gemmAndAllocate(*jac, integ.dPsiDPsiW(), meshj.multDXiDXDXiDXJElement(cache.getGroupId(), derivatorGradSol.compute()));
     }
@@ -90,10 +92,10 @@ class dgGradPsiTerm : public dgTerm {
 };
 
 class dgInterfacePsiTerm:public dgInterfaceTerm {
-  const function *_f;
+  const functor *_f;
   bool _isLinear;
   public:
-  dgInterfacePsiTerm(const function *f, bool isLinear) {
+  dgInterfacePsiTerm(const functor *f, bool isLinear) {
     _isLinear = isLinear;
     _f = f;
   }
@@ -102,24 +104,26 @@ class dgInterfacePsiTerm:public dgInterfaceTerm {
     dataCacheDouble &term = *cache.get(_f);
     const dgGroupOfFaces &group = *cache.getGroupOfInterfaces();
     const dgMeshJacobian &meshj = cache.getJacobians();
+    //printf("%i %s\n", group.nConnection(), group.physicalTag().c_str());
+    //term.get().print();
     gemmAndAllocate(residual, cache.getIntegrationMatrices().psiW(), meshj.multJInterface(group.interfaceVectorId(), group.interfaceId(), term.get()));
     if (jacSol) {
       jacSol->resize(group.nConnection());
       for (int i = 0; i < group.nConnection(); i++) {
-        if (term.doIDependOn(function::getSolution(), i)) {
-          functionDerivator derivator(term, *cache.getSecondaryCache(i)->get(function::getSolution()), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE);
-          const fullMatrix<double> &dSolAtQP = meshj.multJInterface(group.interfaceVectorId(),group.interfaceId(), derivator.compute());
-          gemmAndAllocate((*jacSol)[i], cache.getIntegrationMatrices().psiPsiW(), dSolAtQP);
+        if (cache.dependOn(*_f, *function::getSolution(), -1, i)) {
+          const fullMatrix<double> &d = cache.computeDerivative(*_f, *function::getSolution(), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE, i);
+          const fullMatrix<double> &dSolAtQP = meshj.multJInterface(group.interfaceVectorId(),group.interfaceId(), d);
+          gemmAndAllocate((*jacSol)[i], cache.integrationMatrices().psiPsiW(), dSolAtQP);
         }
       }
     }
     if (jacGradSol) {
       jacGradSol->resize(group.nConnection());
       for (int i = 0; i < group.nConnection(); i++) {
-        if (term.doIDependOn(function::getSolutionGradient(), i)) {
-          functionDerivator derivator(term, *cache.getSecondaryCache(i)->get(function::getSolutionGradient()), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE);
-          const fullMatrix<double> &dGradSolAtQP = meshj.multDXiDXJRightInterface(group.interfaceVectorId(), i, group.interfaceId(), derivator.compute());
-          gemmAndAllocate((*jacGradSol)[i], cache.getSecondaryCache(i)->getIntegrationMatrices().dPsiPsiW(), dGradSolAtQP);
+        if (cache.dependOn(*_f, *function::getSolutionGradient(), -1, i)) {
+          const fullMatrix<double> &d = cache.computeDerivative(*_f, *function::getSolutionGradient(), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE, i);
+          const fullMatrix<double> &dGradSolAtQP = meshj.multDXiDXJRightInterface(group.interfaceVectorId(), i, group.interfaceId(), d);
+          gemmAndAllocate((*jacGradSol)[i], cache.integrationMatrices(i).dPsiPsiW(), dGradSolAtQP);
         }
       }
     }
@@ -172,7 +176,7 @@ class dgLowOrderFPsiTerm : public dgTerm {
   protected:
     bool _onlyVertical;
     int _order;
-    const function *_f;
+    const functor *_f;
     int _nbFields;
     bool _isLinear;
 
@@ -296,7 +300,7 @@ class dgLowOrderFPsiTerm : public dgTerm {
       residual.gemm(cache.getIntegrationMatrices().psiW(), cache.getJacobians().multJElement(cache.getGroupId(), termAtQP));
     }
 
-    if (jac &&term.doIDependOn(function::getSolution())) {
+    if (jac && cache.dependOn(*_f, *function::getSolution())) {
       fullMatrix<double> highIntegToLowIntegPsi = buildLowOrderJacobianMatrix(cache.getIntegrationMatrices(), highIntegToLowInteg);
       functionDerivator derivator(term, *cache.get(function::getSolution()), _isLinear ? 1 : EPSILON_FINITE_DIFFERENCE);
       fullMatrix<double> termAtQP;
@@ -309,12 +313,12 @@ class dgLowOrderFPsiTerm : public dgTerm {
       tmp.reshape(nPsi * nPsi, -1);
       jac->add(tmp);
     }
-    if (term.doIDependOn(function::getSolutionGradient())) {
+    if (cache.dependOn(*_f, *function::getSolutionGradient())) {
       Msg::Fatal("integralLowOrderFTerm should not be used with terms depending on solution gradient");
     }
   }
   public:
-    dgLowOrderFPsiTerm(const function *f, int nbFields, bool isLinear, int order, bool onlyVertical)
+    dgLowOrderFPsiTerm(const functor *f, int nbFields, bool isLinear, int order, bool onlyVertical)
   {
     _order = order;
     _nbFields = nbFields;
@@ -324,32 +328,32 @@ class dgLowOrderFPsiTerm : public dgTerm {
   }
 };
 
-dgTerm *dgIntegralTerm::newLowOrderFPsiTerm(const function *f, int nbFields, bool isLinear, int order, bool onlyVertical){
+dgTerm *dgIntegralTerm::newLowOrderFPsiTerm(const functor *f, int nbFields, bool isLinear, int order, bool onlyVertical){
   return new dgLowOrderFPsiTerm(f, nbFields, isLinear, order, onlyVertical);
 }
 
 
-dgTerm *dgIntegralTerm::newPsiTerm(const function *f, int nbFields, bool isLinear)
+dgTerm *dgIntegralTerm::newPsiTerm(const functor *f, int nbFields, bool isLinear)
 {
   return new dgPsiTerm(f, isLinear);
 }
 
-dgInterfaceTerm *dgIntegralTerm::newInterfacePsiTerm(const function *f, bool isLinear)
+dgInterfaceTerm *dgIntegralTerm::newInterfacePsiTerm(const functor *f, bool isLinear)
 {
   return new dgInterfacePsiTerm(f, isLinear);
 }
 
-dgTerm *dgIntegralTerm::newGradPsiTerm(const function *f, int nbFields, bool isLinear)
+dgTerm *dgIntegralTerm::newGradPsiTerm(const functor *f, int nbFields, bool isLinear)
 {
   return new dgGradPsiTerm(f, isLinear);
 }
 
-dgTerm *dgIntegralTerm::newPsiTermNodal(const function *f, int nbFields, bool isLinear)
+dgTerm *dgIntegralTerm::newPsiTermNodal(const functor *f, int nbFields, bool isLinear)
 {
   return NULL;//new dgPsiTermNodal(f, nbFields, isLinear);
 }
 
-dgTerm *dgIntegralTerm::newGradPsiTermNodal(const function *f, int nbFields, bool isLinear)
+dgTerm *dgIntegralTerm::newGradPsiTermNodal(const functor *f, int nbFields, bool isLinear)
 {
   return NULL;//new dgGradPsiTermNodal(f, nbFields, isLinear);
 }
@@ -424,7 +428,7 @@ void dgFaceTermAssemblerIntegral::compute(const dgGroupOfFaces &faces, double t,
   if (! _term)
     return;
   int nbConnections = faces.nConnection();
-  dataCacheMap cacheMap(dataCacheMap::INTEGRATION_GROUP_MODE, residual.getGroups(), nbConnections - 1);
+  dataCacheMap cacheMap(dataCacheMap::INTEGRATION_GROUP_MODE, residual.getGroups());
   function::getTime()->set(t);
   cacheMap.setInterfaceGroup(&faces);
   cacheMap.setSolutionFunction(solution.getFunction(), solution.getFunctionGradient());

conservationLaw/dgTerm.h

@@ -21,13 +21,13 @@ class dgInterfaceTerm {
 
 class dgIntegralTerm {
  public :
-  static dgTerm *newPsiTerm(const function *f, int nbFields, bool isLinear);
-  static dgTerm *newGradPsiTerm(const function *f, int nbFields, bool isLinear);
+  static dgTerm *newPsiTerm(const functor *f, int nbFields, bool isLinear);
+  static dgTerm *newGradPsiTerm(const functor *f, int nbFields, bool isLinear);
   // linearize f before integration (same as interpolate)
-  static dgTerm *newPsiTermNodal(const function *f, int nbFields, bool isLinear);
-  static dgTerm *newGradPsiTermNodal(const function *f, int nbFields, bool isLinear);
-  static dgTerm *newLowOrderFPsiTerm(const function *f, int nbFields, bool isLinear, int order, bool onlyVertical);
-  static dgInterfaceTerm *newInterfacePsiTerm(const function *f, bool isLinear);
+  static dgTerm *newPsiTermNodal(const functor *f, int nbFields, bool isLinear);
+  static dgTerm *newGradPsiTermNodal(const functor *f, int nbFields, bool isLinear);
+  static dgTerm *newLowOrderFPsiTerm(const functor *f, int nbFields, bool isLinear, int order, bool onlyVertical);
+  static dgInterfaceTerm *newInterfacePsiTerm(const functor *f, bool isLinear);
 };
 
 

dofContainer/dgDofContainer.cpp

@@ -232,7 +232,7 @@ static fullMatrix<double> transposeElementField(const dgFullMatrix<double> &in)
 }
 
 
-void dgDofFunction::call(dataCacheMap *m, dataCacheDouble *d, fullMatrix<double> &sol)
+void dgDofFunction::call(dataCacheMap *m, fullMatrix<double> &sol)
 {
   switch(m->getMode()) {
     case dataCacheMap::INTEGRATION_GROUP_MODE:
@@ -241,6 +241,10 @@ void dgDofFunction::call(dataCacheMap *m, dataCacheDouble *d, fullMatrix<double>
         _dof->_getAtIntegrationPoints(*m->getGroupOfElements(), m->getIntegrationMatrices(), sol);
       }
       else {
+        if (m->connectionId() >= m->getGroupOfInterfaces()->nConnection()) {
+          m->invalidCall();