From 4df00803cdf5e7ee5e5528159014b3d352025e70 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Mon, 15 Mar 2010 14:06:45 +0000
Subject: [PATCH] no more getName() on functions, use functionCoordinates.get()
when you need XYZ
---
Common/LuaBindings.cpp | 1 -
Solver/TESTCASES/Advection1D.lua | 17 +-
Solver/TESTCASES/Advection3D.lua | 19 +--
Solver/TESTCASES/AdvectionDiffusion.lua | 30 +---
Solver/TESTCASES/ConvergenceMultirate.lua | 15 +-
Solver/TESTCASES/Stommel.lua | 13 +-
Solver/TESTCASES/square.geo | 9 +-
Solver/TESTCASES/validation/Advection1D.lua | 14 +-
Solver/TESTCASES/validation/WavePulse.lua | 6 +-
Solver/dgConservationLaw.cpp | 79 ++++-----
Solver/dgConservationLaw.h | 5 +-
Solver/dgConservationLawAdvection.cpp | 71 +++++----
Solver/dgConservationLawAdvection.h | 7 +-
Solver/dgConservationLawMaxwell.cpp | 45 +++---
Solver/dgConservationLawMaxwell.h | 4 +-
Solver/dgConservationLawPerfectGas.cpp | 109 ++++++++-----
Solver/dgConservationLawPerfectGas.h | 36 -----
Solver/dgConservationLawShallowWater1d.cpp | 22 +--
Solver/dgConservationLawShallowWater2d.cpp | 38 ++---
Solver/dgConservationLawWaveEquation.cpp | 14 +-
Solver/dgDofContainer.cpp | 4 +-
Solver/dgDofContainer.h | 3 +-
Solver/dgFunctionIntegrator.cpp | 11 +-
Solver/dgFunctionIntegrator.h | 5 +-
Solver/dgSystemOfEquations.cpp | 4 +-
Solver/dgSystemOfEquations.h | 2 +-
Solver/function.cpp | 167 ++++++++++++++------
Solver/function.h | 32 ++--
Solver/luaFunction.cpp | 12 +-
Solver/luaFunction.h | 4 +-
30 files changed, 410 insertions(+), 388 deletions(-)
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 3b12907cd2..d1523dc3de 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -365,7 +365,6 @@ binding::binding(){
dgFunctionIntegrator::registerBindings(this);
fullMatrix<double>::registerBindings(this);
function::registerBindings(this);
- function::registerDefaultFunctions();
functionLua::registerBindings(this);
gmshOptions::registerBindings(this);
Msg::registerBindings(this);
diff --git a/Solver/TESTCASES/Advection1D.lua b/Solver/TESTCASES/Advection1D.lua
index c6f26566ef..f4a65368f1 100644
--- a/Solver/TESTCASES/Advection1D.lua
+++ b/Solver/TESTCASES/Advection1D.lua
@@ -26,21 +26,10 @@ order=1
dimension=1
-- conservation law
--- advection speed
-v=fullMatrix(3,1);
-v:set(0,0,0.15)
-v:set(1,0,0)
-v:set(2,0,0)
--- diffusivity
-nu=fullMatrix(1,1);
-nu:set(0,0,0)
-
-law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(), '')
+law = dgConservationLawAdvectionDiffusion.advectionLaw(functionConstant({0.15,0,0}))
-- boundary condition
-outside=fullMatrix(1,1)
-outside:set(0,0,0.)
-bndcondition=law:newOutsideValueBoundary(functionConstant(outside):getName())
+bndcondition=law:newOutsideValueBoundary(functionConstant({0}))
law:addBoundaryCondition('Left',bndcondition)
law:addBoundaryCondition('Right',bndcondition)
@@ -53,7 +42,7 @@ limiter = dgSlopeLimiter(law)
rk:setLimiter(limiter)
-- build solution vector
-FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
+FS = functionLua(1, 'initial_condition', {functionCoordinates.get()})
solution = dgDofContainer(groups, law:getNbFields())
solution:L2Projection(FS)
diff --git a/Solver/TESTCASES/Advection3D.lua b/Solver/TESTCASES/Advection3D.lua
index 9c9ab28733..2b5eeccf38 100644
--- a/Solver/TESTCASES/Advection3D.lua
+++ b/Solver/TESTCASES/Advection3D.lua
@@ -6,7 +6,7 @@ dg = dgSystemOfEquations(model)
dg:setOrder(1)
-- initial condition
-function initial_condition( xyz , f )
+function initial_condition(xyz ,f)
for i=0,xyz:size1()-1 do
x = xyz:get(i,0)
y = xyz:get(i,1)
@@ -20,19 +20,14 @@ end
-- conservation law
-- advection speed
-v=fullMatrix(3,1);
-v:set(0,0,0)
-v:set(1,0,0)
-v:set(2,0,0.15)
-
-law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(),'')
+law = dgConservationLawAdvectionDiffusion.advectionLaw(functionConstant({0,0,0.15}))
dg:setConservationLaw(law)
-- boundary condition
-outside=fullMatrix(1,1)
-outside:set(0,0,0.15)
---freestreem=law:newOutsideValueBoundary(functionConstant(outside):getName())
-freestream=law:newOutsideValueBoundary(functionLua(1,'initial_condition',{'XYZ'}):getName())
+--freestreem=law:newOutsideValueBoundary(functionConstant({0,0,0.15}))
+
+xyz = functionCoordinates.get()
+freestream=law:newOutsideValueBoundary(functionLua(1,'initial_condition',{functionCoordinates:get()}))
law:addBoundaryCondition('wall',law:new0FluxBoundary())
law:addBoundaryCondition('inlet',freestream)
@@ -41,7 +36,7 @@ law:addBoundaryCondition('symmetry',law:newSymmetryBoundary())
dg:setup()
-dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
+dg:L2Projection(functionLua(1,'initial_condition',{xyz}))
dg:exportSolution('output/Adv3D-00000')
print'***exporting init solution ***'
diff --git a/Solver/TESTCASES/AdvectionDiffusion.lua b/Solver/TESTCASES/AdvectionDiffusion.lua
index 13b691ddc7..a27a315052 100644
--- a/Solver/TESTCASES/AdvectionDiffusion.lua
+++ b/Solver/TESTCASES/AdvectionDiffusion.lua
@@ -2,35 +2,15 @@ model = GModel ()
model:load ('square.geo')
model:load ('square.msh')
-model2 = GModel()
-model2:load('')
-
-vmodel = {GModel(),GModel()}
-vmodel[1]:load('')
-vmodel[2]:load('')
-
dg = dgSystemOfEquations (model)
-dg:setOrder(3)
+dg:setOrder(4)
-- conservation law
-
--- advection speed
-v=fullMatrix(3,1);
-v:set(0,0,0.15)
-v:set(1,0,0.05)
-v:set(2,0,0)
-
--- diffusivity
-nu=fullMatrix(1,1);
-nu:set(0,0,0.001)
-
-law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(),functionConstant(nu):getName())
+law = dgConservationLawAdvectionDiffusion(functionConstant({0.15,0.05,0}),functionConstant({0.001}))
dg:setConservationLaw(law)
-- boundary condition
-outside=fullMatrix(1,1)
-outside:set(0,0,0.)
-law:addBoundaryCondition('Border',law:newOutsideValueBoundary(functionConstant(outside):getName()))
+law:addBoundaryCondition('Border',law:newOutsideValueBoundary(functionConstant({0,0,0})))
dg:setup()
@@ -43,13 +23,13 @@ function initial_condition( xyz , f )
f:set (i, 0, math.exp(-100*((x-0.2)^2 +(y-0.3)^2)))
end
end
-dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
+dg:L2Projection(functionLua(1,'initial_condition',{functionCoordinates.get()}))
dg:exportSolution('output/Advection_00000')
-- main loop
for i=1,10000 do
- norm = dg:RK44(0.01)
+ norm = dg:RK44(0.001)
if (i % 10 == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/Advection-%05d", i))
diff --git a/Solver/TESTCASES/ConvergenceMultirate.lua b/Solver/TESTCASES/ConvergenceMultirate.lua
index 71ac0ad30d..ad094513a8 100644
--- a/Solver/TESTCASES/ConvergenceMultirate.lua
+++ b/Solver/TESTCASES/ConvergenceMultirate.lua
@@ -50,15 +50,16 @@ elseif (order == 5) then
end
print'*** Create a dg solver ***'
-velF = functionLua(3, 'velocity', {'XYZ'}):getName()
-law=dgConservationLawAdvectionDiffusion(velF,"")
+xyzF = functionCoordinates.get()
+velF = functionLua(3, 'velocity', {xyzF})
+nuF = functionConstant({0})
+law = dgConservationLawAdvectionDiffusion(velF, nuF)
--- law:addBoundaryCondition('Walls',law:newOutsideValueBoundary(functionConstant(g):getName()))
--- law:addBoundaryCondition('Top',law:newOutsideValueBoundary(functionConstant(g):getName()))
+-- law:addBoundaryCondition('Walls',law:newOutsideValueBoundary(functionConstant(g)))
+-- law:addBoundaryCondition('Top',law:newOutsideValueBoundary(functionConstant(g)))
law:addBoundaryCondition('Walls',law:new0FluxBoundary())
law:addBoundaryCondition('Top',law:new0FluxBoundary())
-FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
-
+FS = functionLua(1, 'initial_condition', {xyzF})
GC=dgGroupCollection(myModel,2,order)
solTmp=dgDofContainer(GC,1)
solTmp:L2Projection(FS)
@@ -84,7 +85,7 @@ end
--multirateRK=dgRungeKuttaMultirate43(GC,law)
multirateRK=dgRungeKuttaMultirate22(GC,law)
-integrator=dgFunctionIntegrator(functionLua(1, 'L2Norm', {'Solution'}):getName())
+integrator=dgFunctionIntegrator(functionLua(1, 'L2Norm', {functionSolution.get()}))
oldnorm=1
for n=0,3 do
diff --git a/Solver/TESTCASES/Stommel.lua b/Solver/TESTCASES/Stommel.lua
index 04fba32a3a..26be9b8cf9 100644
--- a/Solver/TESTCASES/Stommel.lua
+++ b/Solver/TESTCASES/Stommel.lua
@@ -29,12 +29,12 @@ Msg.barrier()
claw = dgConservationLawShallowWater2d()
claw:addBoundaryCondition('Wall',claw:newBoundaryWall())
-zero = functionConstant({0}):getName();
-claw:setCoriolisFactor(functionC("tmp.dylib","coriolis",1,{"XYZ"}):getName())
-claw:setQuadraticDissipation(zero)
-claw:setLinearDissipation(functionConstant({1e-6}):getName())
-claw:setSource(functionC("tmp.dylib","wind",2,{"XYZ"}):getName())
-claw:setBathymetry(functionConstant({1000}):getName())
+XYZ = functionCoordinates.get();
+claw:setCoriolisFactor(functionC("tmp.dylib","coriolis",1,{XYZ}))
+claw:setQuadraticDissipation(functionConstant({0}))
+claw:setLinearDissipation(functionConstant({1e-6}))
+claw:setSource(functionC("tmp.dylib","wind",2,{XYZ}))
+claw:setBathymetry(functionConstant({1000}))
groups = dgGroupCollection(model, dimension, order)
groups:buildGroupsOfInterfaces()
solution = dgDofContainer(groups, claw:getNbFields())
@@ -50,3 +50,4 @@ for i=1,60000 do
solution:exportMsh(string.format('output/solution-%06d',i))
end
end
+
diff --git a/Solver/TESTCASES/square.geo b/Solver/TESTCASES/square.geo
index d2bd216e0f..d9302dde17 100644
--- a/Solver/TESTCASES/square.geo
+++ b/Solver/TESTCASES/square.geo
@@ -1,7 +1,8 @@
-Point(1) = {0.0, 0.0, 0, .03};
-Point(2) = {1, 0.0, 0, .03};
-Point(3) = {1, 1, 0, .03};
-Point(4) = {0, 1, 0, .03};
+lc = 0.1;
+Point(1) = {0.0, 0.0, 0, lc};
+Point(2) = {1, 0.0, 0, lc};
+Point(3) = {1, 1, 0, lc};
+Point(4) = {0, 1, 0, lc};
Line(1) = {4, 3};
Line(2) = {3, 2};
Line(3) = {2, 1};
diff --git a/Solver/TESTCASES/validation/Advection1D.lua b/Solver/TESTCASES/validation/Advection1D.lua
index c6f26566ef..a15fdf9616 100644
--- a/Solver/TESTCASES/validation/Advection1D.lua
+++ b/Solver/TESTCASES/validation/Advection1D.lua
@@ -27,20 +27,12 @@ dimension=1
-- conservation law
-- advection speed
-v=fullMatrix(3,1);
-v:set(0,0,0.15)
-v:set(1,0,0)
-v:set(2,0,0)
--- diffusivity
-nu=fullMatrix(1,1);
-nu:set(0,0,0)
-
-law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(), '')
+law = dgConservationLawAdvectionDiffusion:advectionLaw(functionConstant({0.15,0,0}))
-- boundary condition
outside=fullMatrix(1,1)
outside:set(0,0,0.)
-bndcondition=law:newOutsideValueBoundary(functionConstant(outside):getName())
+bndcondition=law:newOutsideValueBoundary(functionConstant({0,0,0}))
law:addBoundaryCondition('Left',bndcondition)
law:addBoundaryCondition('Right',bndcondition)
@@ -53,7 +45,7 @@ limiter = dgSlopeLimiter(law)
rk:setLimiter(limiter)
-- build solution vector
-FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
+FS = functionLua(1, 'initial_condition', {functionCoordinates.get()})
solution = dgDofContainer(groups, law:getNbFields())
solution:L2Projection(FS)
diff --git a/Solver/TESTCASES/validation/WavePulse.lua b/Solver/TESTCASES/validation/WavePulse.lua
index abf15298e4..16cf5b4107 100644
--- a/Solver/TESTCASES/validation/WavePulse.lua
+++ b/Solver/TESTCASES/validation/WavePulse.lua
@@ -22,8 +22,10 @@ myModel = GModel ()
--myModel:load('box.geo')
--myModel:load('box.msh')
--myModel:load('square_quads.msh')
-myModel:load('square_part.msh')
+-- myModel:load('square_part.msh')
+myModel:load('square.geo')
myModel:mesh(2)
+myModel:save('square.msh')
print'*** Create a dg solver ***'
DG = dgSystemOfEquations (myModel)
@@ -33,7 +35,7 @@ DG:setConservationLaw(law)
law:addBoundaryCondition('Border',law:newBoundaryWall())
DG:setup()
-initialCondition = functionLua(3,'initial_condition',{'XYZ'}):getName()
+initialCondition = functionLua(3,'initial_condition',{functionCoordinates.get()})
print'*** setting the initial solution ***'
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index de184a8a2a..2b66374db7 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/Solver/dgConservationLaw.cpp
@@ -1,7 +1,7 @@
#include "dgConservationLaw.h"
#include "function.h"
class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
- std::string _outsideValueFunctionName;
+ const function *_outsideValueFunction;
class term : public dataCacheDouble {
dataCacheMap cacheMapRight; // new cacheMap to pass to the Riemann solver
dataCacheDouble &solutionRight;
@@ -9,10 +9,10 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dataCacheDouble *riemannSolver;
dgConservationLaw *_claw;
public:
- term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
+ term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const function *outsideValueFunction):
dataCacheDouble(cacheMapLeft, 1, claw->getNbFields()),
solutionRight(cacheMapRight.provideSolution(claw->getNbFields())),
- outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
+ outsideValue(cacheMapLeft.get(outsideValueFunction,this)),
_claw(claw)
{
cacheMapRight.setNbEvaluationPoints(cacheMapLeft.getNbEvaluationPoints());
@@ -33,9 +33,9 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
class dirichlet : public dataCacheDouble {
dataCacheDouble &outsideValue;
public:
- dirichlet(dgConservationLaw *claw, dataCacheMap &cacheMap,const std::string outsideValueFunctionName):
+ dirichlet(dgConservationLaw *claw, dataCacheMap &cacheMap,const function *outsideValueFunction):
dataCacheDouble(cacheMap, 1, claw->getNbFields()),
- outsideValue(cacheMap.get(outsideValueFunctionName,this)){}
+ outsideValue(cacheMap.get(outsideValueFunction,this)){}
void _eval () {
for(int i=0;i<_value.size1();i++)
for(int j=0;j<_value.size2();j++)
@@ -49,10 +49,10 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dataCacheDouble *maxDif;
dgConservationLaw *_claw;
public:
- maximumDiffusivity(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
+ maximumDiffusivity(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const function *outsideValueFunction):
dataCacheDouble(cacheMapLeft, 1, claw->getNbFields()),
solutionRight(cacheMapRight.provideSolution(claw->getNbFields())),
- outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
+ outsideValue(cacheMapLeft.get(outsideValueFunction,this)),
_claw(claw)
{
cacheMapRight.setNbEvaluationPoints(cacheMapLeft.getNbEvaluationPoints());
@@ -70,27 +70,27 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
}
};
public:
- dgBoundaryConditionOutsideValue(dgConservationLaw *claw,const std::string outsideValueFunctionName): dgBoundaryCondition(claw),
- _outsideValueFunctionName(outsideValueFunctionName)
+ dgBoundaryConditionOutsideValue(dgConservationLaw *claw,const function *outsideValueFunction): dgBoundaryCondition(claw),
+ _outsideValueFunction(outsideValueFunction)
{ }
dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
- return new term(_claw,cacheMapLeft,_outsideValueFunctionName);
+ return new term(_claw,cacheMapLeft,_outsideValueFunction);
}
dataCacheDouble *newDiffusiveDirichletBC(dataCacheMap &cacheMapLeft) const {
- return new dirichlet(_claw,cacheMapLeft,_outsideValueFunctionName);
+ return new dirichlet(_claw,cacheMapLeft,_outsideValueFunction);
}
dataCacheDouble *newMaximumDiffusivity(dataCacheMap &cacheMapLeft) const {
- return new maximumDiffusivity(_claw,cacheMapLeft,_outsideValueFunctionName);
+ return new maximumDiffusivity(_claw,cacheMapLeft,_outsideValueFunction);
}
};
class dgBoundaryConditionNeumann : public dgBoundaryCondition {
- std::string _fluxFunctionName;
+ const function *_fluxFunction;
class term : public dataCacheDouble {
dataCacheDouble &flux;
public:
- term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string fluxFunctionName):
+ term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const function *fluxFunction):
dataCacheDouble(cacheMapLeft,1 ,claw->getNbFields()),
- flux(cacheMapLeft.get(fluxFunctionName,this))
+ flux(cacheMapLeft.get(fluxFunction,this))
{}
void _eval() {
for(int i=0;i<_value.size1(); i++)
@@ -99,11 +99,11 @@ class dgBoundaryConditionNeumann : public dgBoundaryCondition {
}
};
public:
- dgBoundaryConditionNeumann(dgConservationLaw *claw,const std::string fluxFunctionName): dgBoundaryCondition(claw),
- _fluxFunctionName(fluxFunctionName)
+ dgBoundaryConditionNeumann(dgConservationLaw *claw,const function *fluxFunction): dgBoundaryCondition(claw),
+ _fluxFunction(fluxFunction)
{ }
dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
- return new term(_claw,cacheMapLeft,_fluxFunctionName);
+ return new term(_claw,cacheMapLeft,_fluxFunction);
}
};
@@ -152,51 +152,36 @@ public:
dgBoundaryCondition *dgConservationLaw::newSymmetryBoundary() {
return new dgBoundarySymmetry(this);
}
-dgBoundaryCondition *dgConservationLaw::newOutsideValueBoundary(const std::string outsideValueFunctionName) {
- return new dgBoundaryConditionOutsideValue(this,outsideValueFunctionName);
+dgBoundaryCondition *dgConservationLaw::newOutsideValueBoundary(const function *outsideValueFunction) {
+ return new dgBoundaryConditionOutsideValue(this,outsideValueFunction);
}
-dgBoundaryCondition *dgConservationLaw::newNeumannBoundary(const std::string fluxFunctionName) {
- return new dgBoundaryConditionNeumann(this,fluxFunctionName);
+dgBoundaryCondition *dgConservationLaw::newNeumannBoundary(const function *fluxFunction) {
+ return new dgBoundaryConditionNeumann(this,fluxFunction);
}
dgBoundaryCondition *dgConservationLaw::new0FluxBoundary() {
return new dgBoundaryCondition0Flux(this);
}
-
-class dgBoundaryCondition::dirichlet_ : public dataCacheDouble {
- dataCacheDouble /
-public:
- dirichlet_(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap,1,cacheMap.get("Solution")().size2()),
- sol(cacheMap.get("Solution",this)){}
- void _eval () {
- for(int i=0; i< _value.size1(); i++)
- for (int k=0;k<sol().size2();k++)
- _value(i,k) = sol(i,k);
- }
-};
-
class dgBoundaryCondition::neumann_ : public dataCacheDouble {
- dgConservationLaw *_claw;
- dataCacheDouble &sol,&normals;
+ dataCacheDouble &normals;
dataCacheDouble *diffusiveFlux;
public:
neumann_(dataCacheMap &cacheMap, dgConservationLaw *claw):
- dataCacheDouble(cacheMap,1,cacheMap.get("Solution")().size2()),
- _claw (claw),
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)){
- diffusiveFlux=_claw->newDiffusiveFlux(cacheMap);
+ dataCacheDouble(cacheMap,1,claw->getNbFields()),
+ normals(cacheMap.getNormals(this))
+ {
+ diffusiveFlux=claw->newDiffusiveFlux(cacheMap);
if (diffusiveFlux)diffusiveFlux->addMeAsDependencyOf(this);
}
void _eval () {
+ int nbFields = _value.size2();
const fullMatrix<double> &dfl = (*diffusiveFlux)();
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);
+ dfl(i,k+nbFields*0) *(normals)(0,i) +
+ dfl(i,k+nbFields*1) *(normals)(1,i) +
+ dfl(i,k+nbFields*2) *(normals)(2,i);
}
}
}
@@ -204,7 +189,7 @@ public:
};
dataCacheDouble *dgBoundaryCondition::newDiffusiveDirichletBC(dataCacheMap &cacheMapLeft) const {
- return new dirichlet_(cacheMapLeft);
+ return &cacheMapLeft.getSolution(NULL);
}
dataCacheDouble *dgBoundaryCondition::newDiffusiveNeumannBC(dataCacheMap &cacheMapLeft) const {
return new dgBoundaryCondition::neumann_(cacheMapLeft, _claw);
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 37c46b7d93..2e8c790735 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -9,6 +9,7 @@
#include "fullMatrix.h"
class dataCacheDouble;
class dataCacheMap;
+class function;
class dgConservationLaw;
class binding;
@@ -62,8 +63,8 @@ class dgConservationLaw {
}
//a generic boundary condition using the Riemann solver of the conservation Law
- dgBoundaryCondition *newOutsideValueBoundary(std::string outsideValueFunctionName);
- dgBoundaryCondition *newNeumannBoundary(const std::string fluxFunctionName);
+ dgBoundaryCondition *newOutsideValueBoundary(const function *outsideValueFunction);
+ dgBoundaryCondition *newNeumannBoundary(const function *fluxFunction);
dgBoundaryCondition *new0FluxBoundary();
dgBoundaryCondition *newSymmetryBoundary();
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index e9bc046854..2b32ff8bdc 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -9,10 +9,10 @@
class dgConservationLawAdvectionDiffusion::advection : public dataCacheDouble {
dataCacheDouble &sol, &v;
public:
- advection(std::string vFunctionName, dataCacheMap &cacheMap):
+ advection(const function *vFunction, dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,3),
- sol(cacheMap.get("Solution",this)),
- v(cacheMap.get(vFunctionName,this))
+ sol(cacheMap.getSolution(this)),
+ v(cacheMap.get(vFunction,this))
{};
void _eval () {
for(int i=0; i< sol().size1(); i++) {
@@ -26,9 +26,9 @@ class dgConservationLawAdvectionDiffusion::advection : public dataCacheDouble {
class dgConservationLawAdvectionDiffusion::maxConvectiveSpeed : public dataCacheDouble {
dataCacheDouble &v;
public:
- maxConvectiveSpeed(std::string vFunctionName,dataCacheMap &cacheMap):
+ maxConvectiveSpeed(const function *vFunction,dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,1),
- v(cacheMap.get(vFunctionName,this))
+ v(cacheMap.get(vFunction,this))
{
};
void _eval () {
@@ -42,12 +42,12 @@ class dgConservationLawAdvectionDiffusion::maxConvectiveSpeed : public dataCache
class dgConservationLawAdvectionDiffusion::riemann : public dataCacheDouble {
dataCacheDouble &normals, &solLeft, &solRight,&v;
public:
- riemann(std::string vFunctionName, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+ riemann(const function *vFunction, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
dataCacheDouble(cacheMapLeft,1,2),
- normals(cacheMapLeft.get("Normals", this)),
- solLeft(cacheMapLeft.get("Solution", this)),
- solRight(cacheMapRight.get("Solution", this)),
- v(cacheMapLeft.get(vFunctionName,this))
+ normals(cacheMapLeft.getNormals(this)),
+ solLeft(cacheMapLeft.getSolution(this)),
+ solRight(cacheMapRight.getSolution(this)),
+ v(cacheMapLeft.get(vFunction,this))
{};
void _eval () {
for(int i=0; i< _value.size1(); i++) {
@@ -65,10 +65,10 @@ class dgConservationLawAdvectionDiffusion::riemann : public dataCacheDouble {
class dgConservationLawAdvectionDiffusion::diffusion : public dataCacheDouble {
dataCacheDouble &solgrad, ν
public:
- diffusion(std::string nuFunctionName, dataCacheMap &cacheMap):
+ diffusion(const function * nuFunction, dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,3),
- solgrad(cacheMap.get("SolutionGradient",this)),
- nu(cacheMap.get(nuFunctionName,this))
+ solgrad(cacheMap.getSolutionGradient(this)),
+ nu(cacheMap.get(nuFunction,this))
{};
void _eval () {
for(int i=0; i< solgrad().size1()/3; i++) {
@@ -79,45 +79,48 @@ class dgConservationLawAdvectionDiffusion::diffusion : public dataCacheDouble {
}
};
dataCacheDouble *dgConservationLawAdvectionDiffusion::newConvectiveFlux( dataCacheMap &cacheMap) const {
- if( !_vFunctionName.empty())
- return new advection(_vFunctionName,cacheMap);
- else
- return NULL;
+ return _vFunction ? new advection(_vFunction,cacheMap) : NULL;
}
dataCacheDouble *dgConservationLawAdvectionDiffusion::newMaxConvectiveSpeed( dataCacheMap &cacheMap) const {
- return new maxConvectiveSpeed(_vFunctionName,cacheMap);
+ return _vFunction ? new maxConvectiveSpeed(_vFunction,cacheMap) : NULL ;
+
}
dataCacheDouble *dgConservationLawAdvectionDiffusion::newMaximumDiffusivity( dataCacheMap &cacheMap) const {
- if( !_nuFunctionName.empty())
- return &cacheMap.get(_nuFunctionName);
- else
- return NULL;
+ return _nuFunction ? &cacheMap.get(_nuFunction) : NULL;
}
dataCacheDouble *dgConservationLawAdvectionDiffusion::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
- if( !_vFunctionName.empty())
- return new riemann(_vFunctionName,cacheMapLeft, cacheMapRight);
- else
- return NULL;
+ return _vFunction ? new riemann(_vFunction,cacheMapLeft, cacheMapRight) : NULL;
}
dataCacheDouble *dgConservationLawAdvectionDiffusion::newDiffusiveFlux( dataCacheMap &cacheMap) const {
- if( !_nuFunctionName.empty())
- return new diffusion(_nuFunctionName,cacheMap);
- else
- return NULL;
+ return _nuFunction ? new diffusion(_nuFunction,cacheMap) : NULL;
+}
+
+dgConservationLawAdvectionDiffusion *dgConservationLawAdvectionDiffusion::diffusionLaw(const function *nu) {
+ return new dgConservationLawAdvectionDiffusion(NULL, nu);
}
-dgConservationLawAdvectionDiffusion::dgConservationLawAdvectionDiffusion(std::string vFunctionName, std::string nuFunctionName)
+dgConservationLawAdvectionDiffusion *dgConservationLawAdvectionDiffusion::advectionLaw(const function *v) {
+ return new dgConservationLawAdvectionDiffusion(v,NULL);
+}
+
+dgConservationLawAdvectionDiffusion::dgConservationLawAdvectionDiffusion(const function *vFunction, const function *nuFunction)
{
- _vFunctionName = vFunctionName;
- _nuFunctionName = nuFunctionName;
+ _vFunction = vFunction;
+ _nuFunction = nuFunction;
_nbf = 1;
}
#include "Bindings.h"
void dgConservationLawAdvectionDiffusionRegisterBindings (binding *b){
classBinding *cb = b->addClass<dgConservationLawAdvectionDiffusion>("dgConservationLawAdvectionDiffusion");
methodBinding *cm;
- cm = cb->setConstructor<dgConservationLawAdvectionDiffusion,std::string,std::string>();
+ cm = cb->setConstructor<dgConservationLawAdvectionDiffusion,const function *, const function*>();
cm->setArgNames("v","nu",NULL);
cm->setDescription("A new advection-diffusion law. The advection speed is given by 'v' vector function and the (scalar) diffusivity is given by the function 'nu'.");
+ cm = cb->addMethod("diffusionLaw",dgConservationLawAdvectionDiffusion::diffusionLaw);
+ cm->setDescription("static method to create an instance with diffusion only.");
+ cm->setArgNames("nu",NULL);
+ cm = cb->addMethod("advectionLaw",dgConservationLawAdvectionDiffusion::advectionLaw);
+ cm->setDescription("static method to create an instance with advection only.");
+ cm->setArgNames("v",NULL);
cb->setParentClass<dgConservationLaw>();
cb->setDescription("Advection and diffusion of a scalar, the advection and diffusion are provided by functions.");
}
diff --git a/Solver/dgConservationLawAdvection.h b/Solver/dgConservationLawAdvection.h
index 65316de058..758acf140b 100644
--- a/Solver/dgConservationLawAdvection.h
+++ b/Solver/dgConservationLawAdvection.h
@@ -3,8 +3,9 @@
#include "dgConservationLaw.h"
//AdvectionDiffusion diffusion equation
//dc/dt + v div(c) - nu lapl(c) = 0
+class function;
class dgConservationLawAdvectionDiffusion : public dgConservationLaw {
- std::string _vFunctionName,_nuFunctionName;
+ const function *_vFunction, *_nuFunction;
class advection;
class riemann;
class diffusion;
@@ -15,7 +16,9 @@ class dgConservationLawAdvectionDiffusion : public dgConservationLaw {
dataCacheDouble *newMaxConvectiveSpeed( dataCacheMap &cacheMap) const;
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
- dgConservationLawAdvectionDiffusion(std::string vFunctionName, std::string nuFunctionName);
+ dgConservationLawAdvectionDiffusion(const function *vFunction, const function *nuFunction);
+ static dgConservationLawAdvectionDiffusion *diffusionLaw(const function *nuFunction);
+ static dgConservationLawAdvectionDiffusion *advectionLaw(const function *vFunction);
};
class binding;
void dgConservationLawAdvectionDiffusionRegisterBindings (binding *b);
diff --git a/Solver/dgConservationLawMaxwell.cpp b/Solver/dgConservationLawMaxwell.cpp
index 9ca9c3846e..49743b6da8 100644
--- a/Solver/dgConservationLawMaxwell.cpp
+++ b/Solver/dgConservationLawMaxwell.cpp
@@ -6,10 +6,10 @@
class dgConservationLawMaxwell::advection : public dataCacheDouble {
dataCacheDouble &sol, &mu_eps;
public:
- advection(std::string mu_epsFunctionName,dataCacheMap &cacheMap):
+ advection(const function *mu_epsFunction,dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,18),
- sol(cacheMap.get("Solution",this)),
- mu_eps(cacheMap.get(mu_epsFunctionName,this))
+ sol(cacheMap.getSolution(this)),
+ mu_eps(cacheMap.get(mu_epsFunction,this))
{};
void _eval () {
int nQP = sol().size1();
@@ -49,15 +49,14 @@ class dgConservationLawMaxwell::advection : public dataCacheDouble {
class dgConservationLawMaxwell::source: public dataCacheDouble {
- dataCacheDouble &xyz, &solution;
+ dataCacheDouble &solution;
public :
source(dataCacheMap &cacheMap) :
dataCacheDouble(cacheMap,1,6),
- xyz(cacheMap.get("XYZ",this)),
- solution(cacheMap.get("Solution",this))
+ solution(cacheMap.getSolution(this))
{}
void _eval () {
- int nQP = xyz().size1();
+ int nQP = _value.size1();
for (int i=0; i<nQP; i++) {
_value (i,0) = 0;
@@ -74,12 +73,12 @@ class dgConservationLawMaxwell::source: public dataCacheDouble {
class dgConservationLawMaxwell::riemann:public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR, &mu_eps;
public:
- riemann(std::string mu_epsFunctionName,dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+ riemann(const function *mu_epsFunction,dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
dataCacheDouble(cacheMapLeft,1,12),
- normals(cacheMapLeft.get("Normals", this)),
- solL(cacheMapLeft.get("Solution", this)),
- solR(cacheMapRight.get("Solution", this)),
- mu_eps(cacheMapLeft.get(mu_epsFunctionName,this))
+ normals(cacheMapLeft.getNormals( this)),
+ solL(cacheMapLeft.getSolution( this)),
+ solR(cacheMapRight.getSolution( this)),
+ mu_eps(cacheMapLeft.get(mu_epsFunction,this))
{};
void _eval () {
int nQP = solL().size1();
@@ -154,18 +153,12 @@ class dgConservationLawMaxwell::riemann:public dataCacheDouble {
dataCacheDouble *dgConservationLawMaxwell::newConvectiveFlux( dataCacheMap &cacheMap) const {
- if( !_mu_epsFunctionName.empty())
- return new advection(_mu_epsFunctionName,cacheMap);
- else
- return NULL;
+ return _mu_epsFunction ? new advection(_mu_epsFunction,cacheMap) : NULL;
}
dataCacheDouble *dgConservationLawMaxwell::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
- if( !_mu_epsFunctionName.empty())
- return new riemann(_mu_epsFunctionName,cacheMapLeft, cacheMapRight);
- else
- return NULL;
+ return _mu_epsFunction ? new riemann(_mu_epsFunction,cacheMapLeft, cacheMapRight) : NULL;
}
@@ -173,9 +166,9 @@ dataCacheDouble *dgConservationLawMaxwell::newSourceTerm (dataCacheMap &cacheMap
return new source(cacheMap);
}
-dgConservationLawMaxwell::dgConservationLawMaxwell(std::string mu_epsFunctionName)
+dgConservationLawMaxwell::dgConservationLawMaxwell(const function *mu_epsFunction)
{
- _mu_epsFunctionName = mu_epsFunctionName;
+ _mu_epsFunction = mu_epsFunction;
_nbf = 6;
}
@@ -186,8 +179,8 @@ class dgBoundaryConditionMaxwellWall : public dgBoundaryCondition {
public:
term(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,6),
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this))
+ sol(cacheMap.getSolution(this)),
+ normals(cacheMap.getNormals(this))
{}
void _eval () {
int nQP = sol().size1();
@@ -235,9 +228,9 @@ void dgConservationLawMaxwellRegisterBindings (binding *b){
methodBinding *cm;
cm = cb->addMethod("newBoundaryWall",&dgConservationLawMaxwell::newBoundaryWall);
cm->setDescription("wall boundary");
- cm = cb->setConstructor<dgConservationLawMaxwell, std::string>();
+ cm = cb->setConstructor<dgConservationLawMaxwell, const function *>();
cm->setArgNames("mu_eps",NULL);
cm->setDescription("A new Maxwell conservation law.");
cb->setParentClass<dgConservationLaw>();
- cb->setDescription("Advection for Maxwell equation.");
+ cb->setDescription("Advection for Maxwell equation.");
}
diff --git a/Solver/dgConservationLawMaxwell.h b/Solver/dgConservationLawMaxwell.h
index dd8457e381..aa4e15ac64 100644
--- a/Solver/dgConservationLawMaxwell.h
+++ b/Solver/dgConservationLawMaxwell.h
@@ -5,7 +5,7 @@
class dgConservationLawMaxwell : public dgConservationLaw {
- std::string _mu_epsFunctionName;
+ const function *_mu_epsFunction;
class advection;
class source;
class riemann;
@@ -14,7 +14,7 @@ class dgConservationLawMaxwell : public dgConservationLaw {
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
dgBoundaryCondition *newBoundaryWall () ;
- dgConservationLawMaxwell(std::string mu_epsFunctionName);
+ dgConservationLawMaxwell(const function *mu_epsFunction);
};
class binding;
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index 7dab969778..834d761c59 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -1,5 +1,41 @@
#include "dgConservationLawPerfectGas.h"
#include "function.h"
+class dgPerfectGasLaw2d : public dgConservationLaw {
+ class advection;
+ class diffusion;
+ class riemann;
+ class riemannGodunov;
+ class source;
+ class maxConvectiveSpeed;
+ class maxDiffusivity;
+ class clipToPhysics;
+ // the name of the functions for
+ // viscosity (_muFunctionName)
+ // thermal conductivity (_kappaFunctionName)
+ const function *_kappaFunction, *_muFunction;
+ const function *_sourceFunction;
+
+ public:
+ dataCacheDouble *newMaxConvectiveSpeed( dataCacheMap &cacheMap) const;
+ dataCacheDouble *newMaximumDiffusivity( dataCacheMap &cacheMap) const;
+ dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const;
+ dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
+ dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
+ dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
+ dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const;
+ dgPerfectGasLaw2d();
+ dgBoundaryCondition *newWallBoundary() ;
+ dgBoundaryCondition *newSlipWallBoundary() ;
+ // sets some coefficients
+ void setViscosityAndThermalConductivity (const function *a, const function *b){
+ _muFunction = a;
+ _kappaFunction = b;
+ }
+ // sets some coefficients
+ void setSource(const function *a){
+ _sourceFunction = a;
+ }
+};
/*
\rho \left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v}\right) = -\nabla p + \mu \nabla^2 \mathbf{v} + \left( \tfrac13 \mu + \mu^v) \nabla (\nabla \cdot \mathbf{v} \right) + \rho \mathbf{f}
@@ -329,7 +365,7 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
public:
advection(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,8),
- sol(cacheMap.get("Solution",this))
+ sol(cacheMap.getSolution(this))
{};
void _eval () {
const int nQP = _value.size1();
@@ -362,13 +398,13 @@ class dgPerfectGasLaw2d::diffusion : public dataCacheDouble {
dataCacheDouble &sol,&grad,&mu,κ
public:
diffusion(dataCacheMap &cacheMap,
- const std::string &muFunctionName,
- const std::string &kappaFunctionName ):
+ const function * muFunction,
+ const function * kappaFunction ):
dataCacheDouble(cacheMap,1,12),
- sol(cacheMap.get("Solution",this)),
- grad(cacheMap.get("SolutionGradient",this)),
- mu(cacheMap.get(muFunctionName,this)),
- kappa(cacheMap.get(kappaFunctionName,this))
+ sol(cacheMap.getSolution(this)),
+ grad(cacheMap.getSolutionGradient(this)),
+ mu(cacheMap.get(muFunction,this)),
+ kappa(cacheMap.get(kappaFunction,this))
{};
void _eval () {
const int nQP = sol().size1();
@@ -421,10 +457,10 @@ class dgPerfectGasLaw2d::diffusion : public dataCacheDouble {
class dgPerfectGasLaw2d::source : public dataCacheDouble {
dataCacheDouble &sol,&s;
public:
- source(dataCacheMap &cacheMap, const std::string &sourceFunctionName):
+ source(dataCacheMap &cacheMap, const function *sourceFunction):
dataCacheDouble(cacheMap,1,4),
- sol(cacheMap.get("Solution",this)),
- s(cacheMap.get(sourceFunctionName,this))
+ sol(cacheMap.getSolution(this)),
+ s(cacheMap.get(sourceFunction,this))
{};
void _eval () {
const int nQP = sol().size1();
@@ -443,7 +479,7 @@ class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
public:
clipToPhysics(dataCacheMap &cacheMap, double presMin, double rhoMin):
dataCacheDouble(cacheMap,1,4),
- sol(cacheMap.get("Solution",this))
+ sol(cacheMap.getSolution(this))
{
_presMin=presMin;
_rhoMin=rhoMin;
@@ -476,9 +512,9 @@ class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
public:
riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
dataCacheDouble(cacheMapLeft,1,12),
- normals(cacheMapLeft.get("Normals", this)),
- solL(cacheMapLeft.get("Solution", this)),
- solR(cacheMapRight.get("Solution", this))
+ normals(cacheMapLeft.getNormals( this)),
+ solL(cacheMapLeft.getSolution( this)),
+ solR(cacheMapRight.getSolution( this))
{};
void _eval () {
int nQP = solL().size1();
@@ -507,9 +543,9 @@ class dgPerfectGasLaw2d::riemannGodunov : public dataCacheDouble {
public:
riemannGodunov(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
dataCacheDouble(cacheMapLeft,1,12),
- normals(cacheMapLeft.get("Normals", this)),
- solL(cacheMapLeft.get("Solution", this)),
- solR(cacheMapRight.get("Solution", this))
+ normals(cacheMapLeft.getNormals( this)),
+ solL(cacheMapLeft.getSolution( this)),
+ solR(cacheMapRight.getSolution( this))
{};
void _eval () {
int nQP = solL().size1();
@@ -571,7 +607,7 @@ class dgPerfectGasLaw2d::maxConvectiveSpeed : public dataCacheDouble {
public:
maxConvectiveSpeed(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,1),
- sol(cacheMap.get("Solution",this))
+ sol(cacheMap.getSolution(this))
{
};
void _eval () {
@@ -589,11 +625,11 @@ class dgPerfectGasLaw2d::maxConvectiveSpeed : public dataCacheDouble {
class dgPerfectGasLaw2d::maxDiffusivity : public dataCacheDouble {
dataCacheDouble &sol,&mu,κ
public:
- maxDiffusivity(dataCacheMap &cacheMap, const std::string &muFunctionName, const std::string &kappaFunctionName ):
+ maxDiffusivity(dataCacheMap &cacheMap, const function *muFunction, const function *kappaFunction ):
dataCacheDouble(cacheMap,1,1),
- sol(cacheMap.get("Solution",this)),
- mu(cacheMap.get(muFunctionName,this)),
- kappa(cacheMap.get(kappaFunctionName,this))
+ sol(cacheMap.getSolution(this)),
+ mu(cacheMap.get(muFunction,this)),
+ kappa(cacheMap.get(kappaFunction,this))
{
};
void _eval () {
@@ -606,7 +642,7 @@ class dgPerfectGasLaw2d::maxDiffusivity : public dataCacheDouble {
dataCacheDouble *dgPerfectGasLaw2d::newMaximumDiffusivity( dataCacheMap &cacheMap) const {
- return _muFunctionName.empty() ? NULL : new maxDiffusivity(cacheMap,_muFunctionName,_kappaFunctionName);
+ return _muFunction ? new maxDiffusivity(cacheMap,_muFunction,_kappaFunction) : NULL;
}
dataCacheDouble *dgPerfectGasLaw2d::newMaxConvectiveSpeed( dataCacheMap &cacheMap) const {
@@ -619,16 +655,13 @@ dataCacheDouble *dgPerfectGasLaw2d::newRiemannSolver( dataCacheMap &cacheMapLeft
return new riemannGodunov(cacheMapLeft, cacheMapRight);
}
dataCacheDouble *dgPerfectGasLaw2d::newDiffusiveFlux( dataCacheMap &cacheMap) const {
- if (_muFunctionName.empty() || _kappaFunctionName.empty())
+ if (!_muFunction || !_kappaFunction)
return 0;
else
- return new diffusion(cacheMap,_muFunctionName,_kappaFunctionName);
+ return new diffusion(cacheMap,_muFunction,_kappaFunction);
}
dataCacheDouble *dgPerfectGasLaw2d::newSourceTerm (dataCacheMap &cacheMap) const {
- if (_sourceFunctionName.empty())
- return 0;
- else
- return new source(cacheMap,_sourceFunctionName);
+ return _sourceFunction ? new source(cacheMap, _sourceFunction) : NULL;
}
dataCacheDouble *dgPerfectGasLaw2d::newClipToPhysics( dataCacheMap &cacheMap) const {
return new clipToPhysics(cacheMap,1e-3,1e-3);
@@ -653,8 +686,8 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
//-------------------------------------------------------------------------------
term(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,4),
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)){}
+ sol(cacheMap.getSolution(this)),
+ normals(cacheMap.getNormals(this)){}
void _eval () {
int nQP = sol().size1();
for(int i=0; i< nQP; i++) {
@@ -698,7 +731,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
public:
dirichletNonSlip(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,4),
- sol(cacheMap.get("Solution",this)){}
+ sol(cacheMap.getSolution(this)){}
void _eval () {
int nQP = sol().size1();
for(int i=0; i< nQP; i++) {
@@ -724,8 +757,8 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
neumannNonSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
dataCacheDouble(cacheMap,1,4),
_claw (claw),
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)){
+ sol(cacheMap.getSolution(this)),
+ normals(cacheMap.getNormals(this)){
diffusiveFlux=_claw->newDiffusiveFlux(cacheMap);
if (diffusiveFlux)diffusiveFlux->addMeAsDependencyOf(this);
}
@@ -755,8 +788,8 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
dataCacheDouble /
public:
dirichletSlip(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap,1,cacheMap.get("Solution")().size2()),
- sol(cacheMap.get("Solution",this)){}
+ dataCacheDouble(cacheMap,1,cacheMap.getSolution(NULL)().size2()),
+ sol(cacheMap.getSolution(this)){}
void _eval () {
int nQP = sol().size1();
for(int i=0; i< nQP; i++) {
@@ -778,8 +811,8 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
neumannSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
_claw (claw),
dataCacheDouble(cacheMap,1,4),
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)){
+ sol(cacheMap.getSolution(this)),
+ normals(cacheMap.getNormals(this)){
}
void _eval () {
int nQP = sol().size1();
diff --git a/Solver/dgConservationLawPerfectGas.h b/Solver/dgConservationLawPerfectGas.h
index 3cbe06c1c4..df0c4904f6 100644
--- a/Solver/dgConservationLawPerfectGas.h
+++ b/Solver/dgConservationLawPerfectGas.h
@@ -6,42 +6,6 @@
// d(rhou)/dt + d(rhou^2+p)/dx + d(rhouv)/dy + d(t_xx)/dx + d(t_xy)/dy= 0
// d(rhov)/dt + d(rhouv)/dx + d(rhov^2+p)/dy + d(t_xy)/dx + d(t_yy)/dy= 0
// d(rhoE)/dt + d(rhouH)/dx + d(rhovH)/dy + d(ut_xx+vt_xy-q_x)/dx + d(ut_xy+vt_yy-q_y)/dy= 0
-class dgPerfectGasLaw2d : public dgConservationLaw {
- class advection;
- class diffusion;
- class riemann;
- class riemannGodunov;
- class source;
- class maxConvectiveSpeed;
- class maxDiffusivity;
- class clipToPhysics;
- // the name of the functions for
- // viscosity (_muFunctionName)
- // thermal conductivity (_kappaFunctionName)
- std::string _kappaFunctionName,_muFunctionName;
- std::string _sourceFunctionName;
-
- public:
- dataCacheDouble *newMaxConvectiveSpeed( dataCacheMap &cacheMap) const;
- dataCacheDouble *newMaximumDiffusivity( dataCacheMap &cacheMap) const;
- dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const;
- dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
- dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
- dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
- dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const;
- dgPerfectGasLaw2d();
- dgBoundaryCondition *newWallBoundary() ;
- dgBoundaryCondition *newSlipWallBoundary() ;
- // sets some coefficients
- void setViscosityAndThermalConductivity (std::string a, std::string b){
- _muFunctionName = a;
- _kappaFunctionName = b;
- }
- // sets some coefficients
- void setSource(std::string a){
- _sourceFunctionName = a;
- }
-};
class binding;
void dgPerfectGasLaw2dRegisterBindings(binding *b);
#endif
diff --git a/Solver/dgConservationLawShallowWater1d.cpp b/Solver/dgConservationLawShallowWater1d.cpp
index f342d9119d..88594dc2c3 100644
--- a/Solver/dgConservationLawShallowWater1d.cpp
+++ b/Solver/dgConservationLawShallowWater1d.cpp
@@ -9,7 +9,7 @@ class dgConservationLawShallowWater1d::advection: public dataCacheDouble {
public:
advection(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,9),
- sol(cacheMap.get("Solution",this))
+ sol(cacheMap.getSolution(this))
{};
void _eval () {
int nQP = _value.size1();
@@ -28,7 +28,7 @@ class dgConservationLawShallowWater1d::maxConvectiveSpeed : public dataCacheDoub
public:
maxConvectiveSpeed(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,1),
- sol(cacheMap.get("Solution",this))
+ sol(cacheMap.getSolution(this))
{
};
void _eval () {
@@ -45,9 +45,9 @@ class dgConservationLawShallowWater1d::source: public dataCacheDouble {
public :
source(dataCacheMap &cacheMap) :
dataCacheDouble(cacheMap,1,3),
- xyz(cacheMap.get("XYZ",this)),
- solution(cacheMap.get("Solution",this)),
- solutionGradient(cacheMap.get("SolutionGradient",this))
+ xyz(cacheMap.getParametricCoordinates(this)),
+ solution(cacheMap.getSolution(this)),
+ solutionGradient(cacheMap.getSolutionGradient(this))
{}
void _eval () {
int nQP =_value.size1();
@@ -64,9 +64,9 @@ class dgConservationLawShallowWater1d::riemann:public dataCacheDouble {
public:
riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
dataCacheDouble(cacheMapLeft,1,6),
- normals(cacheMapLeft.get("Normals", this)),
- solL(cacheMapLeft.get("Solution", this)),
- solR(cacheMapRight.get("Solution", this))
+ normals(cacheMapLeft.getNormals(this)),
+ solL(cacheMapLeft.getSolution(this)),
+ solR(cacheMapRight.getSolution(this))
{};
void _eval () {
int nQP = solL().size1();
@@ -102,8 +102,8 @@ class dgConservationLawShallowWater1d::boundaryWall : public dgBoundaryCondition
public:
term(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,3),
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)){}
+ sol(cacheMap.getSolution(this)),
+ normals(cacheMap.getNormals(this)){}
void _eval () {
int nQP = sol().size1();
for(int i=0; i< nQP; i++) {
@@ -130,7 +130,7 @@ class dgConservationLawShallowWater1d::clipToPhysics : public dataCacheDouble {
public:
clipToPhysics(dataCacheMap &cacheMap, double hMin):
dataCacheDouble(cacheMap,1,4),
- sol(cacheMap.get("Solution",this))
+ sol(cacheMap.getSolution(this))
{
_hMin=hMin;
};
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index 6def8ed8a4..b09a5c69aa 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -8,18 +8,18 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
class riemann;
class boundaryWall;
class maxConvectiveSpeed;
- std::string _linearDissipation, _quadraticDissipation, _source, _coriolisFactor, _bathymetry;
+ const function *_linearDissipation, *_quadraticDissipation, *_source, *_coriolisFactor, *_bathymetry;
public:
dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const;
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
dataCacheDouble *newMaxConvectiveSpeed (dataCacheMap &cacheMap) const;
- inline void setCoriolisFactor(std::string coriolisFactor){_coriolisFactor = coriolisFactor;}
- inline void setLinearDissipation(std::string linearDissipation){_linearDissipation = linearDissipation;}
- inline void setQuadraticDissipation(std::string quadraticDissipation){_quadraticDissipation = quadraticDissipation;}
- inline void setSource(std::string source){_source = source;}
- inline void setBathymetry(std::string bathymetry) {_bathymetry = bathymetry;}
+ inline void setCoriolisFactor(const function *coriolisFactor){_coriolisFactor = coriolisFactor;}
+ inline void setLinearDissipation(const function *linearDissipation){_linearDissipation = linearDissipation;}
+ inline void setQuadraticDissipation(const function *quadraticDissipation){_quadraticDissipation = quadraticDissipation;}
+ inline void setSource(const function *source){_source = source;}
+ inline void setBathymetry(const function *bathymetry) {_bathymetry = bathymetry;}
dgConservationLawShallowWater2d()
{
_nbf = 3; // eta u v
@@ -30,9 +30,9 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
class dgConservationLawShallowWater2d::maxConvectiveSpeed: public dataCacheDouble {
dataCacheDouble &sol, &_bathymetry;
public:
- maxConvectiveSpeed(dataCacheMap &cacheMap, std::string bathymetry):
+ maxConvectiveSpeed(dataCacheMap &cacheMap, const function *bathymetry):
dataCacheDouble(cacheMap,1,1),
- sol(cacheMap.get("Solution",this)),
+ sol(cacheMap.getSolution(this)),
_bathymetry(cacheMap.get(bathymetry,this))
{};
void _eval () {
@@ -50,9 +50,9 @@ class dgConservationLawShallowWater2d::maxConvectiveSpeed: public dataCacheDoubl
class dgConservationLawShallowWater2d::advection: public dataCacheDouble {
dataCacheDouble &sol, &_bathymetry;
public:
- advection(dataCacheMap &cacheMap, std::string bathymetry):
+ advection(dataCacheMap &cacheMap, const function *bathymetry):
dataCacheDouble(cacheMap,1,9),
- sol(cacheMap.get("Solution",this)),
+ sol(cacheMap.getSolution(this)),
_bathymetry(cacheMap.get(bathymetry,this))
{};
void _eval () {
@@ -82,14 +82,14 @@ class dgConservationLawShallowWater2d::source: public dataCacheDouble {
dataCacheDouble &solution,&solutionGradient;
dataCacheDouble &_coriolisFactor, &_source, &_linearDissipation, &_quadraticDissipation;
public :
- source(dataCacheMap &cacheMap, std::string linearDissipation, std::string quadraticDissipation, std::string coriolisFactor, std::string source) :
+ source(dataCacheMap &cacheMap, const function *linearDissipation, const function *quadraticDissipation, const function *coriolisFactor, const function *source) :
dataCacheDouble(cacheMap,1,3),
_coriolisFactor(cacheMap.get(coriolisFactor,this)),
_source(cacheMap.get(source,this)),
_linearDissipation(cacheMap.get(linearDissipation,this)),
_quadraticDissipation(cacheMap.get(quadraticDissipation,this)),
- solution(cacheMap.get("Solution",this)),
- solutionGradient(cacheMap.get("SolutionGradient",this))
+ solution(cacheMap.getSolution(this)),
+ solutionGradient(cacheMap.getSolutionGradient(this))
{
}
void _eval () {
@@ -112,11 +112,11 @@ class dgConservationLawShallowWater2d::source: public dataCacheDouble {
class dgConservationLawShallowWater2d::riemann:public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR, &_bathymetry;
public:
- riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, std::string bathymetry):
+ riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, const function *bathymetry):
dataCacheDouble(cacheMapLeft,1,6),
- normals(cacheMapLeft.get("Normals", this)),
- solL(cacheMapLeft.get("Solution", this)),
- solR(cacheMapRight.get("Solution", this)),
+ normals(cacheMapLeft.getNormals( this)),
+ solL(cacheMapLeft.getSolution( this)),
+ solR(cacheMapRight.getSolution( this)),
_bathymetry(cacheMapLeft.get(bathymetry, this))
{};
void _eval () {
@@ -185,8 +185,8 @@ class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition
public:
term(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,3),
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)){}
+ sol(cacheMap.getSolution(this)),
+ normals(cacheMap.getNormals(this)){}
void _eval () {
int nQP = sol().size1();
for(int i=0; i< nQP; i++) {
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index 2ff9fd9978..dc6b7aaacb 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -12,7 +12,7 @@ class dgConservationLawWaveEquation::hyperbolicFlux : public dataCacheDouble {
public:
hyperbolicFlux(dataCacheMap &cacheMap,int DIM):
dataCacheDouble(cacheMap,1,3*(DIM+1)),
- sol(cacheMap.get("Solution",this)),_DIM(DIM),_nbf(DIM+1)
+ sol(cacheMap.getSolution(this)),_DIM(DIM),_nbf(DIM+1)
{};
void _eval () {
int nQP = sol().size1();
@@ -32,7 +32,7 @@ class dgConservationLawWaveEquation::maxConvectiveSpeed : public dataCacheDouble
public:
maxConvectiveSpeed(dataCacheMap &cacheMap):
dataCacheDouble(cacheMap,1,1),
- sol(cacheMap.get("Solution",this))
+ sol(cacheMap.getSolution(this))
{
};
void _eval () {
@@ -46,9 +46,9 @@ class dgConservationLawWaveEquation::riemann : public dataCacheDouble {
public:
riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, int DIM):
dataCacheDouble(cacheMapLeft,1,2*(DIM+1)),
- normals(cacheMapLeft.get("Normals", this)),
- solL(cacheMapLeft.get("Solution", this)),
- solR(cacheMapRight.get("Solution", this)),
+ normals(cacheMapLeft.getNormals( this)),
+ solL(cacheMapLeft.getSolution( this)),
+ solR(cacheMapRight.getSolution( this)),
_DIM(DIM),_nbf(DIM+1)
{};
void _eval () {
@@ -105,8 +105,8 @@ class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition {
public:
term(dataCacheMap &cacheMap, int DIM):
dataCacheDouble(cacheMap,1,DIM+1),
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)),
+ sol(cacheMap.getSolution(this)),
+ normals(cacheMap.getNormals(this)),
_DIM(DIM){}
void _eval () {
int nQP = sol().size1();
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 645b42ec6d..390f7b9f51 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -217,7 +217,7 @@ void dgDofContainer::load(const std::string name) {
fclose(f);
}
-void dgDofContainer::L2Projection(std::string functionName){
+void dgDofContainer::L2Projection(const function *f){
scale(0.);
dgDofContainer rhs(&_groups, _nbFields);
for (int iGroup=0;iGroup<_groups.getNbElementGroups();iGroup++) {
@@ -227,7 +227,7 @@ void dgDofContainer::L2Projection(std::string functionName){
cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
dataCacheElement &cacheElement = cacheMap.getElement();
cacheMap.provideParametricCoordinates().set(group.getIntegrationPointsMatrix());
- dataCacheDouble &sourceTerm = cacheMap.get(functionName);
+ dataCacheDouble &sourceTerm = cacheMap.get(f);
fullMatrix<double> source;
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
cacheElement.set(group.getElement(iElement));
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index e4531cef00..97569432a0 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -6,6 +6,7 @@
#include "linearSystemCSR.h"
class dgGroupCollection;
class dgGroupOfElements;
+class function;
class dgDofContainer {
private:
@@ -42,7 +43,7 @@ public:
void save(const std::string name);
void load(const std::string name);
void setAll(double v);
- void L2Projection(std::string functionName);
+ void L2Projection(const function *f);
void Mesh2Mesh_BuildL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor);
void multAddInverseMassMatrixL2Projection(linearSystemCSRGmm<double> &projector); // this method should be private
void Mesh2Mesh_ApplyL2Projection(linearSystemCSRGmm<double> &projector,dgDofContainer &donor);
diff --git a/Solver/dgFunctionIntegrator.cpp b/Solver/dgFunctionIntegrator.cpp
index f04677d184..64bf289f52 100644
--- a/Solver/dgFunctionIntegrator.cpp
+++ b/Solver/dgFunctionIntegrator.cpp
@@ -6,7 +6,7 @@
#include <stdio.h>
-dgFunctionIntegrator::dgFunctionIntegrator(std::string fName):_fName(fName){}
+dgFunctionIntegrator::dgFunctionIntegrator(const function *f):_f(f){}
void dgFunctionIntegrator::compute(dgDofContainer *sol,fullMatrix<double> &result){
int nbFields=sol->getNbFields();
@@ -14,8 +14,7 @@ void dgFunctionIntegrator::compute(dgDofContainer *sol,fullMatrix<double> &resul
dataCacheDouble &UVW=cacheMap.provideParametricCoordinates();
dataCacheDouble &solutionQPe=cacheMap.provideSolution(nbFields);
dataCacheElement &cacheElement=cacheMap.getElement();
- function *f=function::get(_fName);
- dataCacheDouble *F=f->newDataCache(&cacheMap);
+ dataCacheDouble &F=cacheMap.get(_f);
int nbRowResult=result.size1();
result.scale(0.0);
for(int iGroup=0;iGroup<sol->getGroups()->getNbElementGroups();iGroup++){
@@ -32,7 +31,7 @@ void dgFunctionIntegrator::compute(dgDofContainer *sol,fullMatrix<double> &resul
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iElement, iPt);
for (int k=0;k<nbRowResult;k++){
- result(0,k) += (*F)(iPt,k)*detJ*IPMatrix(iPt,3);
+ result(0,k) += F(iPt,k)*detJ*IPMatrix(iPt,3);
}
}
}
@@ -44,8 +43,8 @@ void dgFunctionIntegrator::registerBindings(binding *b){
classBinding *cb = b->addClass<dgFunctionIntegrator>("dgFunctionIntegrator");
cb->setDescription("Integrates a function, using the compute function");
methodBinding *mb;
- mb = cb->setConstructor<dgFunctionIntegrator,std::string>();
- mb->setArgNames("functionName",NULL);
+ mb = cb->setConstructor<dgFunctionIntegrator,const function*>();
+ mb->setArgNames("function",NULL);
mb->setDescription("a new dgFunctionIntegrator, get the solution using the compute method");
mb = cb->addMethod("compute", &dgFunctionIntegrator::compute);
mb->setArgNames("dofContainer","result",NULL);
diff --git a/Solver/dgFunctionIntegrator.h b/Solver/dgFunctionIntegrator.h
index ebbb803764..3c3d0cbf1a 100644
--- a/Solver/dgFunctionIntegrator.h
+++ b/Solver/dgFunctionIntegrator.h
@@ -1,10 +1,11 @@
#include <string>
class dgDofContainer;
class binding;
+class function;
class dgFunctionIntegrator{
- std::string _fName;
+ const function *_f;
public:
- dgFunctionIntegrator(std::string fName);
+ dgFunctionIntegrator(const function *f);
void compute(dgDofContainer *sol,fullMatrix<double> &result);
static void registerBindings(binding *b);
};
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index d0e4fa66bf..5debcc304d 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -92,8 +92,8 @@ cm = cb->addMethod("RK44_TransformNodalValue",&dgSystemOfEquations::RK44_Transfo
}
// do a L2 projection
-void dgSystemOfEquations::L2Projection (std::string functionName){
- _solution->L2Projection(functionName);
+void dgSystemOfEquations::L2Projection (const function *f){
+ _solution->L2Projection(f);
}
// ok, we can setup the groups and create solution vectors
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index abc1a92abc..292515d244 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -41,7 +41,7 @@ public:
void setup (); // setup the groups and allocate
void exportSolution (std::string filename); // export the solution
void limitSolution (); // apply the limiter on the solution
- void L2Projection (std::string functionName); // assign the solution to a given function
+ void L2Projection (const function *f); // assign the solution to a given function
double RK44(double dt);
double RK44_limiter(double dt);
double RK44_TransformNodalValue(double dt);
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 9d4ce36541..ccec03aad6 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -53,7 +53,6 @@ function *function::get(std::string functionName, bool acceptNull)
}
return it->second;
}
-
//dataCacheMap members
dataCacheElement &dataCacheMap::getElement(dataCache *caller)
@@ -63,37 +62,60 @@ dataCacheElement &dataCacheMap::getElement(dataCache *caller)
return *_cacheElement;
}
-dataCacheDouble &dataCacheMap::get(const std::string &functionName, dataCache *caller)
+static dataCacheDouble &returnDataCacheDouble(dataCacheDouble *data, dataCache *caller)
{
- dataCacheDouble *&r= _cacheDoubleMap[functionName];
- if(r==NULL)
- r = function::get(functionName)->newDataCache(this);
+ if(data==NULL) throw;
if(caller)
- r->addMeAsDependencyOf(caller);
- return *r;
+ data->addMeAsDependencyOf(caller);
+ return *data;
+}
+dataCacheDouble &dataCacheMap::get(const function *f, dataCache *caller)
+{
+ dataCacheDouble *&r= _cacheDoubleMap[f];
+ if(r==NULL)
+ r = const_cast<function*>(f)->newDataCache(this);
+ return returnDataCacheDouble(r,caller);
+}
+dataCacheDouble &dataCacheMap::getSolution(dataCacheDouble *caller)
+{
+ return returnDataCacheDouble(_solution,caller);
}
+dataCacheDouble &dataCacheMap::getSolutionGradient(dataCacheDouble *caller)
+{
+ return returnDataCacheDouble(_solutionGradient,caller);
+}
+dataCacheDouble &dataCacheMap::getParametricCoordinates(dataCacheDouble *caller)
+{
+ return returnDataCacheDouble(_parametricCoordinates,caller);
+}
+dataCacheDouble &dataCacheMap::getNormals(dataCacheDouble *caller)
+{
+ returnDataCacheDouble(_normals,caller);
+ return *_normals;
+}
+
+
+
+
+
dataCacheDouble &dataCacheMap::provideSolution(int nbFields)
{
- dataCacheDouble *r = new providedDataDouble(*this,1, nbFields);
- _cacheDoubleMap["Solution"] = r;
- return *r;
+ _solution = new providedDataDouble(*this,1, nbFields);
+ return *_solution;
}
dataCacheDouble &dataCacheMap::provideSolutionGradient(int nbFields){
- dataCacheDouble *r = new providedDataDouble(*this,3, nbFields);
- _cacheDoubleMap["SolutionGradient"] = r;
- return *r;
+ _solutionGradient = new providedDataDouble(*this,3, nbFields);
+ return *_solutionGradient;
}
dataCacheDouble &dataCacheMap::provideParametricCoordinates()
{
- dataCacheDouble *r = new providedDataDouble(*this,1, 3);
- _cacheDoubleMap["UVW"] = r;
- return *r;
+ _parametricCoordinates = new providedDataDouble(*this,1, 3);
+ return *_parametricCoordinates;
}
dataCacheDouble &dataCacheMap::provideNormals()
{
- dataCacheDouble *r = new providedDataDouble(*this,1, 3);
- _cacheDoubleMap["Normals"] = r;
- return *r;
+ _normals = new providedDataDouble(*this,1, 3);
+ return *_normals;
}
dataCacheMap::~dataCacheMap()
@@ -107,17 +129,16 @@ dataCacheMap::~dataCacheMap()
// now some example of functions
// get XYZ coordinates
-class functionXYZ : public function {
- private:
+class functionCoordinates : public function {
+ static functionCoordinates *_instance;
class data : public dataCacheDouble{
- private:
dataCacheElement &_element;
dataCacheDouble &_uvw;
- int count;
- public:
+ int count;
+ public:
data(dataCacheMap *m) :
dataCacheDouble(*m, 1,3),
- _element(m->getElement(this)), _uvw(m->get("UVW", this))
+ _element(m->getElement(this)), _uvw(m->getParametricCoordinates(this))
{
}
void _eval()
@@ -133,17 +154,56 @@ class functionXYZ : public function {
~data(){
}
};
+ functionCoordinates(){};// constructor is private only 1 instance can exists, call get to access the instance
public:
dataCacheDouble *newDataCache(dataCacheMap *m)
{
return new data(m);
}
+ static function *get() {
+ if(!_instance)
+ _instance = new functionCoordinates();
+ return _instance;
+ }
+};
+functionCoordinates *functionCoordinates::_instance = NULL;
+
+class functionSolution : public function {
+ static functionSolution *_instance;
+ functionSolution(){};// constructor is private only 1 instance can exists, call get to access the instance
+ public:
+ dataCacheDouble *newDataCache(dataCacheMap *m)
+ {
+ return &m->getSolution(NULL);
+ }
+ static function *get() {
+ if(!_instance)
+ _instance = new functionSolution();
+ return _instance;
+ }
};
+functionSolution *functionSolution::_instance = NULL;
+
+class functionSolutionGradient : public function {
+ static functionSolutionGradient *_instance;
+ functionSolutionGradient(){};// constructor is private only 1 instance can exists, call get to access the instance
+ public:
+ dataCacheDouble *newDataCache(dataCacheMap *m)
+ {
+ return &m->getSolutionGradient(NULL);
+ }
+ static function *get() {
+ if(!_instance)
+ _instance = new functionSolutionGradient();
+ return _instance;
+ }
+};
+functionSolutionGradient *functionSolutionGradient::_instance = NULL;
class functionStructuredGridFile : public function {
public:
class data;
- std::string _coordFunction;
+ const function *_coordFunction;
int n[3];
double d[3],o[3];
double get(int i,int j, int k){
@@ -186,7 +246,7 @@ class functionStructuredGridFile : public function {
dataCacheDouble *newDataCache(dataCacheMap* m) {
return new data(this,m);
}
- functionStructuredGridFile(const std::string filename, const std::string coordFunction){
+ functionStructuredGridFile(const std::string filename, const function *coordFunction){
_coordFunction=coordFunction;
std::ifstream input(filename.c_str());
if(!input)
@@ -258,6 +318,7 @@ functionConstant::functionConstant(std::vector<double> source){
function::add(_name,this);
}
+
// function that enables to interpolate a DG solution using
// geometrical search in a mesh
@@ -277,7 +338,7 @@ class functionMesh2Mesh::data : public dataCacheDouble {
public:
data(dataCacheMap &m, dgDofContainer *sys) :
_dofContainer(sys),
- xyz(m.get("XYZ",this)),
+ xyz(m.get(functionCoordinates::get(),this)),
dataCacheDouble(m,1, sys->getNbFields())
{
}
@@ -311,6 +372,12 @@ public:
}
}
};
+
+dataCacheDouble *functionMesh2Mesh::newDataCache(dataCacheMap *m)
+{
+ return new data(*m,_dofContainer);
+}
+
void dataCacheMap::setNbEvaluationPoints(int nbEvaluationPoints) {
_nbEvaluationPoints = nbEvaluationPoints;
for(std::set<dataCacheDouble*>::iterator it = _toResize.begin(); it!= _toResize.end(); it++){
@@ -319,12 +386,6 @@ void dataCacheMap::setNbEvaluationPoints(int nbEvaluationPoints) {
}
}
-dataCacheDouble *functionMesh2Mesh::newDataCache(dataCacheMap *m)
-{
- 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;
@@ -338,7 +399,7 @@ void dataCacheDouble::resize() {
//functionC
class functionC : public function {
void (*callback)(void);
- std::vector<std::string> _dependenciesName;
+ std::vector<const function*> _dependenciesF;
int _nbCol;
class data : public dataCacheDouble{
const functionC *_function;
@@ -348,9 +409,9 @@ class functionC : public function {
dataCacheDouble(*m,1,f->_nbCol)
{
_function = f;
- _dependencies.resize ( _function->_dependenciesName.size());
- for (int i=0;i<_function->_dependenciesName.size();i++)
- _dependencies[i] = &m->get(_function->_dependenciesName[i],this);
+ _dependencies.resize ( _function->_dependenciesF.size());
+ for (int i=0;i<_function->_dependenciesF.size();i++)
+ _dependencies[i] = &m->get(_function->_dependenciesF[i],this);
}
void _eval()
@@ -394,8 +455,8 @@ class functionC : public function {
}
};
public:
- functionC (std::string file, std::string symbol, int nbCol, std::vector<std::string> dependencies):
- _dependenciesName(dependencies),_nbCol(nbCol)
+ functionC (std::string file, std::string symbol, int nbCol, std::vector<const function *> dependencies):
+ _dependenciesF(dependencies),_nbCol(nbCol)
{
void *dlHandler;
dlHandler = dlopen(file.c_str(),RTLD_NOW);
@@ -419,10 +480,6 @@ class functionC : public function {
#include "Bindings.h"
-void function::registerDefaultFunctions()
-{
- function::add("XYZ", new functionXYZ);
-}
void function::registerBindings(binding *b){
classBinding *cb = b->addClass<function>("function");
cb->setDescription("A generic function that can be evaluated on a set of points. Functions can call other functions and their values are cached so that if two different functions call the same function f, f is only evaluated once.");
@@ -437,10 +494,28 @@ void function::registerBindings(binding *b){
mb->setDescription("A new constant function wich values 'v' everywhere. v can be a row-vector.");
cb->setParentClass<function>();
+ cb = b->addClass<functionCoordinates>("functionCoordinates");
+ cb->setDescription("A function to access the coordinates (xyz). This is a single-instance class, use the 'get' member to access the instance.");
+ mb = cb->addMethod("get",&functionCoordinates::get);
+ mb->setDescription("return the unique instance of this class");
+ cb->setParentClass<function>();
+
+ cb = b->addClass<functionSolution>("functionSolution");
+ cb->setDescription("A function to access the solution. This is a single-instance class, use the 'get' member to access the instance.");
+ mb = cb->addMethod("get",&functionCoordinates::get);
+ mb->setDescription("return the unique instance of this class");
+ cb->setParentClass<function>();
+
+ cb = b->addClass<functionSolutionGradient>("functionSolutionGradient");
+ cb->setDescription("A function to access the gradient of the solution. This is a single-instance class, use the 'get' member to access the instance.");
+ mb = cb->addMethod("get",&functionCoordinates::get);
+ mb->setDescription("return the unique instance of this class");
+ cb->setParentClass<function>();
+
cb = b->addClass<functionStructuredGridFile>("functionStructuredGridFile");
cb->setParentClass<function>();
cb->setDescription("A function to interpolate through data given on a structured grid");
- mb = cb->setConstructor<functionStructuredGridFile,std::string, std::string>();
+ mb = cb->setConstructor<functionStructuredGridFile,std::string, const function*>();
mb->setArgNames("fileName","coordinateFunction",NULL);
mb->setDescription("Tri-linearly interpolate through data in file 'fileName' at coordinate given by 'coordinateFunction'.\nThe file format is :\nx0 y0 z0\ndx dy dz\nnx ny nz\nv(0,0,0) v(0,0,1) v(0 0 2) ...");
@@ -454,7 +529,7 @@ void function::registerBindings(binding *b){
#if defined(HAVE_DLOPEN)
cb = b->addClass<functionC>("functionC");
cb->setDescription("A function that compile a C code");
- mb = cb->setConstructor<functionC,std::string, std::string,int,std::vector<std::string> >();
+ mb = cb->setConstructor<functionC,std::string, std::string,int,std::vector<const function*> >();
mb->setArgNames("file", "symbol", "nbCol", "arguments",NULL);
mb->setDescription(" ");
cb->setParentClass<function>();
diff --git a/Solver/function.h b/Solver/function.h
index f87ae98e78..55ba50b137 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -12,24 +12,24 @@ class binding;
class dgDofContainer;
// those classes manage complex function dependencies and keep their values in cache so that they are not recomputed when it is not necessary. To do this, we use three classes : function, dataCache and dataCacheMap. The workflow is :
+// a) while parsing the input file and during the initialisation of the conservationLaw : all user-defined instance of function are inserted in the function map. (for example an user can create a function named "wind" of the class functionMathex with parameter "0.1*sin(xyz(0)/1e6); 0" and then give the string "wind" as parameter to it's conservation law to let the law know that this is the function to use as wind forcing)
//
-// a) before parsing the input file : a few general function that are always available are created (xyz for example).
+// b) before starting to iterate over the element :
+// b1) the algorithm create a new dataCacheMap instance.
+// b2) The algo insert in the dataCacheMap the dataCache it will provide (ie the quadrature points, the solution value,...)
+// b3) Then the algo give this dataCacheMap to the law and ask the source term of the law. This source term will be given as a dataCache, i.e. a cached value, a function to update it and a node in the dependency tree. In this example, to create this datacache two new dataCache will be required : one for the "wind" and one for "xyz". So, the source dataCache will ask the dataCacheMap for the dataCache called "wind". If this one does not already exist, the dataCacheMap will ask the function called "wind" to create it. In turn, the constructor of the dataCache associated with the "wind" will ask for the for the dataCache "xyz" wich will be created by the function "xyz" if it does not already exist and so on. The dataCacheMap contain also a special dataCache associated with the current element. During it's construction it's dataCache keep references to all the dataCache needed to compute it. When they request for the dataCache they need, the dataCacheMap is informed and the dependency tree is built.
//
-// b) while parsing the input file and during the initialisation of the conservationLaw : all user-defined instance of function are inserted in the function map. (for example an user can create a function named "wind" of the class functionMathex with parameter "0.1*sin(xyz(0)/1e6); 0" and then give the string "wind" as parameter to it's conservation law to let the law know that this is the function to use as wind forcing)
+// c) When iterating over elements, for each element : the Algo will change the dataCache element so that it point to the current element. This will invalidate everything that depend on the current on the current element (in this case, xyz, wind (because it depends on xyz) and the source dataCache (because it depends on wind). So, in this case everything (exept the current element dataCache) is marked as invalid. Now the algo want to access the value of the 'source' dataCache this one is invalid, so it is recomputed. To recompute it, the value of the wind is accessed and it's value is recomputed, and so on. Now, if some other function access the wind or the position (xyz), they are already marked as valid and are not re-computed.
//
-// c) before starting to iterate over the element :
-// c1) the algorithm create a new dataCacheMap instance.
-// c2) The algo insert in the dataCacheMap the dataCache it will provide (ie the quadrature points, the solution value,...)
-// c3) Then the algo give this dataCacheMap to the law and ask the source term of the law. This source term will be given as a dataCache, i.e. a cached value, a function to update it and a node in the dependency tree. In this example, to create this datacache two new dataCache will be required : one for the "wind" and one for "xyz". So, the source dataCache will ask the dataCacheMap for the dataCache called "wind". If this one does not already exist, the dataCacheMap will ask the function called "wind" to create it. In turn, the constructor of the dataCache associated with the "wind" will ask for the for the dataCache "xyz" wich will be created by the function "xyz" if it does not already exist and so on. The dataCacheMap contain also a special dataCache associated with the current element. During it's construction it's dataCache keep references to all the dataCache needed to compute it. When they request for the dataCache they need, the dataCacheMap is informed and the dependency tree is built.
-//
-// d) When iterating over elements, for each element : the Algo will change the dataCache element so that it point to the current element. This will invalidate everything that depend on the current on the current element (in this case, xyz, wind (because it depends on xyz) and the source dataCache (because it depends on wind). So, in this case everything (exept the current element dataCache) is marked as invalid. Now the algo want to access the value of the 'source' dataCache this one is invalid, so it is recomputed. To recompute it, the value of the wind is accessed and it's value is recomputed, and so on. Now, if some other function access the wind or the position (xyz), they are already marked as valid and are not re-computed.
-//
-// e) after the loop over the elements the dataCacheMap is deleted and it delete all its dataCaches.
+// d) after the loop over the elements the dataCacheMap is deleted and it delete all its dataCaches.
//
// NB) in the case of integration over faces, there is two dataCacheMap : one for the Right side and one for the left side i.e. two sets of cached values and two dependency trees. Some dataCache (ie the value of the flux) have dependencies in both tree. This is not a problem, the only difference with other dataCaches is that they are not owned by the dataCacheMap and have to be deleted manually.
//
// a node in the dependency tree. The usefull field is _dependOnMe which is the set of every other nodes that depend on me. When the value of this node change all nodes depending on this one are marked as "invalid" and will be recomputed the next time their data are accessed. To be able to maintain _dependOnMe up to date when a new node is inserted in the tree, we need _iDependOn list. So we do not really store a tree but instead each node contain a complete list of all it's parents and all it's children (and the parents of the parents of ... of its parents and the children of the children of ... of it's children). This way invalidate all the dependencies of a node is really fast and does not involve a complex walk accross the tree structure.
+
+class function;
+
class dataCache {
friend class dataCacheMap;
// pointers to all of the dataCache depending on me
@@ -131,7 +131,6 @@ class function {
protected:
std::string _name;
public:
- static void registerDefaultFunctions();
static bool add(const std::string functionName, function *f);
static function *get(std::string functionName, bool acceptNull=false);
@@ -166,7 +165,7 @@ class dataCacheMap {
// keep track of the current element and all the dataCaches that
// depend on it
dataCacheElement *_cacheElement;
- std::map<std::string, dataCacheDouble*> _cacheDoubleMap;
+ std::map<const function*, dataCacheDouble*> _cacheDoubleMap;
class providedDataDouble : public dataCacheDouble
// for data provided by the algorithm and that does not have an _eval function
// (typically "UVW")
@@ -194,15 +193,21 @@ class dataCacheMap {
//list of dgDofContainer whom gradient are needed
std::map<dgDofContainer*,dataCacheDouble*> gradientFields;
// end dg Part
+ dataCacheDouble *_solution, *_solutionGradient, *_parametricCoordinates, *_normals;
+ dataCacheDouble &getSolution(dataCacheDouble *caller);
+ dataCacheDouble &getSolutionGradient(dataCacheDouble *caller);
+ dataCacheDouble &getParametricCoordinates(dataCacheDouble *caller);
+ dataCacheDouble &getNormals(dataCacheDouble *caller);
dataCacheDouble &provideSolution(int nbFields);
dataCacheDouble &provideSolutionGradient(int nbFields);
dataCacheDouble &provideParametricCoordinates();
dataCacheDouble &provideNormals();
- dataCacheDouble &get(const std::string &functionName, dataCache *caller=0);
+ dataCacheDouble &get(const function *f, dataCache *caller=0);
dataCacheElement &getElement(dataCache *caller=0);
dataCacheMap(){
_cacheElement = new dataCacheElement(this);
+ _normals = _solution = _solutionGradient = _parametricCoordinates = 0;
_nbEvaluationPoints = 0;
}
void setNbEvaluationPoints(int nbEvaluationPoints);
@@ -226,5 +231,4 @@ public:
dataCacheDouble *newDataCache(dataCacheMap *m);
functionMesh2Mesh(dgDofContainer *dofc) ;
};
-
#endif
diff --git a/Solver/luaFunction.cpp b/Solver/luaFunction.cpp
index f53b0a979c..f52cbc968d 100644
--- a/Solver/luaFunction.cpp
+++ b/Solver/luaFunction.cpp
@@ -17,9 +17,9 @@ class functionLua::data : public dataCacheDouble{
dataCacheDouble(*m,1,f->_nbCol),
_function(f)
{
- _dependencies.resize ( _function->_dependenciesName.size());
- for (int i=0;i<_function->_dependenciesName.size();i++)
- _dependencies[i] = &m->get(_function->_dependenciesName[i],this);
+ _dependencies.resize ( _function->_dependenciesF.size());
+ for (int i=0;i<_function->_dependenciesF.size();i++)
+ _dependencies[i] = &m->get(_function->_dependenciesF[i],this);
}
void _eval()
{
@@ -32,8 +32,8 @@ class functionLua::data : public dataCacheDouble{
lua_call(_function->_L,_dependencies.size()+1,0); /* call Lua function */
}
};
-functionLua::functionLua (int nbCol, std::string luaFunctionName, std::vector<std::string> dependencies, lua_State *L)
- : _luaFunctionName(luaFunctionName), _dependenciesName(dependencies),_L(L),_nbCol(nbCol)
+functionLua::functionLua (int nbCol, std::string luaFunctionName, std::vector<const function*> dependencies, lua_State *L)
+ : _luaFunctionName(luaFunctionName), _dependenciesF(dependencies),_L(L),_nbCol(nbCol)
{
static int c=0;
std::ostringstream oss;
@@ -50,7 +50,7 @@ void functionLua::registerBindings(binding *b){
classBinding *cb= b->addClass<functionLua>("functionLua");
cb->setDescription("A function (see the 'function' documentation entry) defined in LUA.");
methodBinding *mb;
- mb = cb->setConstructor<functionLua,int,std::string,std::vector<std::string>,lua_State*>();
+ mb = cb->setConstructor<functionLua,int,std::string,std::vector<const function*>,lua_State*>();
mb->setArgNames("d", "f", "dep", NULL);
mb->setDescription("A new functionLua which evaluates a vector of dimension 'd' using the lua function 'f'. This function can take other functions as arguments listed by the 'dep' vector.");
cb->setParentClass<function>();
diff --git a/Solver/luaFunction.h b/Solver/luaFunction.h
index 0ffb9456dd..4dd73bceba 100644
--- a/Solver/luaFunction.h
+++ b/Solver/luaFunction.h
@@ -10,11 +10,11 @@ class binding;
class functionLua : public function {
lua_State *_L;
std::string _luaFunctionName;
- std::vector<std::string> _dependenciesName;
+ std::vector<const function*> _dependenciesF;
int _nbCol;
class data;
public:
- functionLua (int nbCol, std::string luaFunctionName, std::vector<std::string> dependenciesName, lua_State *L);
+ functionLua (int nbCol, std::string luaFunctionName, std::vector<const function*> dependencies, lua_State *L);
dataCacheDouble *newDataCache(dataCacheMap *m);
static void registerBindings(binding *b);
--
GitLab