diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 694438eda9a5a984a2d31052d0c93a38f64d8db4..6d74439f89799caab25ef00ef62645192e83d27e 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -13,6 +13,7 @@ set(SRC
multiscaleLaplace.cpp
dgGroupOfElements.cpp
dgAlgorithm.cpp
+ dgConservationLaw.cpp
dgConservationLawAdvection.cpp
function.cpp
)
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index ffef08e52305f9c26abbd03127684a4379ccfe60..eb97d9578aef0d1641f07bec07d095480a6df0fc 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -25,13 +25,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
// ----- 1.2 --- multiply the solution by the collocation matrix
group.getCollocationMatrix().mult(solution , solutionQP);
- // ----- 1.3 --- if the conservation law is diffusive, compute the gradients too
- fullMatrix<double> *gradientSolutionQP= 0;
- if (claw.diffusiveFlux()){
- gradientSolutionQP = new fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields * 3);
- //group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP);
- }
-
+
// ----- 2 ---- compute all fluxes (diffusive and convective) at integration points
// ----- 2.1 --- allocate elementary fluxes (compued in XYZ coordinates)
fullMatrix<double> fConv[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), nbFields ),
@@ -47,22 +41,30 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields)};
dataCacheMap cacheMap;
dataCacheElement &cacheElement = cacheMap.getElement();
+ // provided dataCache
cacheMap.provideData("UVW").set(group.getIntegrationPointsMatrix());
+ dataCacheDouble &solutionQPe = cacheMap.provideData("Solution");
+ dataCacheDouble &gradSolutionQPe = cacheMap.provideData("gradSolution");
+ // needed dataCache
dataCacheDouble *sourceTerm = claw.newSourceTerm(cacheMap);
+ // fluxes in XYZ coordinates;
+ dataCacheDouble *convectiveFlux = claw.newConvectiveFlux(cacheMap);
+ dataCacheDouble *diffusiveFlux = claw.newDiffusiveFlux(cacheMap);
+ // compute the gradient of solution if necessary
+ fullMatrix<double> *gradientSolutionQP= 0;
+ if (gradSolutionQPe.somethingDependOnMe()) {
+ gradientSolutionQP = new fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields * 3);
+ //group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP);
+ }
// ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
- fullMatrix<double> solutionQPe (solutionQP, iElement*nbFields, nbFields );
- fullMatrix<double> *gradSolutionQPe;
- if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradientSolutionQP, 3*iElement*nbFields,3*nbFields );
- else gradSolutionQPe = new fullMatrix<double>;
- dgElement DGE( group.getElement(iElement), solutionQPe, *gradSolutionQPe, group.getIntegrationPointsMatrix());
+ solutionQPe.set(fullMatrix<double>(solutionQP, iElement*nbFields, nbFields ));
+ if (gradSolutionQPe.somethingDependOnMe())
+ gradSolutionQPe.set(fullMatrix<double>(*gradientSolutionQP, 3*iElement*nbFields,3*nbFields ));
cacheElement.set(group.getElement(iElement));
- if(claw.convectiveFlux() || claw.diffusiveFlux()) {
- // ----- 2.3.2 --- compute fluxes in XYZ coordinates
- if (claw.convectiveFlux()) (*claw.convectiveFlux())(DGE,fConv);
- if (claw.diffusiveFlux()) (*claw.diffusiveFlux())(DGE,fDiff);
+ if(convectiveFlux || diffusiveFlux) {
// ----- 2.3.3 --- compute fluxes in UVW coordinates
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++) {
// ----- 2.3.3.1 --- get a proxy on the big local flux matrix
@@ -76,8 +78,12 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
const double detJ = group.getDetJ (iElement, iPt);
const double factor = invJ * detJ;
// compute fluxes in the reference coordinate system at this integration point
- for (int k=0;k<nbFields;k++)
- fuvwe(iPt,k) += ( fConv[iXYZ](iPt,k) + fDiff[iXYZ](iPt,k)) * factor;
+ for (int k=0;k<nbFields;k++){
+ if(convectiveFlux)
+ fuvwe(iPt,k) += ( (*convectiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
+ if(diffusiveFlux)
+ fuvwe(iPt,k) += ( (*diffusiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
+ }
}
}
}
@@ -91,13 +97,18 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
}
}
}
- delete gradSolutionQPe;
+ if(gradientSolutionQP)
+ delete gradientSolutionQP;
}
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
- if(claw.convectiveFlux() || claw.diffusiveFlux()){
+ if(convectiveFlux || diffusiveFlux){
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){
residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
}
+ if(convectiveFlux)
+ delete convectiveFlux;
+ if(diffusiveFlux)
+ delete diffusiveFlux;
}
if(sourceTerm){
residual.gemm(group.getSourceRedistributionMatrix(),Source);
@@ -118,23 +129,34 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
group.getCollocationMatrix().mult(solution, solutionQP);
// ----- 2 ---- compute normal fluxes at integration points
fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields*2);
+ dataCacheMap cacheMapLeft, cacheMapRight;
+ // provided dataCache
+ cacheMapLeft.provideData("UVW").set(group.getIntegrationPointsMatrix());
+ dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
+ dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution");
+ dataCacheDouble &normals = cacheMapLeft.provideData("Normals");
+ dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
+ dataCacheElement &cacheElementRight = cacheMapRight.getElement();
+ // required data
+ dataCacheDouble *riemannSolver = claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
- fullMatrix<double> solutionQPLeft (solutionQP, iFace*2*nbFields, nbFields );
- fullMatrix<double> solutionQPRight (solutionQP, (iFace*2+1)*nbFields, nbFields );
- fullMatrix<double> normalFluxQP (NormalFluxQP, iFace*2*nbFields, nbFields*2);
- dgFace DGF( group.getFace(iFace), group.getElementLeft(iFace), group.getElementRight(iFace),
- solutionQPLeft, solutionQPRight, group.getIntegrationPointsMatrix(),group.getNormals(iFace));
+ solutionQPLeft.set(fullMatrix<double>(solutionQP, iFace*2*nbFields, nbFields));
+ solutionQPRight.set(fullMatrix<double>(solutionQP, (iFace*2+1)*nbFields, nbFields));
+ normals.set(group.getNormals(iFace));
+ cacheElementLeft.set(group.getElementLeft(iFace));
+ cacheElementRight.set(group.getElementRight(iFace));
+ fullMatrix<double> normalFluxQP (NormalFluxQP, iFace*nbFields*2, nbFields*2);
// ----- 2.3.2 --- compute fluxes
- (*claw.riemannSolver())(DGF,&normalFluxQP);
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
for (int k=0;k<nbFields*2;k++)
- normalFluxQP(iPt,k) *= detJ;
+ normalFluxQP(iPt,k) = (*riemannSolver)(iPt,k)*detJ;
}
}
// ----- 3 ---- do the redistribution at face nodes using BLAS3
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
+ delete riemannSolver;
}
void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
@@ -157,40 +179,39 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
)
{
int nbFields = claw.nbFields();
- const dgBoundaryCondition *boundaryCondition = claw.boundaryCondition(group.getBoundaryTag());
+ const dgBoundaryCondition *boundaryCondition = claw.getBoundaryCondition(group.getBoundaryTag());
// ----- 1 ---- get the solution at quadrature points
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
group.getCollocationMatrix().mult(solution, solutionQP);
// ----- 2 ---- compute normal fluxes at integration points
fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields);
- fullMatrix<double> solutionQPRight(group.getNbIntegrationPoints(), nbFields);
+
+ dataCacheMap cacheMapLeft;
+ // provided dataCache
+ cacheMapLeft.provideData("UVW").set(group.getIntegrationPointsMatrix());
+ dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
+ dataCacheDouble &normals = cacheMapLeft.provideData("Normals");
+ dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
+ // required data
+ dataCacheDouble *boundaryTerm = boundaryCondition->newBoundaryTerm(cacheMapLeft);
+
for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
- // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
- fullMatrix<double> solutionQPLeft (solutionQP, iFace*nbFields, nbFields);
- fullMatrix<double> normalFluxQP (NormalFluxQP, iFace*nbFields, nbFields);
- dgFace DGF( group.getFace(iFace), group.getElementLeft(iFace), NULL,
- solutionQPLeft, solutionQPRight, group.getIntegrationPointsMatrix(),group.getNormals(iFace));
+ // ----- 2.3.1 --- provide the data to the cacheMap
+ solutionQPLeft.set(fullMatrix<double>(solutionQP, iFace*nbFields, nbFields));
+ normals.set(group.getNormals(iFace));
+ cacheElementLeft.set(group.getElementLeft(iFace));
+
// ----- 2.3.2 --- compute fluxes
- switch (boundaryCondition->type()){
- case dgBoundaryCondition::EXTERNAL_VALUES : {
- (*boundaryCondition)(DGF, &solutionQPRight);
- (*claw.riemannSolver())(DGF, &normalFluxQP);
- break;
- }
- case dgBoundaryCondition::FLUX :{
- (*boundaryCondition)(DGF, &normalFluxQP);
- }
- default :
- throw;
- }
+ fullMatrix<double> normalFluxQP (NormalFluxQP, iFace*nbFields, nbFields);
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
for (int k=0;k<nbFields;k++)
- normalFluxQP(iPt,k) *= detJ;
+ normalFluxQP(iPt,k) = (*boundaryTerm)(iPt,k)*detJ;
}
}
// ----- 3 ---- do the redistribution at face nodes using BLAS3
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
+ delete boundaryTerm;
}
void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f979df2e5826d0756e801527321050de0f2d078c
--- /dev/null
+++ b/Solver/dgConservationLaw.cpp
@@ -0,0 +1,65 @@
+#include "dgConservationLaw.h"
+#include "function.h"
+class dgBoundaryCondition0Out : public dgBoundaryCondition {
+ dgConservationLaw &_claw;
+ class term : public dataCacheDouble {
+ dataCacheMap cacheMapRight; // new cacheMap to pass to the Riemann solver
+ dataCacheDouble &solutionRight;
+ dataCacheDouble &solutionLeft;
+ dataCacheDouble *riemannSolver;
+ dgConservationLaw &_claw;
+ public:
+ term(dgConservationLaw &claw, dataCacheMap &cacheMapLeft):
+ solutionRight(cacheMapRight.provideData("Solution")),
+ solutionLeft(cacheMapLeft.get("Solution",this)),
+ _claw(claw)
+ {
+ riemannSolver=_claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
+ }
+
+ void _eval() {
+ if(_value.size1()!=solutionLeft().size1()){
+ //adjust sizes
+ solutionRight.set(fullMatrix<double>(solutionLeft().size1(),_claw.nbFields()));
+ _value = fullMatrix<double>(solutionLeft().size1(),_claw.nbFields());
+ }
+ for(int i=0;i<_value.size1(); i++)
+ for(int j=0;j<_value.size2(); j++)
+ _value(i,j) = (*riemannSolver)(i,j*2);
+ }
+ };
+ public:
+ dgBoundaryCondition0Out(dgConservationLaw &claw): _claw(claw) {}
+ dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+ return new term(_claw,cacheMapLeft);
+ }
+};
+
+class dgBoundaryCondition0Flux : public dgBoundaryCondition {
+ dgConservationLaw &_claw;
+ class term : public dataCacheDouble {
+ dataCacheDouble &UVW;
+ dgConservationLaw &_claw;
+ public:
+ term(dgConservationLaw &claw,dataCacheMap &cacheMapLeft): UVW(cacheMapLeft.get("UVW",this)),_claw(claw) {}
+ void _eval() {
+ if(_value.size1()!=UVW().size1()){
+ //adjust sizes
+ _value = fullMatrix<double>(UVW().size1(),_claw.nbFields());
+ }
+ }
+ };
+ public:
+ dgBoundaryCondition0Flux(dgConservationLaw &claw): _claw(claw) {}
+ dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+ return new term(_claw,cacheMapLeft);
+ }
+};
+
+dgBoundaryCondition *dgBoundaryCondition::new0OutCondition(dgConservationLaw &claw) {
+ return new dgBoundaryCondition0Out(claw);
+}
+dgBoundaryCondition *dgBoundaryCondition::new0FluxCondition(dgConservationLaw &claw) {
+ return new dgBoundaryCondition0Flux(claw);
+}
+
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 481d125f23aaed14e338942f17df178b2f642864..9e72ab9a3745e85d501afb94d67b92df1c131655 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -7,59 +7,50 @@
+ r(u,forcings) -> source term r
*/
#include "fullMatrix.h"
-class dgElement;
-class dgFace;
class dataCacheDouble;
+class dataCacheMap;
-class dgTerm{
- public:
- virtual ~dgTerm () {}
- virtual void operator () (const dgElement &, fullMatrix<double> fcx[]) const = 0;
-};
-
-class dgFaceTerm{
- public:
- virtual ~dgFaceTerm () {}
- virtual void operator () (const dgFace &, fullMatrix<double> fcx[]) const = 0;
-};
+class dgConservationLaw;
class dgBoundaryCondition {
-public:
+ public:
virtual ~dgBoundaryCondition () {}
- typedef enum {FLUX=0,EXTERNAL_VALUES}boundaryType;
- virtual void operator () (const dgFace&, fullMatrix<double> fcx[]) const = 0;
- virtual boundaryType type() const =0;
+ virtual dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const = 0;
+ //a generic boundary condition using the Riemann solver of the conservation Law
+ static dgBoundaryCondition *new0OutCondition(dgConservationLaw &claw);
+ static dgBoundaryCondition *new0FluxCondition(dgConservationLaw &claw);
};
-class dataCacheMap;
class dgConservationLaw {
- protected :
- int _nbf;
- dgTerm *_diffusive, *_convective, *_source, *_maxConvectiveSpeed;
- dgFaceTerm *_riemannSolver;
std::map<const std::string,dgBoundaryCondition*> _boundaryConditions;
-public:
- dgConservationLaw () : _diffusive(0), _convective(0), _source (0),
- _riemannSolver(0),_maxConvectiveSpeed (0) {}
- ~dgConservationLaw () {}
+ protected :
+ int _nbf;
+ public:
+ virtual ~dgConservationLaw () {}
int nbFields() const {return _nbf;}
- virtual dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const { return NULL; }
+ virtual dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const {return NULL;}
+ virtual dataCacheDouble *newDiffusiveFlux (dataCacheMap &cacheMap) const {return NULL;}
+ virtual dataCacheDouble *newConvectiveFlux (dataCacheMap &cacheMap) const {return NULL;}
+ virtual dataCacheDouble *newMaxConvectiveSpeed (dataCacheMap &cacheMap) const {return NULL;}
+ virtual dataCacheDouble *newRiemannSolver (dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {return NULL;}
- inline const dgTerm * convectiveFlux () const {return _convective;}
- inline const dgTerm * diffusiveFlux () const {return _diffusive;}
- inline const dgFaceTerm * riemannSolver () const {return _riemannSolver;}
- inline const dgTerm * maxConvectiveSpeed () const {return _maxConvectiveSpeed;}
- inline const dgBoundaryCondition *boundaryCondition(const std::string tag) const {
+ inline const dgBoundaryCondition *getBoundaryCondition(const std::string tag) const {
std::map<const std::string,dgBoundaryCondition*>::const_iterator it = _boundaryConditions.find(tag);
if(it==_boundaryConditions.end())
throw;
return it->second;
}
+ inline void addBoundaryCondition(const std::string tag, dgBoundaryCondition * condition) {
+ if(_boundaryConditions.find(tag)!=_boundaryConditions.end())
+ throw;
+ _boundaryConditions[tag]=condition;
+ }
+
};
-dgConservationLaw *dgNewConservationLawAdvection();
+dgConservationLaw *dgNewConservationLawAdvection(const std::string vname);
#endif
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index 8f4677c9bc248eba70297de1fae432bcf19b7531..a886ef5bf8aeb7c0eb2c3e1088eb1485b8de8d7d 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -3,95 +3,70 @@
#include "dgGroupOfElements.h"
#include "SPoint3.h"
#include "MElement.h"
-void getV(SPoint3 &p, double (&v)[3]){
- /*double r=hypot(p.x(),p.y());
- double alpha=atan2(p.y(),p.x());
- v[0]=r*sin(alpha);
- v[1]=-r*cos(alpha);
- v[2]=0;*/
- v[0]=0.15;
- v[1]=0;
- v[2]=0;
-}
+#include "function.h"
-class dgAdvectionFluxTerm : public dgTerm {
- void operator () (const dgElement &el, fullMatrix<double> fcx[]) const
- {
- const fullMatrix<double> &sol = el.solution();
- const fullMatrix<double> &QP = el.integration();
- for(int i=0; i< sol.size1(); i++) {
- SPoint3 p;
- el.element()->pnt(QP(i,0),QP(i,1),QP(i,2),p);
- double v[3];
- getV(p,v);
- fcx[0](i,0) = sol(i,0)*v[0];
- fcx[1](i,0) = sol(i,0)*v[1];
- fcx[2](i,0) = sol(i,0)*v[2];
+class dgConservationLawAdvection : public dgConservationLaw {
+ std::string _vFunctionName;
+ class advection : public dataCacheDouble {
+ dataCacheDouble &sol, &uvw, &v;
+ public:
+ advection(std::string vFunctionName, dataCacheMap &cacheMap):
+ uvw(cacheMap.get("UVW",this)),
+ sol(cacheMap.get("Solution",this)),
+ v(cacheMap.get(vFunctionName,this))
+ {};
+ void _eval () {
+ if(_value.size1() !=sol().size1())
+ _value=fullMatrix<double>(sol().size1(),3);
+ for(int i=0; i< sol().size1(); i++) {
+ _value(i,0) = sol(i,0)*v(i,0);
+ _value(i,1) = sol(i,0)*v(i,1);
+ _value(i,2) = sol(i,0)*v(i,2);
+ }
}
- }
-};
-
-class dgAdvectionRiemanTerm : public dgFaceTerm{
- public:
- void operator () (const dgFace &face, fullMatrix<double> fcx[]) const
- {
- const fullMatrix<double> &solLeft = face.solutionLeft();
- const fullMatrix<double> &solRight = face.solutionRight();
- const fullMatrix<double> &normals = face.normals();
- const fullMatrix<double> &QP = face.integration();
- for(int i=0; i< solLeft.size1(); i++) {
- SPoint3 p;
- face.face()->pnt(QP(i,0),QP(i,1),QP(i,2),p);
- double v[3];
- getV(p,v);
- double un=v[0]*normals(0,i)+v[1]*normals(1,i)+v[2]*normals(2,i);
- if(un>0){
- fcx[0](i,0) = -solLeft(i,0)*un;
- if(fcx[0].size2()==2)
- fcx[0](i,1) = solLeft(i,0)*un;
- }else{
- fcx[0](i,0) = -solRight(i,0)*un;
- if(fcx[0].size2()==2)
- fcx[0](i,1) = solRight(i,0)*un;
+ };
+ class riemann : public dataCacheDouble {
+ dataCacheDouble &uvw, &normals, &solLeft, &solRight,&v;
+ public:
+ riemann(std::string vFunctionName, dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+ uvw(cacheMapLeft.get("UVW", this)),
+ normals(cacheMapLeft.get("Normals", this)),
+ solLeft(cacheMapLeft.get("Solution", this)),
+ solRight(cacheMapRight.get("Solution", this)),
+ v(cacheMapLeft.get(vFunctionName,this))
+ {};
+ void _eval () {
+ if(_value.size1() !=solLeft().size1())
+ _value = fullMatrix<double>(solLeft().size1(),2);
+ for(int i=0; i< _value.size1(); i++) {
+ double un=v(i,0)*normals(0,i)+v(i,1)*normals(1,i)+v(i,2)*normals(2,i);
+ if(un>0){
+ _value(i,0) = -solLeft(i,0)*un;
+ _value(i,1) = solLeft(i,0)*un;
+ }else{
+ _value(i,0) = -solRight(i,0)*un;
+ _value(i,1) = solRight(i,0)*un;
+ }
}
}
+ };
+ public:
+ dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const {
+ return new advection(_vFunctionName,cacheMap);
}
-};
-class dgAdvectionBoundary0External : public dgBoundaryCondition {
- boundaryType type() const {
- return EXTERNAL_VALUES;
+ dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
+ return new riemann(_vFunctionName,cacheMapLeft, cacheMapRight);
}
- void operator () (const dgFace &face, fullMatrix<double> fcx[]) const{
- const fullMatrix<double> &solLeft = face.solutionLeft();
- for(int i=0; i< solLeft.size1(); i++) {
- fcx[0](i,0) = 0;
- }
+ dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const {
+ return 0;
}
-};
-
-class dgConservationLawAdvection : public dgConservationLaw {
- public:
- dgConservationLawAdvection()
+ dgConservationLawAdvection(std::string vFunctionName)
{
+ _vFunctionName = vFunctionName;
_nbf = 1;
- _convective = new dgAdvectionFluxTerm;
- _riemannSolver = new dgAdvectionRiemanTerm;
- _boundaryConditions["Left"] = new dgAdvectionBoundary0External();
- _boundaryConditions["Right"] = new dgAdvectionBoundary0External();
- _boundaryConditions["Top"] = new dgAdvectionBoundary0External();
- _boundaryConditions["Bottom"] = new dgAdvectionBoundary0External();
- }
- ~dgConservationLawAdvection()
- {
- delete _convective;
- delete _riemannSolver;
- delete _boundaryConditions["Left"];
- delete _boundaryConditions["Right"];
- delete _boundaryConditions["Top"];
- delete _boundaryConditions["Bottom"];
}
};
-dgConservationLaw *dgNewConservationLawAdvection() {
- return new dgConservationLawAdvection;
+dgConservationLaw *dgNewConservationLawAdvection(std::string vFunctionName) {
+ return new dgConservationLawAdvection(vFunctionName);
}
diff --git a/Solver/dgMainTest.cpp b/Solver/dgMainTest.cpp
index 3d75e92b0e62eb7bf32d1ed670e5fe02e3ef3667..8d085e0e6dc502a75b20fdda42437e8c8f28b183 100644
--- a/Solver/dgMainTest.cpp
+++ b/Solver/dgMainTest.cpp
@@ -44,7 +44,7 @@ int main(int argc, char **argv){
int order=1;
int dimension=2;
dgAlgorithm algo;
- function::registerAllFunctions();
+ function::registerDefaultFunctions();
algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
//for now, we suppose there is only one group of elements
@@ -58,8 +58,21 @@ int main(int argc, char **argv){
algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
}
print("init.pos",*elementGroups[0],&sol(0,0));
+
//advection
- dgConservationLaw *law = dgNewConservationLawAdvection();
+ fullMatrix<double> advectionSpeed(3,1);
+ advectionSpeed(0,0)=0.15;
+ advectionSpeed(1,0)=0.05;
+ advectionSpeed(2,0)=0.;
+
+ function::add("advectionSpeed",function::newFunctionConstant(advectionSpeed));
+
+ dgConservationLaw *law = dgNewConservationLawAdvection("advectionSpeed");
+ law->addBoundaryCondition("Left",dgBoundaryCondition::new0OutCondition(*law));
+ law->addBoundaryCondition("Right",dgBoundaryCondition::new0OutCondition(*law));
+ law->addBoundaryCondition("Top",dgBoundaryCondition::new0FluxCondition(*law));
+ law->addBoundaryCondition("Bottom",dgBoundaryCondition::new0FluxCondition(*law));
+
for(int iT=0; iT<1000; iT++) {
algo.residual(*law,elementGroups,faceGroups,boundaryGroups,sol,residu);
residu.scale(0.01);
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 977ed647982820a215a9f79ab8c2a34fd9c409f4..d83696e91dd1df28162d31ca9522f4e9e4069460 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -77,6 +77,7 @@ dataCacheMap::~dataCacheMap()
// now some example of functions
+// get XYZ coordinates
class functionXYZ : public function {
private:
class data : public dataCacheDouble{
@@ -107,8 +108,41 @@ class functionXYZ : public function {
}
};
+// constant values copied over each line
+class functionConstant : public function {
+ private :
+ class data : public dataCacheDouble {
+ const functionConstant *_function;
+ dataCacheDouble &_uvw;
+ public:
+ data(const functionConstant * function,dataCacheMap *m) :_uvw(m->get("UVW",this)){
+ _function = function;
+ }
+ void _eval() {
+ if(_value.size1()!=_uvw().size1()){
+ _value=fullMatrix<double>(_uvw().size1(),_function->_source.size1());
+ }
+ for(int i=0;i<_value.size1();i++)
+ for(int j=0;j<_function->_source.size1();j++)
+ _value(i,j)=_function->_source(j,0);
+ }
+ };
+ fullMatrix<double> _source;
+ public:
+ dataCacheDouble *newDataCache(dataCacheMap *m)
+ {
+ return new data(this,m);
+ }
+ functionConstant(const fullMatrix<double> &source){
+ _source = source;
+ }
+};
+
+function *function::newFunctionConstant(const fullMatrix<double> &source){
+ return new functionConstant(source);
+}
-void function::registerAllFunctions()
+void function::registerDefaultFunctions()
{
function::add("XYZ", new functionXYZ);
}
diff --git a/Solver/function.h b/Solver/function.h
index bc756ea772cd913e3d8c7ae5869d3965451335be..4f64dfe1fcc1f5ca5f124c505aa8ba03d2685d94 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -28,6 +28,10 @@ protected :
void addMeAsDependencyOf (dataCache *newDep);
dataCache() : _valid(false) {}
virtual ~dataCache(){};
+public :
+ inline bool somethingDependOnMe() {
+ return !_dependOnMe.empty();
+ }
};
// A node in the dependency tree for which all the leafs depend on the
@@ -70,11 +74,14 @@ class function {
static std::map<std::string, function*> _allFunctions;
std::string _name;
public:
- static void registerAllFunctions();
+ static void registerDefaultFunctions();
static bool add(const std::string functionName, function *f);
static function *get(std::string functionName, bool acceptNull=false);
virtual dataCacheDouble *newDataCache(dataCacheMap *m) =0;
+
+ //we need a parser for this
+ static function *newFunctionConstant(const fullMatrix<double> &source);
};
// A special node in the dependency tree for which all the leafs
@@ -97,15 +104,19 @@ class dataCacheMap {
dataCacheElement _cacheElement;
std::map<std::string, dataCacheDouble*> _cacheDoubleMap;
class providedDataDouble : public dataCacheDouble
+ // for data provided by the algorithm and that does not have an _eval function
+ // (typically "UVW") this class is not stricly necessary, we could write
+ // a function for each case but I think it's more practical like this
{
void _eval() {throw;};
+ public:
+ providedDataDouble() {
+ _valid=true;
+ }
};
public:
dataCacheDouble &get(const std::string &functionName, dataCache *caller=0);
dataCacheElement &getElement(dataCache *caller=0);
- // for data provided by the user and that does not have an _eval function
- // (typically "UVW") this class is not stricly necessary, we could write
- // a function for each case but I think it's more practical like this
dataCacheDouble &provideData(std::string name);
~dataCacheMap();
};