diff --git a/Solver/TESTCASES/Advection.lua b/Solver/TESTCASES/Advection.lua
index d66da7451f3bf76a2f5302b844f596646b47c26e..5d4f9f02047361edb2b3df09038e1a968f7aebbe 100644
--- a/Solver/TESTCASES/Advection.lua
+++ b/Solver/TESTCASES/Advection.lua
@@ -2,7 +2,7 @@ model = GModel ()
model:load ('square.geo')
model:load ('square.msh')
dg = dgSystemOfEquations (model)
-dg:setOrder(1)
+dg:setOrder(5)
-- conservation law
@@ -36,9 +36,9 @@ dg:L2Projection(createFunction.lua(1,'initial_condition','XYZ'))
dg:exportSolution('output/Advection_00000')
-- main loop
-for i=1,800 do
- norm = dg:RK44(0.01)
- if (i % 10 == 0) then
+for i=1,10000 do
+ norm = dg:RK44(0.001)
+ if (i % 50 == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/Advection-%05d", i))
end
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 5a15521e0b0d5c0a86cacd82dbd50f0f87ea2c90..718d2cd8f410289624f8f2bc67b8d169099d7423 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -122,7 +122,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
- const dgGroupOfFaces &group,
+ dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
fullMatrix<double> &solutionLeft,
fullMatrix<double> &solutionRight,
@@ -154,7 +154,8 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
dataCacheDouble &normals = cacheMapLeft.provideData("Normals");
dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
dataCacheElement &cacheElementRight = cacheMapRight.getElement();
- cacheMapLeft.provideData("UVW").set(group.getIntegrationPointsMatrix());
+ dataCacheDouble &uvwLeft = cacheMapLeft.provideData("UVW");
+ dataCacheDouble &uvwRight = cacheMapRight.provideData("UVW");
dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution");
dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("GradientSolution");
@@ -185,6 +186,8 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
dPsiRightDx.setAsProxy(DPsiLeftDx,iFace*nbNodesRight,nbNodesRight);
dofsLeft.setAsProxy(solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
dofsRight.setAsProxy(solutionRight, nbFields*group.getElementRightId(iFace), nbFields);
+ uvwLeft.setAsProxy(group.getIntegrationOnElementLeft(iFace));
+ uvwRight.setAsProxy(group.getIntegrationOnElementRight(iFace));
// proxies for the flux
normalFluxQP.setAsProxy(NormalFluxQP, iFace*nbFields*2, nbFields*2);
@@ -297,7 +300,7 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
- const dgGroupOfFaces &group,
+ dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
@@ -313,7 +316,7 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
dataCacheMap cacheMapLeft(group.getNbIntegrationPoints());
// provided dataCache
- cacheMapLeft.provideData("UVW").set(group.getIntegrationPointsMatrix());
+ dataCacheDouble &uvw=cacheMapLeft.provideData("UVW");
dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
dataCacheDouble &normals = cacheMapLeft.provideData("Normals");
dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
@@ -326,6 +329,7 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
// ----- 2.3.1 --- provide the data to the cacheMap
solutionQPLeft.setAsProxy(solutionQP, iFace*nbFields, nbFields);
normals.setAsProxy(group.getNormals(),iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
+ uvw.setAsProxy(group.getIntegrationOnElementLeft(iFace));
cacheElementLeft.set(group.getElementLeft(iFace));
// ----- 2.3.2 --- compute fluxes
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 566cf906d0c268872dcebe91a3867f6284790bf9..fd9c321490ff310f3c859b259f24ae6de9738332 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -19,14 +19,14 @@ class dgAlgorithm {
fullMatrix<double> &residual);
static void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
- const dgGroupOfFaces & group,
+ dgGroupOfFaces & group,
const fullMatrix<double> &solution, // ! at face nodes
fullMatrix<double> &solutionRight,
fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual); // ! at face nodes
static void residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
- const dgGroupOfFaces &group,
+ dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index d0fd5b45e0880c9e19646c3484d2b48787809e28..e4c7d1cb6d77a306aee59fe0b6c58b0c90e083ab 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/Solver/dgConservationLaw.cpp
@@ -6,7 +6,6 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
class term : public dataCacheDouble {
dataCacheMap cacheMapRight; // new cacheMap to pass to the Riemann solver
dataCacheDouble &solutionRight;
- dataCacheDouble &solutionLeft;
dataCacheDouble &outsideValue;
dataCacheDouble *riemannSolver;
dgConservationLaw &_claw;
@@ -15,11 +14,11 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw.nbFields()),
cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
solutionRight(cacheMapRight.provideData("Solution")),
- solutionLeft(cacheMapLeft.get("Solution",this)),
outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
_claw(claw)
{
riemannSolver=_claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
+ riemannSolver->addMeAsDependencyOf(this);
}
void _eval() {
@@ -43,15 +42,10 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
class dgBoundaryCondition0Flux : public dgBoundaryCondition {
dgConservationLaw &_claw;
class term : public dataCacheDouble {
- dataCacheDouble &UVW;
- dgConservationLaw &_claw;
public:
- term(dgConservationLaw &claw,dataCacheMap &cacheMapLeft): UVW(cacheMapLeft.get("UVW",this)),_claw(claw) {}
+ term(dgConservationLaw &claw,dataCacheMap &cacheMapLeft):
+ dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw.nbFields()) {}
void _eval() {
- if(_value.size1()!=UVW().size1()){
- //adjust sizes
- _value = fullMatrix<double>(UVW().size1(),_claw.nbFields());
- }
}
};
public:
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index ce654004be7002c474f720a3a34b58ca59ccdaf3..b29bdb89ce85ef68b149c071beeff8ff42428cac 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -8,16 +8,14 @@
class dgConservationLawAdvection : public dgConservationLaw {
std::string _vFunctionName;
class advection : public dataCacheDouble {
- dataCacheDouble &sol, &uvw, &v;
+ dataCacheDouble &sol, &v;
public:
advection(std::string vFunctionName, dataCacheMap &cacheMap):
- uvw(cacheMap.get("UVW",this)),
+ dataCacheDouble(cacheMap.getNbEvaluationPoints(),3),
sol(cacheMap.get("Solution",this)),
v(cacheMap.get(vFunctionName,this))
{};
void _eval () {
- if(_value.size1() != sol().size1())
- _value=fullMatrix<double>(sol().size1(),3);
for(int i=0; i< sol().size1(); i++) {
_value(i,0) = sol(i,0)*v(i,0);
_value(i,1) = sol(i,0)*v(i,1);
@@ -26,18 +24,16 @@ class dgConservationLawAdvection : public dgConservationLaw {
}
};
class riemann : public dataCacheDouble {
- dataCacheDouble &uvw, &normals, &solLeft, &solRight,&v;
+ dataCacheDouble &normals, &solLeft, &solRight,&v;
public:
riemann(std::string vFunctionName, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
- uvw(cacheMapLeft.get("UVW", this)),
+ dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),2),
normals(cacheMapLeft.get("Normals", this)),
solLeft(cacheMapLeft.get("Solution", this)),
solRight(cacheMapRight.get("Solution", this)),
v(cacheMapLeft.get(vFunctionName,this))
{};
void _eval () {
- if(_value.size1() !=solLeft().size1())
- _value = fullMatrix<double>(solLeft().size1(),2);
for(int i=0; i< _value.size1(); i++) {
double un=v(i,0)*normals(0,i)+v(i,1)*normals(1,i)+v(i,2)*normals(2,i);
if(un>0){
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index fdeed3c2e28d829eef5f211aa11a0e4e74831737..e3bd1a6fdcc56352243ad6cf4247eee51db38932 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -21,24 +21,24 @@ static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
return m;
}
-static fullMatrix<double> * dgGetFaceIntegrationRuleOnElement (
+static fullMatrix<double> dgGetFaceIntegrationRuleOnElement (
const polynomialBasis *fsFace,
const fullMatrix<double> &intgFace,
const polynomialBasis *fsElement,
- const std::vector <int> *closure) {
+ const std::vector <int> &closure) {
int npts=intgFace.size1();
- fullMatrix<double> *m = new fullMatrix<double> (npts, 4);
+ fullMatrix<double> intgElement(npts, 4);
double f[256];
for (int i=0;i<npts;i++){
fsFace->f(intgFace(i,0),intgFace(i,1),intgFace(i,2),f);
- for(size_t j=0; j<closure->size();j++){
- int jNod=(*closure)[j];
+ for(size_t j=0; j<closure.size();j++){
+ int jNod=closure[j];
for(int k=0;k<fsElement->points.size2();k++)
- (*m)(i,k) += f[j] * fsElement->points(jNod,k);
- (*m)(i,3) = intgFace(i,3);
+ intgElement(i,k) += f[j] * fsElement->points(jNod,k);
+ intgElement(i,3) = intgFace(i,3);
}
}
- return m;
+ return intgElement;
}
@@ -204,6 +204,10 @@ void dgGroupOfFaces::init(int pOrder) {
_redistribution = new fullMatrix<double> (_fsFace->coefficients.size1(),_integration->size1());
_collocation = new fullMatrix<double> (_integration->size1(), _fsFace->coefficients.size1());
_detJac = new fullMatrix<double> (_integration->size1(), _faces.size());
+ for (size_t i=0; i<_closuresLeft.size(); i++)
+ _integrationPointsLeft.push_back(dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,_closuresLeft[i]));
+ for (size_t i=0; i<_closuresRight.size(); i++)
+ _integrationPointsRight.push_back(dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsRight,_closuresRight[i]));
double f[256];
for (int j=0;j<_integration->size1();j++) {
_fsFace->f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
@@ -231,11 +235,11 @@ void dgGroupOfFaces::init(int pOrder) {
int index = 0;
for (size_t i=0; i<_faces.size();i++){
const std::vector<int> &closureLeft=getClosureLeft(i);
- fullMatrix<double> *intLeft=dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,&closureLeft);
+ fullMatrix<double> &intLeft=_integrationPointsLeft[_closuresIdLeft[i]];
double jac[3][3],ijac[3][3];
- for (int j=0; j<intLeft->size1(); j++){
- _fsLeft->df((*intLeft)(j,0),(*intLeft)(j,1),(*intLeft)(j,2),g);
- getElementLeft(i)->getJacobian ((*intLeft)(j,0), (*intLeft)(j,1), (*intLeft)(j,2), jac);
+ for (int j=0; j<intLeft.size1(); j++){
+ _fsLeft->df((intLeft)(j,0),(intLeft)(j,1),(intLeft)(j,2),g);
+ getElementLeft(i)->getJacobian ((intLeft)(j,0), (intLeft)(j,1), (intLeft)(j,2), jac);
inv3x3(jac,ijac);
//compute dPsiLeftDxOnQP
//(iPsi*3+iXYZ,iQP+iFace*NQP);
@@ -264,13 +268,12 @@ void dgGroupOfFaces::init(int pOrder) {
nz/=norm;
index++;
}
- delete intLeft;
// there is nothing on the right for boundary groups
if(_fsRight){
- fullMatrix<double> *intRight=dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsRight,&getClosureRight(i));
- for (int j=0; j<intRight->size1(); j++){
- _fsRight->df((*intRight)(j,0),(*intRight)(j,1),(*intRight)(j,2),g);
- getElementRight(i)->getJacobian ((*intRight)(j,0), (*intRight)(j,1), (*intRight)(j,2), jac);
+ fullMatrix<double> &intRight=_integrationPointsLeft[_closuresIdRight[i]];
+ for (int j=0; j<intRight.size1(); j++){
+ _fsRight->df((intRight)(j,0),(intRight)(j,1),(intRight)(j,2),g);
+ getElementRight(i)->getJacobian ((intRight)(j,0), (intRight)(j,1), (intRight)(j,2), jac);
inv3x3(jac,ijac);
//compute dPsiRightDxOnQP
// (iQP*3+iXYZ , iFace*NPsi+iPsi)
@@ -281,7 +284,6 @@ void dgGroupOfFaces::init(int pOrder) {
(*_dPsiRightDxOnQP)(j*3+2,i*nPsi+iPsi) = g[iPsi][0]*ijac[2][0]+g[iPsi][1]*ijac[2][1]+g[iPsi][2]*ijac[2][2];
}
}
- delete intRight;
}
}
}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 2836e546e7e5bc8c76787c072525686c61a0ed50..793eb98bf04f9babf9cde73cf9b360f0bf27a2c6 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -111,6 +111,9 @@ class dgGroupOfFaces {
std::vector<std::vector<int> > _closuresRight;
std::vector<int> _closuresIdLeft;
std::vector<int> _closuresIdRight;
+ // face integration point in the coordinate of the left and right element (one fullMatrix per closure)
+ std::vector<fullMatrix<double> > _integrationPointsLeft;
+ std::vector<fullMatrix<double> > _integrationPointsRight;
// XYZ gradient of the shape functions of both elements on the integrations points of the face
// (iQP*3+iXYZ , iFace*NPsi+iPsi)
fullMatrix<double> *_dPsiLeftDxOnQP;
@@ -134,8 +137,11 @@ public:
inline MElement* getElementLeft (int i) const {return _groupLeft.getElement(_left[i]);}
inline MElement* getElementRight (int i) const {return _groupRight.getElement(_right[i]);}
inline MElement* getFace (int iElement) const {return _faces[iElement];}
- const std::vector<int> &getClosureLeft(int iFace) const{ return _closuresLeft[_closuresIdLeft[iFace]];}
- const std::vector<int> &getClosureRight(int iFace) const{ return _closuresRight[_closuresIdRight[iFace]];}
+ inline const std::vector<int> &getClosureLeft(int iFace) const{ return _closuresLeft[_closuresIdLeft[iFace]];}
+ inline const std::vector<int> &getClosureRight(int iFace) const{ return _closuresRight[_closuresIdRight[iFace]];}
+ inline fullMatrix<double> &getIntegrationOnElementLeft(int iFace) { return _integrationPointsLeft[_closuresIdLeft[iFace]];}
+ inline fullMatrix<double> &getIntegrationOnElementRight(int iFace) { return _integrationPointsRight[_closuresIdRight[iFace]];}
+
inline fullMatrix<double> &getNormals () const {return *_normals;}
dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices);
diff --git a/Solver/function.h b/Solver/function.h
index c9d2b6ba32e553e43486dace44b90504a1296e45..89d775ebc23070aa2820942d150adec34aa40014 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -42,11 +42,11 @@ protected :
it!=_dependOnMe.end(); it++)
**it=false;
}
- // dataCacheMap is the only one supposed to call this
- void addMeAsDependencyOf (dataCache *newDep);
dataCache() : _valid(false) {}
virtual ~dataCache(){};
public :
+ // dataCacheMap is the only one supposed to call this
+ void addMeAsDependencyOf (dataCache *newDep);
inline bool somethingDependOnMe() {
return !_dependOnMe.empty();
}
@@ -72,6 +72,13 @@ class dataCacheDouble : public dataCache {
_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(fullMatrix<double> &mat) {
+ _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