From 369ed8b483581a4da8a51f477cb21fa17c1615f2 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Fri, 26 Feb 2010 10:40:00 +0000
Subject: [PATCH] resize dataCacheMap automatically
---
Solver/TESTCASES/ForwardFacingStep.lua | 2 +-
Solver/dgAlgorithm.cpp | 52 +++++++++----------
Solver/dgConservationLaw.cpp | 48 ++++++++----------
Solver/dgConservationLawAdvection.cpp | 8 ++-
Solver/dgConservationLawPerfectGas.cpp | 58 +++++-----------------
Solver/dgConservationLawShallowWater2d.cpp | 12 ++---
Solver/dgConservationLawWaveEquation.cpp | 17 ++-----
Solver/dgDofContainer.cpp | 5 +-
Solver/dgGroupOfElements.cpp | 4 --
Solver/dgLimiter.cpp | 5 +-
Solver/function.cpp | 32 ++++++++----
Solver/function.h | 30 ++++++++---
Solver/luaFunction.cpp | 2 +-
13 files changed, 122 insertions(+), 153 deletions(-)
diff --git a/Solver/TESTCASES/ForwardFacingStep.lua b/Solver/TESTCASES/ForwardFacingStep.lua
index 4c9165ce23..0507416ca7 100644
--- a/Solver/TESTCASES/ForwardFacingStep.lua
+++ b/Solver/TESTCASES/ForwardFacingStep.lua
@@ -41,7 +41,7 @@ FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
law=dgPerfectGasLaw2d()
DG:setConservationLaw(law)
-law:addBoundaryCondition('Walls',law:newWallBoundary())
+law:addBoundaryCondition('Walls',law:newBoundaryWall())
law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS))
--law:addBoundaryCondition('Walls',law:newOutsideValueBoundary(FS))
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 99fc7bce66..5521ba2b01 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -43,13 +43,13 @@ 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(group.getNbIntegrationPoints());
+ dataCacheMap cacheMap;
+ cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
dataCacheElement &cacheElement = cacheMap.getElement();
// provided dataCache
- cacheMap.provideData("UVW").set(group.getIntegrationPointsMatrix());
- dataCacheDouble &solutionQPe = cacheMap.provideData("Solution");
- dataCacheDouble &gradientSolutionQPe = cacheMap.provideData("SolutionGradient");
- gradientSolutionQPe.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
+ cacheMap.provideData("UVW",1,3).set(group.getIntegrationPointsMatrix());
+ dataCacheDouble &solutionQPe = cacheMap.provideData("Solution",1,nbFields);
+ dataCacheDouble &gradientSolutionQPe = cacheMap.provideData("SolutionGradient",3,nbFields);
// needed dataCache
dataCacheDouble *sourceTerm = claw.newSourceTerm(cacheMap);
// fluxes in XYZ coordinates;
@@ -142,23 +142,22 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), nbFaces*nbFields*2);
// create one dataCache for each side
- dataCacheMap cacheMapLeft(group.getNbIntegrationPoints());
- dataCacheMap cacheMapRight(group.getNbIntegrationPoints());
+ dataCacheMap cacheMapLeft;
+ cacheMapLeft.setNbEvaluationPoints(group.getNbIntegrationPoints());
+ dataCacheMap cacheMapRight;
+ cacheMapRight.setNbEvaluationPoints(group.getNbIntegrationPoints());
// data this algorithm provide to the cache map (we can maybe move each data to a separate function but
// I think It's easier like this)
- dataCacheDouble &normals = cacheMapLeft.provideData("Normals");
+ dataCacheDouble &normals = cacheMapLeft.provideData("Normals",1,1);
dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
dataCacheElement &cacheElementRight = cacheMapRight.getElement();
- dataCacheDouble &uvwLeft = cacheMapLeft.provideData("UVW");
- dataCacheDouble &uvwRight = cacheMapRight.provideData("UVW");
- dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
- dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution");
- dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient");
- dataCacheDouble &gradientSolutionRight = cacheMapRight.provideData("SolutionGradient");
- //resize gradientSolutionLeft and Right
- gradientSolutionLeft.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
- gradientSolutionRight.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
+ dataCacheDouble &uvwLeft = cacheMapLeft.provideData("UVW",1,3);
+ dataCacheDouble &uvwRight = cacheMapRight.provideData("UVW",1,3);
+ dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution",1,nbFields);
+ dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution",1,nbFields);
+ dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient",3,nbFields);
+ dataCacheDouble &gradientSolutionRight = cacheMapRight.provideData("SolutionGradient",3,nbFields);
// A2 ) ask the law to provide the fluxes (possibly NULL)
dataCacheDouble *riemannSolver = claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
@@ -425,15 +424,15 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
// ----- 2 ---- compute normal fluxes at integration points
fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields);
- dataCacheMap cacheMapLeft(group.getNbIntegrationPoints());
+ dataCacheMap cacheMapLeft;
+ cacheMapLeft.setNbEvaluationPoints(group.getNbIntegrationPoints());
fullMatrix<double> &DPsiLeftDx = group.getDPsiLeftDxMatrix();
// provided dataCache
- dataCacheDouble &uvw=cacheMapLeft.provideData("UVW");
- dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
- dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient");
- dataCacheDouble &normals = cacheMapLeft.provideData("Normals");
+ dataCacheDouble &uvw=cacheMapLeft.provideData("UVW",1,3);
+ dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution",1,nbFields);
+ dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient",3,nbFields);
+ dataCacheDouble &normals = cacheMapLeft.provideData("Normals",1,1);
dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
- gradientSolutionLeft.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
// required data
// inviscid
dataCacheDouble *boundaryTerm = boundaryCondition->newBoundaryTerm(cacheMapLeft);
@@ -553,14 +552,15 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
std::vector<double> & DT
)
{
- dataCacheMap cacheMap(group.getNbNodes());
- dataCacheDouble &sol = cacheMap.provideData("Solution");
+ const int nbFields = claw.getNbFields();
+ dataCacheMap cacheMap;
+ cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
+ dataCacheDouble &sol = cacheMap.provideData("Solution",1,nbFields);
dataCacheElement &cacheElement = cacheMap.getElement();
// provided dataCache
dataCacheDouble *maxConvectiveSpeed = claw.newMaxConvectiveSpeed(cacheMap);
dataCacheDouble *maximumDiffusivity = claw.newMaximumDiffusivity(cacheMap);
- const int nbFields = claw.getNbFields();
/* This is an estimate on how lengths changes with p
It is merely the smallest distance between gauss
points at order p + 1 */
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index eb6c83db84..85f74c9f60 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/Solver/dgConservationLaw.cpp
@@ -10,12 +10,12 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dgConservationLaw *_claw;
public:
term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
- dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()),
- cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
- solutionRight(cacheMapRight.provideData("Solution")),
+ dataCacheDouble(cacheMapLeft, 1, claw->getNbFields()),
+ solutionRight(cacheMapRight.provideData("Solution",1,claw->getNbFields())),
outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
_claw(claw)
{
+ cacheMapRight.setNbEvaluationPoints(cacheMapLeft.getNbEvaluationPoints());
riemannSolver=_claw->newRiemannSolver(cacheMapLeft,cacheMapRight);
if(riemannSolver)
riemannSolver->addMeAsDependencyOf(this);
@@ -34,7 +34,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dataCacheDouble &outsideValue;
public:
dirichlet(dgConservationLaw *claw, dataCacheMap &cacheMap,const std::string outsideValueFunctionName):
- dataCacheDouble(cacheMap, cacheMap.getNbEvaluationPoints(),claw->getNbFields()),
+ dataCacheDouble(cacheMap, 1, claw->getNbFields()),
outsideValue(cacheMap.get(outsideValueFunctionName,this)){}
void _eval () {
for(int i=0;i<_value.size1();i++)
@@ -50,12 +50,12 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dgConservationLaw *_claw;
public:
maximumDiffusivity(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
- dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()),
- cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
- solutionRight(cacheMapRight.provideData("Solution")),
+ dataCacheDouble(cacheMapLeft, 1, claw->getNbFields()),
+ solutionRight(cacheMapRight.provideData("Solution",1,claw->getNbFields())),
outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
_claw(claw)
{
+ cacheMapRight.setNbEvaluationPoints(cacheMapLeft.getNbEvaluationPoints());
maxDif = _claw->newMaximumDiffusivity(cacheMapRight);
if(maxDif)
maxDif->addMeAsDependencyOf(this);
@@ -89,7 +89,7 @@ class dgBoundaryConditionNeumann : public dgBoundaryCondition {
dataCacheDouble &flux;
public:
term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string fluxFunctionName):
- dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()),
+ dataCacheDouble(cacheMapLeft,1 ,claw->getNbFields()),
flux(cacheMapLeft.get(fluxFunctionName,this))
{}
void _eval() {
@@ -113,7 +113,7 @@ class dgBoundarySymmetry : public dgBoundaryCondition {
dgConservationLaw *_claw;
public:
term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft):
- dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()), _claw(claw)
+ dataCacheDouble(cacheMapLeft, 1,claw->getNbFields()), _claw(claw)
{
riemannSolver=_claw->newRiemannSolver(cacheMapLeft,cacheMapLeft);
riemannSolver->addMeAsDependencyOf(this);
@@ -138,7 +138,7 @@ class dgBoundaryCondition0Flux : public dgBoundaryCondition {
class term : public dataCacheDouble {
public:
term(dgConservationLaw *claw,dataCacheMap &cacheMapLeft):
- dataCacheDouble(cacheMapLeft,cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()) {}
+ dataCacheDouble(cacheMapLeft,1,claw->getNbFields()) {}
void _eval() {
}
};
@@ -167,15 +167,12 @@ class dgBoundaryCondition::dirichlet_ : public dataCacheDouble {
dataCacheDouble /
public:
dirichlet_(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,cacheMap.get("Solution")().size2()),
sol(cacheMap.get("Solution",this)){}
void _eval () {
- int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,sol().size2());
- for(int i=0; i< nQP; i++)
+ for(int i=0; i< _value.size1(); i++)
for (int k=0;k<sol().size2();k++)
- _value(i,k) = sol(i,k);
+ _value(i,k) = sol(i,k);
}
};
@@ -185,7 +182,7 @@ class dgBoundaryCondition::neumann_ : public dataCacheDouble {
dataCacheDouble *diffusiveFlux;
public:
neumann_(dataCacheMap &cacheMap, dgConservationLaw *claw):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,cacheMap.get("Solution")().size2()),
_claw (claw),
sol(cacheMap.get("Solution",this)),
normals(cacheMap.get("Normals",this)){
@@ -193,18 +190,13 @@ public:
if (diffusiveFlux)diffusiveFlux->addMeAsDependencyOf(this);
}
void _eval () {
- int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,sol().size2());
-
const fullMatrix<double> &dfl = (*diffusiveFlux)();
-
- for(int i=0; i< nQP; i++) {
- for (int k=0;k<sol().size2();k++) {
- _value(i,k) =
- dfl(i,k+sol().size2()*0) *normals(0,i) +
- dfl(i,k+sol().size2()*1) *normals(1,i) +
- dfl(i,k+sol().size2()*2) *normals(2,i);
+ for(int i=0; i< _value.size1(); i++) {
+ for (int k=0;k<_value.size2();k++) {
+ _value(i,k) =
+ dfl(i,k+sol().size2()*0) *normals(0,i) +
+ dfl(i,k+sol().size2()*1) *normals(1,i) +
+ dfl(i,k+sol().size2()*2) *normals(2,i);
}
}
}
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index 45e6176814..06674dd29f 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -10,7 +10,7 @@ class dgConservationLawAdvectionDiffusion::advection : public dataCacheDouble {
dataCacheDouble &sol, &v;
public:
advection(std::string vFunctionName, dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap,cacheMap.getNbEvaluationPoints(),3),
+ dataCacheDouble(cacheMap,1,3),
sol(cacheMap.get("Solution",this)),
v(cacheMap.get(vFunctionName,this))
{};
@@ -26,7 +26,7 @@ class dgConservationLawAdvectionDiffusion::riemann : public dataCacheDouble {
dataCacheDouble &normals, &solLeft, &solRight,&v;
public:
riemann(std::string vFunctionName, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
- dataCacheDouble(cacheMapLeft,cacheMapLeft.getNbEvaluationPoints(),2),
+ dataCacheDouble(cacheMapLeft,1,2),
normals(cacheMapLeft.get("Normals", this)),
solLeft(cacheMapLeft.get("Solution", this)),
solRight(cacheMapRight.get("Solution", this)),
@@ -49,13 +49,11 @@ class dgConservationLawAdvectionDiffusion::diffusion : public dataCacheDouble {
dataCacheDouble &solgrad, ν
public:
diffusion(std::string nuFunctionName, dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,3),
solgrad(cacheMap.get("SolutionGradient",this)),
nu(cacheMap.get(nuFunctionName,this))
{};
void _eval () {
- if(_value.size1() != solgrad().size1()/3)
- _value=fullMatrix<double>(solgrad().size1()/3,3);
for(int i=0; i< solgrad().size1()/3; i++) {
_value(i,0) = -solgrad(i*3,0)*nu(i,0);
_value(i,1) = -solgrad(i*3+1,0)*nu(i,0);
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index 9eac309d43..1622dc7c81 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -328,13 +328,11 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
dataCacheDouble /
public:
advection(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,8),
sol(cacheMap.get("Solution",this))
{};
void _eval () {
const int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,8);
const double GM1 = GAMMA - 1.0;
for (size_t k = 0 ; k < nQP; k++ ){
const double invrho = 1./sol(k,0);
@@ -366,7 +364,7 @@ class dgPerfectGasLaw2d::diffusion : public dataCacheDouble {
diffusion(dataCacheMap &cacheMap,
const std::string &muFunctionName,
const std::string &kappaFunctionName ):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,12),
sol(cacheMap.get("Solution",this)),
grad(cacheMap.get("SolutionGradient",this)),
mu(cacheMap.get(muFunctionName,this)),
@@ -374,9 +372,6 @@ class dgPerfectGasLaw2d::diffusion : public dataCacheDouble {
{};
void _eval () {
const int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,12);
-
for (size_t k = 0 ; k < nQP; k++ ){
const double muk = mu(k,0);
const double kappak = kappa(k,0);
@@ -427,14 +422,12 @@ class dgPerfectGasLaw2d::source : public dataCacheDouble {
dataCacheDouble &sol,&s;
public:
source(dataCacheMap &cacheMap, const std::string &sourceFunctionName):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,4),
sol(cacheMap.get("Solution",this)),
s(cacheMap.get(sourceFunctionName,this))
{};
void _eval () {
const int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,4);
for (size_t k = 0 ; k < nQP; k++ ){
_value(k,0) = sol(k,0)*s(0,0);
_value(k,1) = sol(k,0)*s(0,1);
@@ -449,12 +442,11 @@ class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
double _presMin, _rhoMin;
public:
clipToPhysics(dataCacheMap &cacheMap, double presMin, double rhoMin):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,4),
sol(cacheMap.get("Solution",this))
{
_presMin=presMin;
_rhoMin=rhoMin;
- _value=fullMatrix<double>(cacheMap.getNbEvaluationPoints(),4);
};
void _eval () {
const int nDofs = sol().size1();
@@ -483,16 +475,13 @@ class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
public:
riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
- dataCacheDouble(cacheMapLeft),
+ dataCacheDouble(cacheMapLeft,1,12),
normals(cacheMapLeft.get("Normals", this)),
solL(cacheMapLeft.get("Solution", this)),
solR(cacheMapRight.get("Solution", this))
{};
void _eval () {
int nQP = solL().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,12);
-
for(int i=0; i< nQP; i++) {
const double solLeft [4] = {solL(i,0),solL(i,1),solL(i,2),solL(i,3)};
const double solRight[4] = {solR(i,0),solR(i,1),solR(i,2),solR(i,3)};
@@ -517,16 +506,13 @@ class dgPerfectGasLaw2d::riemannGodunov : public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
public:
riemannGodunov(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
- dataCacheDouble(cacheMapLeft),
+ dataCacheDouble(cacheMapLeft,1,12),
normals(cacheMapLeft.get("Normals", this)),
solL(cacheMapLeft.get("Solution", this)),
solR(cacheMapRight.get("Solution", this))
{};
void _eval () {
int nQP = solL().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,12);
-
for(int i=0; i< nQP; i++) {
const double nx = normals(0,i);
const double ny = normals(1,i);
@@ -584,14 +570,12 @@ class dgPerfectGasLaw2d::maxConvectiveSpeed : public dataCacheDouble {
dataCacheDouble /
public:
maxConvectiveSpeed(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,1),
sol(cacheMap.get("Solution",this))
{
};
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,1);
for(int i=0; i< nQP; i++) {
const double rhov2 = (sol(i,1)*sol(i,1) + sol(i,2)*sol(i,2))/sol(i,0);
const double p = (GAMMA-1.0)*(sol(i,3) - 0.5 * rhov2);
@@ -606,7 +590,7 @@ class dgPerfectGasLaw2d::maxDiffusivity : public dataCacheDouble {
dataCacheDouble &sol,&mu,κ
public:
maxDiffusivity(dataCacheMap &cacheMap, const std::string &muFunctionName, const std::string &kappaFunctionName ):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,1),
sol(cacheMap.get("Solution",this)),
mu(cacheMap.get(muFunctionName,this)),
kappa(cacheMap.get(kappaFunctionName,this))
@@ -614,8 +598,6 @@ class dgPerfectGasLaw2d::maxDiffusivity : public dataCacheDouble {
};
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,1);
for(int i=0; i< nQP; i++) {
_value(i,0) = std::max(mu(i,0),kappa(i,0));
}
@@ -670,14 +652,11 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
// NON VISCOUS BOUNDARY FLUX
//-------------------------------------------------------------------------------
term(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,4),
sol(cacheMap.get("Solution",this)),
normals(cacheMap.get("Normals",this)){}
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,4);
-
for(int i=0; i< nQP; i++) {
const double nx = normals(0,i);
const double ny = normals(1,i);
@@ -718,13 +697,10 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
dataCacheDouble /
public:
dirichletNonSlip(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,4),
sol(cacheMap.get("Solution",this)){}
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,4);
-
for(int i=0; i< nQP; i++) {
_value(i,0) = sol(i,0);
_value(i,1) = 0.0;
@@ -746,7 +722,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
dataCacheDouble *diffusiveFlux;
public:
neumannNonSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,4),
_claw (claw),
sol(cacheMap.get("Solution",this)),
normals(cacheMap.get("Normals",this)){
@@ -755,9 +731,6 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
}
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,4);
-
const fullMatrix<double> &dfl = (*diffusiveFlux)();
for(int i=0; i< nQP; i++) {
@@ -782,13 +755,10 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
dataCacheDouble /
public:
dirichletSlip(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,cacheMap.get("Solution")().size2()),
sol(cacheMap.get("Solution",this)){}
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,sol().size2());
-
for(int i=0; i< nQP; i++) {
for(int k=0; k< sol().size2(); k++) {
_value(i,k) = sol(i,k);
@@ -807,14 +777,12 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
public:
neumannSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
_claw (claw),
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,4),
sol(cacheMap.get("Solution",this)),
normals(cacheMap.get("Normals",this)){
}
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,4);
_value.setAll(0.0);
}
};
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index 4873c3b9d0..572c219ac2 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -8,11 +8,11 @@ class dgConservationLawShallowWater2d::advection: public dataCacheDouble {
dataCacheDouble /
public:
advection(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap,cacheMap.getNbEvaluationPoints(),9),
+ dataCacheDouble(cacheMap,1,9),
sol(cacheMap.get("Solution",this))
{};
void _eval () {
- int nQP = sol().size1();
+ int nQP = _value.size1();
for(int i=0; i< nQP; i++) {
double eta = sol(i,0);
double u = sol(i,1);
@@ -37,7 +37,7 @@ class dgConservationLawShallowWater2d::source: public dataCacheDouble {
dataCacheDouble &xyz, &solution,&solutionGradient;
public :
source(dataCacheMap &cacheMap) :
- dataCacheDouble(cacheMap,cacheMap.getNbEvaluationPoints(),3),
+ dataCacheDouble(cacheMap,1,3),
xyz(cacheMap.get("XYZ",this)),
solution(cacheMap.get("Solution",this)),
solutionGradient(cacheMap.get("SolutionGradient",this))
@@ -64,7 +64,7 @@ class dgConservationLawShallowWater2d::riemann:public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
public:
riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
- dataCacheDouble(cacheMapLeft,cacheMapLeft.getNbEvaluationPoints(),6),
+ dataCacheDouble(cacheMapLeft,1,6),
normals(cacheMapLeft.get("Normals", this)),
solL(cacheMapLeft.get("Solution", this)),
solR(cacheMapRight.get("Solution", this))
@@ -133,13 +133,11 @@ class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition
dataCacheDouble &sol,&normals;
public:
term(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,3),
sol(cacheMap.get("Solution",this)),
normals(cacheMap.get("Normals",this)){}
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,3);
for(int i=0; i< nQP; i++) {
double nx = normals(0,i);
double ny = normals(1,i);
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index 9a8c9d8c63..370bb2508b 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -10,13 +10,11 @@ class dgConservationLawWaveEquation::hyperbolicFlux : public dataCacheDouble {
const int _DIM,_nbf;
public:
hyperbolicFlux(dataCacheMap &cacheMap,int DIM):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,3*(DIM+1)),
sol(cacheMap.get("Solution",this)),_DIM(DIM),_nbf(DIM+1)
{};
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,3*_nbf);
_value.scale(0.);
for(int i=0; i< nQP; i++) {
const double p = sol(i,0);
@@ -32,14 +30,11 @@ class dgConservationLawWaveEquation::maxConvectiveSpeed : public dataCacheDouble
dataCacheDouble /
public:
maxConvectiveSpeed(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,1),
sol(cacheMap.get("Solution",this))
{
};
void _eval () {
- int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value=fullMatrix<double>(nQP,1);
_value.setAll(c);
}
};
@@ -49,7 +44,7 @@ class dgConservationLawWaveEquation::riemann : public dataCacheDouble {
const int _DIM,_nbf;
public:
riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, int DIM):
- dataCacheDouble(cacheMapLeft),
+ dataCacheDouble(cacheMapLeft,1,2*(DIM+1)),
normals(cacheMapLeft.get("Normals", this)),
solL(cacheMapLeft.get("Solution", this)),
solR(cacheMapRight.get("Solution", this)),
@@ -57,8 +52,6 @@ class dgConservationLawWaveEquation::riemann : public dataCacheDouble {
{};
void _eval () {
int nQP = solL().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,2*_nbf);
for(int i=0; i< nQP; i++) {
const double n[3] = {normals(0,i),normals(1,i),normals(2,i)};
double unL=0,unR=0;
@@ -110,14 +103,12 @@ class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition {
dataCacheDouble &sol,&normals;
public:
term(dataCacheMap &cacheMap, int DIM):
- dataCacheDouble(cacheMap),
+ dataCacheDouble(cacheMap,1,DIM+1),
sol(cacheMap.get("Solution",this)),
normals(cacheMap.get("Normals",this)),
_DIM(DIM){}
void _eval () {
int nQP = sol().size1();
- if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,_DIM+1);
for(int i=0; i< nQP; i++) {
const double n[3] = {normals(0,i),normals(1,i),normals(2,i)};
double p = sol(i,0);
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 7eb8fa5c82..7a2b0d0cf6 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -193,9 +193,10 @@ void dgDofContainer::L2Projection(std::string functionName){
for (int iGroup=0;iGroup<_groups.getNbElementGroups();iGroup++) {
const dgGroupOfElements &group = *_groups.getElementGroup(iGroup);
fullMatrix<double> Source = fullMatrix<double> (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
- dataCacheMap cacheMap(group.getNbIntegrationPoints());
+ dataCacheMap cacheMap;
+ cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
dataCacheElement &cacheElement = cacheMap.getElement();
- cacheMap.provideData("UVW").set(group.getIntegrationPointsMatrix());
+ cacheMap.provideData("UVW",1,3).set(group.getIntegrationPointsMatrix());
dataCacheDouble &sourceTerm = cacheMap.get(functionName);
fullMatrix<double> source;
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 551752deec..824d280832 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -674,8 +674,6 @@ void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
}
}
-
- Msg::Info("%d groups of elements",localElements.size());
// do a group of elements for every element type in the mesh
for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
_elementGroups.push_back(new dgGroupOfElements(it->second,order));
@@ -911,8 +909,6 @@ dgGroupCollection::dgGroupCollection(GModel *model, int dimension, int order)
buildGroupsOfElements(model,dimension,order);
}
-
-
dgGroupCollection::~dgGroupCollection() {
for (int i=0; i< _elementGroups.size(); i++)
delete _elementGroups[i];
diff --git a/Solver/dgLimiter.cpp b/Solver/dgLimiter.cpp
index 45c7b55e99..6b9c9fb294 100644
--- a/Solver/dgLimiter.cpp
+++ b/Solver/dgLimiter.cpp
@@ -107,8 +107,9 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution)
dgGroupOfElements* egroup = groups->getElementGroup(iG);
fullMatrix<double> &solGroup = solution->getGroupProxy(iG);
- dataCacheMap cacheMap(egroup->getNbNodes());//nbdofs for each element
- dataCacheDouble &solutionE = cacheMap.provideData("Solution");
+ dataCacheMap cacheMap;
+ cacheMap.setNbEvaluationPoints(egroup->getNbNodes());//nbdofs for each element
+ dataCacheDouble &solutionE = cacheMap.provideData("Solution",1,nbFields);
dataCacheElement &cacheElement = cacheMap.getElement();
dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
if (solutionEClipped){
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 15afbfa243..62b88d983d 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -68,12 +68,12 @@ dataCacheDouble &dataCacheMap::get(const std::string &functionName, dataCache *c
return *r;
}
-dataCacheDouble &dataCacheMap::provideData(std::string name)
+dataCacheDouble &dataCacheMap::provideData(std::string name,int nRowByPoints, int nCol)
{
- dataCacheDouble *&r= _cacheDoubleMap[name];
- if(r!=NULL)
+ if (_cacheDoubleMap[name] != NULL)
throw;
- r = new providedDataDouble(*this);
+ dataCacheDouble *r = new providedDataDouble(*this,nRowByPoints, nCol);
+ _cacheDoubleMap[name] = r;
return *r;
}
@@ -97,7 +97,7 @@ class functionXYZ : public function {
int count;
public:
data(dataCacheMap *m) :
- dataCacheDouble(*m, m->getNbEvaluationPoints(),3),
+ dataCacheDouble(*m, 1,3),
_element(m->getElement(this)), _uvw(m->get("UVW", this))
{
}
@@ -136,7 +136,7 @@ class functionStructuredGridFile : public function {
functionStructuredGridFile *_fun;
public:
data(functionStructuredGridFile *fun, dataCacheMap *m):
- dataCacheDouble(*m,m->getNbEvaluationPoints(),1),
+ dataCacheDouble(*m,1,1),
coord(m->get(fun->_coordFunction,this))
{
_fun=fun;
@@ -204,7 +204,7 @@ class functionConstant::data : public dataCacheDouble {
const functionConstant *_function;
public:
data(const functionConstant * function,dataCacheMap *m):
- dataCacheDouble(*m,m->getNbEvaluationPoints(),function->_source.size1()){
+ dataCacheDouble(*m,1,function->_source.size1()){
_function = function;
}
void _eval() {
@@ -247,13 +247,11 @@ public:
data(dataCacheMap &m, dgDofContainer *sys) :
_dofContainer(sys),
xyz(m.get("XYZ",this)),
- dataCacheDouble(m,m.getNbEvaluationPoints(), sys->getNbFields())
+ dataCacheDouble(m,1, sys->getNbFields())
{
}
void _eval() {
int nP =xyz().size1();
- if(_value.size1() != nP)
- _value = fullMatrix<double>(nP,_dofContainer->getNbFields());
_value.setAll(0.0);
double fs[256];
fullMatrix<double> solEl;
@@ -282,14 +280,26 @@ public:
}
}
};
+void dataCacheMap::setNbEvaluationPoints(int nbEvaluationPoints) {
+ _nbEvaluationPoints = nbEvaluationPoints;
+ for(std::map<std::string, dataCacheDouble*>::iterator it = _cacheDoubleMap.begin(); it!= _cacheDoubleMap.end(); it++)
+ it->second->resize();
+}
dataCacheDouble *functionMesh2Mesh::newDataCache(dataCacheMap *m)
{
- printf("coussdo %d %d\n",m->getNbEvaluationPoints(),_dofContainer->getNbFields());
return new data(*m,_dofContainer);
}
+dataCacheDouble::dataCacheDouble(dataCacheMap &map,int nRowByPoint, int nCol):
+ dataCache(&map),_cacheMap(map),_value(nRowByPoint==0?1:nRowByPoint*map.getNbEvaluationPoints(),nCol){
+ _nRowByPoint=nRowByPoint;
+};
+void dataCacheDouble::resize() {
+ _value = fullMatrix<double>(_nRowByPoint==0?1:_nRowByPoint*_cacheMap.getNbEvaluationPoints(),_value.size2());
+}
+
#include "Bindings.h"
void function::registerDefaultFunctions()
diff --git a/Solver/function.h b/Solver/function.h
index 333dca348a..607135a53a 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -61,8 +61,12 @@ public :
};
-// dataCache when the value is a double
+// dataCache when the value is a matrix of double
+// the user should provide the number of rows by evaluating points and the number of columns
+// then the size of the matrix is automatically adjusted
class dataCacheDouble : public dataCache {
+ int _nRowByPoint;
+ dataCacheMap &_cacheMap;
protected:
fullMatrix<double> _value;
// do the actual computation and put the result into _value
@@ -113,8 +117,8 @@ class dataCacheDouble : public dataCache {
}
return _value;
}
- dataCacheDouble(dataCacheMap &map):dataCache(&map){};
- dataCacheDouble(dataCacheMap &map,int size1, int size2):dataCache(&map),_value(size1,size2){};
+ void resize();
+ dataCacheDouble(dataCacheMap &map,int nRowByPoints, int nCol);
virtual ~dataCacheDouble(){};
};
@@ -167,22 +171,32 @@ class dataCacheMap {
{
void _eval() {throw;};
public:
- providedDataDouble(dataCacheMap &map):dataCacheDouble(map) {
+ providedDataDouble(dataCacheMap &map, int nRowByPoints, int ncol):dataCacheDouble(map,nRowByPoints,ncol) {
_valid=true;
}
};
std::set<dataCache*> _toDelete;
+
protected:
void addDataCache(dataCache *data){
_toDelete.insert(data);
}
public:
+ // dg Part
+ //list of dgDofContainer that have to be evaluated, this not work as a function as it depends on the algorithm
+ std::map<dgDofContainer*,dataCacheDouble*> fields;
+ //list of dgDofContainer whom gradient are needed
+ std::map<dgDofContainer*,dataCacheDouble*> gradientFields;
+ // end dg Part
+
dataCacheDouble &get(const std::string &functionName, dataCache *caller=0);
dataCacheElement &getElement(dataCache *caller=0);
- dataCacheDouble &provideData(std::string name);
- dataCacheMap(int nbEvaluationPoints):_nbEvaluationPoints(nbEvaluationPoints){
- _cacheElement= new dataCacheElement(this);
-}
+ dataCacheDouble &provideData(std::string name, int nRowByPoints, int nCol);
+ dataCacheMap(){
+ _cacheElement = new dataCacheElement(this);
+ _nbEvaluationPoints = 0;
+ }
+ void setNbEvaluationPoints(int nbEvaluationPoints);
inline int getNbEvaluationPoints(){return _nbEvaluationPoints;}
~dataCacheMap();
diff --git a/Solver/luaFunction.cpp b/Solver/luaFunction.cpp
index ba1225216a..f53b0a979c 100644
--- a/Solver/luaFunction.cpp
+++ b/Solver/luaFunction.cpp
@@ -14,7 +14,7 @@ class functionLua::data : public dataCacheDouble{
std::vector<dataCacheDouble *> _dependencies;
public:
data(const functionLua *f, dataCacheMap *m):
- dataCacheDouble(*m,m->getNbEvaluationPoints(),f->_nbCol),
+ dataCacheDouble(*m,1,f->_nbCol),
_function(f)
{
_dependencies.resize ( _function->_dependenciesName.size());
--
GitLab