diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 3b12907cd20ebffe709244c831cd892b94201c67..d1523dc3deb62845e0cc6e8d7003890407b42c9a 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 c6f26566ef96a27d19d65cef281fd370597d8130..f4a65368f1fc563af44bfdba8d08dd94eed913db 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 9c9ab2873352e4b7e54bf44922f6026f138fc431..2b5eeccf38d078828d96b2c4440bbb35ef3e7a6a 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 13b691ddc757c4359d269d27cea0d262e210c657..a27a3150521c1eb1ee455ae0465c51f4789de175 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 71ac0ad30ddc68453fa16f87d614e2e2553f3fb0..ad094513a8ac4ab756293ecc31c171ce588d431d 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 04fba32a3aa2e241aba8c45bf9ddb18a9e86f1e8..26be9b8cf9ab24689b5f55bcb80ee44d350dc5c0 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 d2bd216e0f9333a241778610a256ad5dde9b4df1..d9302dde173ceace7f1cd33eb744cd609610c585 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 c6f26566ef96a27d19d65cef281fd370597d8130..a15fdf96163a41802d6ba7ba95b8713fe7ec9425 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 abf15298e4c4b6917b57aaa5f6432b9d35e86e51..16cf5b4107977a298cd1138d72b2f5cd5dd4a0b6 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 de184a8a2a32eb85745e31984303e2b4daadc128..2b66374db7030e0adf49e3a39325566e09840d80 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 37c46b7d931ffa182a8597fac4407e6399bf4469..2e8c790735b1e0507e0b9b4f84904253befa8076 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 e9bc0468546ee6280693f5b8da738e761891eb67..2b32ff8bdcac6d58aa71275e6e2d314c6aeae56e 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 65316de0587f55601f4b1ed7252539eb89448231..758acf140b6dc3ec115e3cee0da97458d42aee4a 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 9ca9c3846e83e65f5e592790bcb3ab25d7069bc9..49743b6da85efb5799895c9a0841676573c512c7 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 dd8457e3811eb38276eec524a86012594bb4ada8..aa4e15ac6456b49f3c148faea8d21c20a16c7e4d 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 7dab96977826844e784de4c7e960d4854e1ca0d8..834d761c59d0a4a5d5483edc7c3bf615f75a2472 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 3cbe06c1c44f0dfd55e7f596e2c006fb78f255fe..df0c4904f6986d049de5450ec37306d4fb5e3de2 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 f342d9119dee04f58675bb1980dc6c7b177c6e9d..88594dc2c3e5d06ef65c852837aa8d9117a1fc7c 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 6def8ed8a4b106c43ab30a6b15d07e14caefd379..b09a5c69aadea1ed2c6a0e0c19ddfd6fd20d137c 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 2ff9fd997803ee835381fec8973e5140a8be6e2d..dc6b7aaacb44621fb3b2580f3c52a85f6a9c7301 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 645b42ec6de422c776977f78139474e9e2410f20..390f7b9f51e9ab0b2320ac316632c08c055e5ee9 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 e4531cef00a6a2135582e533061c13b959d63d46..97569432a0465ac03718cae7b29296475a89c3a5 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 f04677d184206d411e0cd165aaf1d206d6157ec4..64bf289f52375474d395615a130afbf73bfe8b8f 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 ebbb803764c6968e0b004025ca2b4974313547aa..3c3d0cbf1a024f00426137fd3faab7c9127e996e 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 d0e4fa66bfeb6c7b312314cfd3055daeb09950fe..5debcc304d61940db20345a5519a9927d1f70c4e 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 abc1a92abc79bff359b69b0923e00b5899d4ef0b..292515d2449b66e9f991294328a76cff76b6c448 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 9d4ce365416b58d775ec25a2265308a95ecbc8b0..ccec03aad6df2786d33064006b09844eb1bfebb6 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 f87ae98e78302c4213d58943c0bc8c9fbdad24a6..55ba50b13778ea5ee368df31d69a5c4c43389f13 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 f53b0a979caaf470677442c13e025e1b3e2bf44a..f52cbc968d9831a6e180882412a89afa2ed68e93 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 0ffb9456dddf08b75755990e3e56150881624019..4dd73bceba6a415318d4ed33e2fb04ce41ef8865 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);