diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 2e8c790735b1e0507e0b9b4f84904253befa8076..b1890cd9ba0720ef1b53456c7a46f82734e923f1 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -45,6 +45,7 @@ class dgConservationLaw {
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;}
+ virtual dataCacheDouble *newBifurcationRiemannSolver (dataCacheMap &cacheMap0, dataCacheMap &cacheMap1, dataCacheMap &cacheMap2) const {return NULL;}
virtual dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const {return NULL;}
inline const dgBoundaryCondition *getBoundaryCondition(const std::string tag) const {
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index f9c34e96527994b17b9a3e9f6a4a6959863a2607..c98b9c8b3f32088eb63f0ef0e5a85c70e36a5a3d 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -196,59 +196,32 @@ dgGroupOfConnections::dgGroupOfConnections(const dgGroupOfElements &elementGroup
_fs = _elementGroup.getElement(0)->getFunctionSpace(pOrder);
}
-void dgGroupOfFaces::mapToInterface ( int nFields,
- const fullMatrix<double> &vLeft,
- const fullMatrix<double> &vRight,
- fullMatrix<double> &v)
+void dgGroupOfFaces::mapToInterface (int nFields, std::vector<const fullMatrix<double> *> &proxies, fullMatrix<double> &v)
{
- if(_connections.size()==1){
- for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureLeft = _connections[0]->getClosure(i);
- for (int iField=0; iField<nFields; iField++){
- for(size_t j =0; j < closureLeft.size(); j++){
- v(j, i*nFields + iField) = vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField);
- }
- }
- }
-
- }else{
- for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureLeft = _connections[0]->getClosure(i);
- const std::vector<int> &closureRight = _connections[1]->getClosure(i);
- for (int iField=0; iField<nFields; iField++){
- for(size_t j =0; j < closureLeft.size(); j++){
- v(j, i*2*nFields + iField) = vLeft(closureLeft[j],_connections[0]->getElementId(i)*nFields + iField);
+ int nbConnections = getNbGroupOfConnections();
+ for (int i=0; i<getNbElements(); i++) {
+ for (int iConnection = 0; iConnection<getNbGroupOfConnections(); iConnection++) {
+ const std::vector<int> &closure = _connections[iConnection]->getClosure(i);
+ for (int iField=0; iField<nFields; iField++) {
+ for(size_t j =0; j < closure.size(); j++) {
+ v(j, (i*nbConnections+iConnection)*nFields + iField) =
+ (*proxies[iConnection])(closure[j],_connections[iConnection]->getElementId(i)*nFields + iField);
}
- for(size_t j =0; j < closureRight.size(); j++)
- v(j, (i*2+1)*nFields + iField) = vRight(closureRight[j],_connections[1]->getElementId(i)*nFields + iField);
}
}
}
}
-void dgGroupOfFaces::mapFromInterface ( int nFields,
- const fullMatrix<double> &v,
- fullMatrix<double> &vLeft,
- fullMatrix<double> &vRight
- )
+void dgGroupOfFaces::mapFromInterface (int nFields, const fullMatrix<double> &v, std::vector< fullMatrix<double> *> &proxies)
{
- if(_connections.size()==1){
- for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureLeft = _connections[0]->getClosure(i);
- for (int iField=0; iField<nFields; iField++){
- for(size_t j =0; j < closureLeft.size(); j++)
- vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField) += v(j, i*nFields + iField);
- }
- }
- }else{
- for(int i=0; i<getNbElements(); i++) {
- const std::vector<int> &closureLeft = _connections[0]->getClosure(i);
- const std::vector<int> &closureRight = _connections[1]->getClosure(i);
+ int nbConnections = getNbGroupOfConnections();
+ for (int i=0; i<getNbElements(); i++) {
+ for (int iConnection = 0; iConnection<getNbGroupOfConnections(); iConnection++) {
+ const std::vector<int> &closure = _connections[iConnection]->getClosure(i);
for (int iField=0; iField<nFields; iField++){
- for(size_t j =0; j < closureLeft.size(); j++)
- vLeft(closureLeft[j], _connections[0]->getElementId(i)*nFields + iField) += v(j, i*2*nFields + iField);
- for(size_t j =0; j < closureRight.size(); j++)
- vRight(closureRight[j], _connections[1]->getElementId(i)*nFields + iField) += v(j, (i*2+1)*nFields + iField);
+ for(size_t j =0; j < closure.size(); j++)
+ (*proxies[iConnection])(closure[j], _connections[iConnection]->getElementId(i)*nFields + iField) +=
+ v(j, (i*nbConnections+iConnection)*nFields + iField);
}
}
}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 42c8e1c218cc47d496275b340089284f242f3201..0791de9c301ced2193e6e49661ba93106629c012 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -209,8 +209,8 @@ public:
const polynomialBasis * getPolynomialBasis() const {return _fsFace;}
inline MElement* getFace (int iElement) const {return _faces[iElement];}
public: // should be more generic on the number of connections
- void mapToInterface(int nFields, const fullMatrix<double> &vLeft, const fullMatrix<double> &vRight, fullMatrix<double> &v);
- void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight);
+ void mapToInterface (int nFields, std::vector<const fullMatrix<double> *> &proxies, fullMatrix<double> &v);
+ void mapFromInterface (int nFields, const fullMatrix<double> &v, std::vector<fullMatrix<double> *> &proxies);
void mapLeftFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft);
void mapRightFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vRight);
};
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index 02cecb65bd58f96f0ac6096f4413165c01733efb..26d825742b1f1d21edd241e1146ad00fd903b49f 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -261,134 +261,127 @@ void dgResidualVolume::computeAndMap1Group(dgGroupOfElements &group, dgDofContai
}
-dgResidualInterface::dgResidualInterface (const dgConservationLaw &claw) :
- _cacheMapLeft (new dataCacheMap),
- _cacheMapRight (new dataCacheMap),
- _claw(claw),
- _nbFields(claw.getNbFields()),
- _cacheElementLeft(_cacheMapLeft->getElement()),
- _cacheElementRight(_cacheMapRight->getElement()),
- _uvwLeft(_cacheMapLeft->provideParametricCoordinates()),
- _uvwRight(_cacheMapRight->provideParametricCoordinates()),
- _solutionQPLeft(_cacheMapLeft->provideSolution(_nbFields)),
- _solutionQPRight(_cacheMapRight->provideSolution(_nbFields)),
- _gradientSolutionLeft(_cacheMapLeft->provideSolutionGradient(_nbFields)),
- _gradientSolutionRight(_cacheMapRight->provideSolutionGradient(_nbFields)),
- _normals(_cacheMapLeft->provideNormals()),
- _riemannSolver(claw.newRiemannSolver(*_cacheMapLeft,*_cacheMapRight)),
- _diffusiveFluxLeft(claw.newDiffusiveFlux(*_cacheMapLeft)),
- _diffusiveFluxRight(claw.newDiffusiveFlux(*_cacheMapRight)),
- _maximumDiffusivityLeft(claw.newMaximumDiffusivity(*_cacheMapLeft)),
- _maximumDiffusivityRight(claw.newMaximumDiffusivity(*_cacheMapRight))
+dgResidualInterface::dgResidualInterface (const dgConservationLaw &claw):
+ _claw (claw)
{
+ _nbFields = _claw.getNbFields();
}
void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager (maybe useless here)
dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &solutionLeft,
- fullMatrix<double> &solutionRight,
+ const std::vector<const fullMatrix<double>*> & solutionOnElements,
fullMatrix<double> &residual // residual !! at faces nodes
)
{
// A) global operations before entering the loop over the faces
// A1 ) copy locally some references from the group for the sake of lisibility
- const dgGroupOfConnections &left = group.getGroupOfConnections(0);
- const dgGroupOfConnections &right = group.getGroupOfConnections(1);
- const dgGroupOfElements &groupLeft = left.getGroupOfElements();
- const dgGroupOfElements &groupRight = right.getGroupOfElements();
- int nbFaces = group.getNbElements();
- //get matrices needed to compute the gradient on both sides
- const fullMatrix<double> &DPsiLeftDx = left.getDPsiDx();
- const fullMatrix<double> &DPsiRightDx = right.getDPsiDx();
+ int nbConnections = group.getNbGroupOfConnections();
+ std::vector<dataCacheMap> caches(nbConnections);
+ std::vector<const dgGroupOfConnections*> connections(nbConnections);
+ std::vector<const dgGroupOfElements*> elementGroups(nbConnections);
+ std::vector<dataCacheDouble*> diffusiveFluxes(nbConnections), maximumDiffusivities(nbConnections);
+ for (int i=0; i<nbConnections; i++) {
+ connections[i] = &group.getGroupOfConnections(i);
+ elementGroups[i] = &connections[i]->getGroupOfElements();
+ caches[i].provideParametricCoordinates();
+ caches[i].provideSolution(_nbFields);
+ caches[i].provideSolutionGradient(_nbFields);
+ caches[i].provideNormals();
+ diffusiveFluxes[i] =_claw.newDiffusiveFlux(caches[i]);
+ maximumDiffusivities[i] = _claw.newMaximumDiffusivity(caches[i]);
+ caches[i].setNbEvaluationPoints(group.getNbIntegrationPoints());
+ }
+ dataCacheDouble *riemannSolver = NULL;
+ if (nbConnections == 2)
+ riemannSolver = _claw.newRiemannSolver(caches[0], caches[1]);
+ if (nbConnections == 3)
+ riemannSolver = _claw.newBifurcationRiemannSolver(caches[0], caches[1], caches[2]);
+
// ----- 1 ---- get the solution at quadrature points
- fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),nbFaces * _nbFields*2);
- group.getCollocationMatrix().mult(solution, solutionQP);
+ int nbFaces = group.getNbElements();
+ fullMatrix<double> solutionQP (group.getNbIntegrationPoints(), nbFaces * _nbFields*nbConnections);
+ group.getCollocationMatrix().mult (solution, solutionQP);
// needed tocompute normal fluxes at integration points
- fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), nbFaces*_nbFields*2);
- // create one dataCache for each side
- _cacheMapLeft->setNbEvaluationPoints(group.getNbIntegrationPoints());
- _cacheMapRight->setNbEvaluationPoints(group.getNbIntegrationPoints());
- fullMatrix<double> dPsiLeftDx,dPsiRightDx,dofsLeft,dofsRight,normalFluxQP;
- int p = groupLeft.getOrder();
- int dim = left.getElement(0)->getDim();
+ fullMatrix<double> NormalFluxQP (group.getNbIntegrationPoints(), nbFaces*_nbFields*nbConnections);
+ fullMatrix<double> normalFluxQP, dofs, dPsiDx;
+ int p = elementGroups[0]->getOrder();
+ int dim = connections[0]->getElement(0)->getDim();
for (int iFace=0 ; iFace < nbFaces ; ++iFace) {
- // B1 ) adjust the proxies for this element
- // NB : the B1 section will almost completely disapear
- // if we write a function (from the function class) for moving proxies across big matrices
- // give the current elements to the dataCacheMap
- _cacheElementLeft.set(left.getElement(iFace));
- _cacheElementRight.set(left.getElement(iFace));
- // proxies for the solution
- _solutionQPLeft.setAsProxy(solutionQP, iFace*2*_nbFields, _nbFields);
- _solutionQPRight.setAsProxy(solutionQP, (iFace*2+1)*_nbFields, _nbFields);
- _normals.setAsProxy(left.getNormals(), iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
- // proxies needed to compute the gradient of the solution and the IP term
- dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*groupLeft.getNbNodes(),groupLeft.getNbNodes());
- dPsiRightDx.setAsProxy(DPsiRightDx,iFace*groupRight.getNbNodes(),groupRight.getNbNodes());
- dofsLeft.setAsProxy(solutionLeft, _nbFields*left.getElementId(iFace), _nbFields);
- dofsRight.setAsProxy(solutionRight, _nbFields*right.getElementId(iFace), _nbFields);
- _uvwLeft.setAsProxy(left.getIntegrationPointsOnElement(iFace));
- _uvwRight.setAsProxy(right.getIntegrationPointsOnElement(iFace));
- // proxies for the flux
- normalFluxQP.setAsProxy(NormalFluxQP, iFace*_nbFields*2, _nbFields*2);
- // B2 ) compute the gradient of the solution
- if(_gradientSolutionLeft.somethingDependOnMe()){
- dPsiLeftDx.mult(dofsLeft, _gradientSolutionLeft.set());
- dPsiRightDx.mult(dofsRight, _gradientSolutionRight.set());
+ for (int i=0; i<nbConnections; i++) {
+ // B1 ) adjust the proxies for this element
+ caches[i].getElement().set(connections[i]->getElement(iFace));
+ caches[i].getParametricCoordinates(NULL).setAsProxy(connections[i]->getIntegrationPointsOnElement(iFace));
+ caches[i].getSolution(NULL).setAsProxy(solutionQP, (iFace*nbConnections+i)*_nbFields, _nbFields);
+ caches[i].getNormals(NULL).setAsProxy(connections[i]->getNormals(), iFace*group.getNbIntegrationPoints(), group.getNbIntegrationPoints());
+ // B2 ) compute the gradient of the solution
+ if(caches[i].getSolutionGradient(NULL).somethingDependOnMe()) {
+ dofs.setAsProxy(*solutionOnElements[i], _nbFields*connections[i]->getElementId(iFace), _nbFields);
+ dPsiDx.setAsProxy(connections[i]->getDPsiDx(),iFace*elementGroups[i]->getNbNodes(),elementGroups[i]->getNbNodes());
+ dPsiDx.mult(dofs, caches[i].getSolutionGradient(NULL).set());
+ }
}
+ // proxies for the flux
+ normalFluxQP.setAsProxy(NormalFluxQP, iFace*_nbFields*nbConnections, _nbFields*nbConnections);
// B3 ) compute fluxes
- if (_diffusiveFluxLeft) {
+ if (diffusiveFluxes[0] && nbConnections==2) { // IP only implemented for 2-sides faces (not for bifurcations)
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
//just for the lisibility :
- const fullMatrix<double> &dfl = (*_diffusiveFluxLeft)();
- const fullMatrix<double> &dfr = (*_diffusiveFluxRight)();
+ const fullMatrix<double> &dfl = (*diffusiveFluxes[0])();
+ const fullMatrix<double> &dfr = (*diffusiveFluxes[1])();
+ const fullMatrix<double> &n = (caches[0].getNormals(NULL))();
for (int k=0;k<_nbFields;k++) {
double meanNormalFlux =
- ((dfl(iPt,k+_nbFields*0 )+dfr(iPt,k+_nbFields*0)) * _normals(0,iPt)
- +(dfl(iPt,k+_nbFields*1 )+dfr(iPt,k+_nbFields*1)) * _normals(1,iPt)
- +(dfl(iPt,k+_nbFields*2 )+dfr(iPt,k+_nbFields*2)) * _normals(2,iPt))/2;
- double minl = std::min(groupLeft.getElementVolume(left.getElementId(iFace)),
- groupRight.getElementVolume(right.getElementId(iFace))
+ ((dfl(iPt,k+_nbFields*0)+dfr(iPt,k+_nbFields*0)) * n(0,iPt)
+ +(dfl(iPt,k+_nbFields*1)+dfr(iPt,k+_nbFields*1)) * n(1,iPt)
+ +(dfl(iPt,k+_nbFields*2)+dfr(iPt,k+_nbFields*2)) * n(2,iPt))/2;
+ double minl = std::min(elementGroups[0]->getElementVolume(connections[0]->getElementId(iFace)),
+ elementGroups[1]->getElementVolume(connections[1]->getElementId(iFace))
)/group.getInterfaceSurface(iFace);
- double nu = std::max((*_maximumDiffusivityRight)(iPt,0),(*_maximumDiffusivityLeft)(iPt,0));
+ double nu = std::max((*maximumDiffusivities[1])(iPt,0),(*maximumDiffusivities[0])(iPt,0));
double mu = (p+1)*(p+dim)/dim*nu/minl;
- double solutionJumpPenalty = (_solutionQPLeft(iPt,k)-_solutionQPRight(iPt,k))/2*mu;
+ double solutionJumpPenalty = (caches[0].getSolution(NULL)(iPt,k)-caches[1].getSolution(NULL)(iPt,k))/2*mu;
normalFluxQP(iPt,k) -= (meanNormalFlux+solutionJumpPenalty)*detJ;
normalFluxQP(iPt,k+_nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ;
}
}
}
- if (_riemannSolver) {
+ if (riemannSolver) {
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) += (*_riemannSolver)(iPt,k)*detJ;
+ for (int k=0;k<_nbFields*nbConnections;k++)
+ normalFluxQP(iPt,k) += (*riemannSolver)(iPt,k)*detJ;
}
}
}
// C ) redistribute the flux to the residual (at Faces nodes)
- if(_riemannSolver || _diffusiveFluxLeft)
+ if(riemannSolver || diffusiveFluxes[0])
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
}
void dgResidualInterface::computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
{
- fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*_nbFields);
- fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*_nbFields);
- const dgGroupOfElements *groupLeft = &faces.getGroupOfConnections(0).getGroupOfElements();
- const dgGroupOfElements *groupRight = &faces.getGroupOfConnections(1).getGroupOfElements();
- faces.mapToInterface(_nbFields, solution.getGroupProxy(groupLeft), solution.getGroupProxy(groupRight), solInterface);
- compute1Group(faces,solInterface,solution.getGroupProxy(groupLeft), solution.getGroupProxy(groupRight),residuInterface);
- faces.mapFromInterface(_nbFields, residuInterface,residual.getGroupProxy(groupLeft), residual.getGroupProxy(groupRight));
+ int nbConnections = faces.getNbGroupOfConnections();
+ fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*_nbFields*nbConnections);
+ fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*_nbFields*nbConnections);
+
+ std::vector<const dgGroupOfElements *> elGroups(nbConnections);
+ std::vector<const fullMatrix<double> *> solutionProxies(nbConnections);
+ std::vector<fullMatrix<double> *> residualProxies(nbConnections);
+ for (int i=0; i<nbConnections; i++) {
+ elGroups[i] = &faces.getGroupOfConnections(i).getGroupOfElements();
+ solutionProxies[i] = &solution.getGroupProxy(elGroups[i]);
+ residualProxies[i] = &residual.getGroupProxy(elGroups[i]);
+ }
+ faces.mapToInterface(_nbFields, solutionProxies, solInterface);
+ compute1Group(faces, solInterface, solutionProxies, residuInterface);
+ faces.mapFromInterface(_nbFields, residuInterface, residualProxies);
}
dgResidualInterface::~dgResidualInterface ()
{
- delete _cacheMapLeft;
- delete _cacheMapRight;
}
@@ -398,7 +391,7 @@ dgResidualBoundary::dgResidualBoundary(const dgConservationLaw &claw):_claw(claw
void dgResidualBoundary::compute1Group(
dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &solutionLeft,
+ const fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
)
{
@@ -486,12 +479,21 @@ void dgResidualBoundary::compute1Group(
void dgResidualBoundary::computeAndMap1Group(dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
{
int _nbFields = _claw.getNbFields();
- fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*_nbFields);
- fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*_nbFields);
- const dgGroupOfElements &groupLeft = faces.getGroupOfConnections(0).getGroupOfElements();
- faces.mapToInterface(_nbFields, solution.getGroupProxy(&groupLeft), solution.getGroupProxy(&groupLeft), solInterface);
- compute1Group(faces,solInterface,solution.getGroupProxy(&groupLeft),residuInterface);
- faces.mapFromInterface(_nbFields, residuInterface, residual.getGroupProxy(&groupLeft), residual.getGroupProxy(&groupLeft));
+ int nbConnections = faces.getNbGroupOfConnections();
+ fullMatrix<double> solInterface(faces.getNbNodes(), faces.getNbElements()*_nbFields*nbConnections);
+ fullMatrix<double> residuInterface(faces.getNbNodes(), faces.getNbElements()*_nbFields*nbConnections);
+
+ std::vector<const dgGroupOfElements *> elGroups(nbConnections);
+ std::vector<const fullMatrix<double> *> solutionProxies(nbConnections);
+ std::vector<fullMatrix<double> *> residualProxies(nbConnections);
+ for (int i=0; i<nbConnections; i++) {
+ elGroups[i] = &faces.getGroupOfConnections(i).getGroupOfElements();
+ solutionProxies[i] = &solution.getGroupProxy(elGroups[i]);
+ residualProxies[i] = &residual.getGroupProxy(elGroups[i]);
+ }
+ faces.mapToInterface(_nbFields, solutionProxies, solInterface);
+ compute1Group(faces, solInterface, *solutionProxies[0],residuInterface);
+ faces.mapFromInterface(_nbFields, residuInterface, residualProxies);
}
void dgResidual::compute(dgGroupCollection &groups, dgDofContainer &solution, dgDofContainer &residual)
@@ -540,9 +542,9 @@ void dgResidual::registerBindings (binding *b)
mb = cb->addMethod("computeAndMap1Group",&dgResidualInterface::computeAndMap1Group);
mb->setDescription("compute the residual for one group given a dgDofContainer solution");
mb->setArgNames("group", "solution", "residual", NULL);
- mb = cb->addMethod("compute1Group", &dgResidualInterface::compute1Group);
- mb->setDescription("compute the residual for one group given fullMatrices with the solution at faces nodes and at element nodes");
- mb->setArgNames("group", "solutionFaces", "solutionGroupLeft", "solutionGroupRight", "residual", NULL);
+// mb = cb->addMethod("compute1Group", &dgResidualInterface::compute1Group);
+ //mb->setDescription("compute the residual for one group given fullMatrices with the solution at faces nodes and at element nodes");
+ //mb->setArgNames("group", "solutionFaces", "solutionOnElements", "residual", NULL);
cb = b->addClass<dgResidualBoundary>("dgResidualBoundary");
cb->setDescription("compute the residual for one boundary dgGroupOfFaces");
diff --git a/Solver/dgResidual.h b/Solver/dgResidual.h
index 47ab9226d2e93fb4f83377dd47c71d67754a8aa2..3b59fecfb9a7decd9a08727fc555c0ea48a7011f 100644
--- a/Solver/dgResidual.h
+++ b/Solver/dgResidual.h
@@ -10,6 +10,7 @@ class dgGroupOfElements;
class dgGroupOfFaces;
class binding;
#include "fullMatrix.h"
+#include <vector>
class dgResidualVolume {
dataCacheMap *_cacheMap;
@@ -27,20 +28,14 @@ class dgResidualVolume {
};
class dgResidualInterface {
- dataCacheMap *_cacheMapLeft, *_cacheMapRight;
const dgConservationLaw &_claw;
int _nbFields;
- dataCacheElement &_cacheElementLeft, &_cacheElementRight;
- dataCacheDouble &_uvwLeft, &_uvwRight, &_solutionQPLeft, &_solutionQPRight, &_gradientSolutionLeft, &_gradientSolutionRight;
- dataCacheDouble &_normals;
- dataCacheDouble *_riemannSolver, *_maximumDiffusivityLeft,*_maximumDiffusivityRight, *_diffusiveFluxLeft, *_diffusiveFluxRight;
public:
dgResidualInterface (const dgConservationLaw &claw);
void compute1Group ( //dofManager &dof, // the DOF manager (maybe useless here)
dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &solutionLeft,
- fullMatrix<double> &solutionRight,
+ const std::vector<const fullMatrix<double>*> & solutionOnElements,
fullMatrix<double> &residual // residual !! at faces nodes
);
void computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual);
@@ -53,7 +48,7 @@ class dgResidualBoundary {
void compute1Group ( //dofManager &dof, // the DOF manager (maybe useless here)
dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &solutionLeft,
+ const fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
);
void computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual);
diff --git a/Solver/dgRungeKuttaMultirate.cpp b/Solver/dgRungeKuttaMultirate.cpp
index 9f8238e04418894ee554c1a0226208fcd75e3610..23f41da85b4186e6d29cbc602df10f5c2c5143e3 100644
--- a/Solver/dgRungeKuttaMultirate.cpp
+++ b/Solver/dgRungeKuttaMultirate.cpp
@@ -198,8 +198,11 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
const dgGroupOfElements *gR = &faces->getGroupOfConnections(1).getGroupOfElements();
fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
- faces->mapToInterface(_law->getNbFields(), _currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR), solInterface);
- _residualInterface->compute1Group(*faces,solInterface,_currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR),residuInterface);
+ std::vector<const fullMatrix<double>*> proxies(2);
+ proxies[0] = &_currentInput->getGroupProxy(gL);
+ proxies[1] = &_currentInput->getGroupProxy(gR);
+ faces->mapToInterface(_law->getNbFields(), proxies, solInterface);
+ _residualInterface->compute1Group(*faces,solInterface, proxies, residuInterface);
if(gL->getMultirateExponent()==exponent && gL->getIsMultirateBuffer()){
faces->mapLeftFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gL));
}
@@ -247,8 +250,11 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
const dgGroupOfElements *gR = &faces->getGroupOfConnections(1).getGroupOfElements();
fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
- faces->mapToInterface(_law->getNbFields(), _currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR), solInterface);
- _residualInterface->compute1Group(*faces,solInterface,_currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR),residuInterface);
+ std::vector<const fullMatrix<double>*> proxies(2);
+ proxies[0] = &_currentInput->getGroupProxy(gL);
+ proxies[1] = &_currentInput->getGroupProxy(gR);
+ faces->mapToInterface(_law->getNbFields(), proxies, solInterface);
+ _residualInterface->compute1Group(*faces,solInterface, proxies, residuInterface);
if(gL->getMultirateExponent()==exponent && !gL->getIsMultirateBuffer()){
faces->mapLeftFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gL));
}
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 6472f00a7abaaab611cf124e86128dfe17227037..db348626dd7c8118a2085453c7cf5c06079c0506 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -58,8 +58,8 @@ function *function::get(std::string functionName, bool acceptNull)
dataCacheElement &dataCacheMap::getElement(dataCache *caller)
{
if(caller)
- _cacheElement->addMeAsDependencyOf(caller);
- return *_cacheElement;
+ _cacheElement.addMeAsDependencyOf(caller);
+ return _cacheElement;
}
static dataCacheDouble &returnDataCacheDouble(dataCacheDouble *data, dataCache *caller)
diff --git a/Solver/function.h b/Solver/function.h
index 55ba50b13778ea5ee368df31d69a5c4c43389f13..74764b99cb53cc2fdce355a4fe9845cfe14478b5 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -164,7 +164,6 @@ class dataCacheMap {
int _nbEvaluationPoints;
// keep track of the current element and all the dataCaches that
// depend on it
- dataCacheElement *_cacheElement;
std::map<const function*, dataCacheDouble*> _cacheDoubleMap;
class providedDataDouble : public dataCacheDouble
// for data provided by the algorithm and that does not have an _eval function
@@ -179,6 +178,7 @@ class dataCacheMap {
std::set<dataCache*> _toDelete;
std::set<dataCacheDouble*> _toResize;
+ dataCacheElement _cacheElement;
protected:
void addDataCache(dataCache *data){
_toDelete.insert(data);
@@ -205,8 +205,8 @@ class dataCacheMap {
dataCacheDouble &get(const function *f, dataCache *caller=0);
dataCacheElement &getElement(dataCache *caller=0);
- dataCacheMap(){
- _cacheElement = new dataCacheElement(this);
+ dataCacheMap():_cacheElement(this){
+ _toDelete.erase(&_cacheElement);
_normals = _solution = _solutionGradient = _parametricCoordinates = 0;
_nbEvaluationPoints = 0;
}