From e9522f0139e4c4dcb26e7ba61910652a358bd485 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Wed, 17 Mar 2010 14:29:51 +0000
Subject: [PATCH] dg : merge residualBoundary into residualInterface : same
code for boundaries, regular 2-sided-faces and bifurcations
---
Solver/TESTCASES/AdvectionDiffusion.lua | 2 +-
Solver/TESTCASES/Diffusion.lua | 17 +--
Solver/TESTCASES/supra.lua | 12 +-
Solver/dgGroupOfElements.cpp | 5 -
Solver/dgGroupOfElements.h | 2 -
Solver/dgLimiter.cpp | 2 +
Solver/dgResidual.cpp | 179 ++++++------------------
Solver/dgResidual.h | 13 --
Solver/dgRungeKuttaMultirate.cpp | 37 ++---
Solver/dgRungeKuttaMultirate.h | 1 -
10 files changed, 63 insertions(+), 207 deletions(-)
diff --git a/Solver/TESTCASES/AdvectionDiffusion.lua b/Solver/TESTCASES/AdvectionDiffusion.lua
index a27a315052..981314d0cd 100644
--- a/Solver/TESTCASES/AdvectionDiffusion.lua
+++ b/Solver/TESTCASES/AdvectionDiffusion.lua
@@ -3,7 +3,7 @@ model:load ('square.geo')
model:load ('square.msh')
dg = dgSystemOfEquations (model)
-dg:setOrder(4)
+dg:setOrder(1)
-- conservation law
law = dgConservationLawAdvectionDiffusion(functionConstant({0.15,0.05,0}),functionConstant({0.001}))
diff --git a/Solver/TESTCASES/Diffusion.lua b/Solver/TESTCASES/Diffusion.lua
index 2b2800dcf5..c826589a6d 100644
--- a/Solver/TESTCASES/Diffusion.lua
+++ b/Solver/TESTCASES/Diffusion.lua
@@ -7,15 +7,12 @@ dg:setOrder(1)
-- conservation law
-- advection speed
-nu=fullMatrix(1,1);
-nu:set(0,0,0.01)
-law = dgConservationLawAdvectionDiffusion('',functionConstant(nu):getName())
+law = dgConservationLawAdvectionDiffusion.diffusionLaw(functionConstant({0.01}))
dg:setConservationLaw(law)
-- boundary condition
-outside=fullMatrix(1,1)
-outside:set(0,0,0.)
-law:addBoundaryCondition('Border',law:new0FluxBoundary())
+law:addBoundaryCondition('Border',law:newOutsideValueBoundary(functionConstant({1})))
+--law:addBoundaryCondition('Border',law:new0FluxBoundary())
dg:setup()
@@ -28,14 +25,14 @@ 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/Diffusion_00000')
-- main loop
-for i=1,10000 do
- --norm = dg:RK44(0.001)
- norm = dg:RK44_TransformNodalValue(0.00001)
+for i=1,100 do
+ norm = dg:RK44(0.01)
+ --norm = dg:RK44_TransformNodalValue(0.00001)
if (i % 50 == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/Diffusion-%05d", i))
diff --git a/Solver/TESTCASES/supra.lua b/Solver/TESTCASES/supra.lua
index 9187b68ad3..1e828128ae 100644
--- a/Solver/TESTCASES/supra.lua
+++ b/Solver/TESTCASES/supra.lua
@@ -28,10 +28,8 @@ function valueLeft(f)
end
end
-law = dgConservationLawAdvectionDiffusion('',functionLua(1,'diffusivity',{'Solution'}):getName())
-print(law)
+law = dgConservationLawAdvectionDiffusion.diffusionLaw(functionLua(1,'diffusivity',{functionSolution.get()}))
dg:setConservationLaw(law)
-print(law)
-- boundary condition
--[[
@@ -42,8 +40,8 @@ outsideRight:set(0,0,1)
--]]
law:addBoundaryCondition('Top',law:new0FluxBoundary())
law:addBoundaryCondition('Bottom',law:new0FluxBoundary())
-leftFunction=functionLua(1,'valueLeft',{}):getName()
-rightFunction=functionLua(1,'valueRight',{}):getName()
+leftFunction=functionLua(1,'valueLeft',{})
+rightFunction=functionLua(1,'valueRight',{})
--[[
law:addBoundaryCondition('Left',law:newNeumannBoundary(leftFunction))
law:addBoundaryCondition('Right',law:newNeumannBoundary(rightFunction))
@@ -62,7 +60,7 @@ function initial_condition( xyz , f )
f:set (i, 0,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/supra-00000')
@@ -70,7 +68,7 @@ dg:exportSolution('output/supra-00000')
for i=1,150000 do
norm = dg:RK44_limiter(dt)
t=t+dt
- if (i % 100 == 0) then
+ if (i % 1 == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/supra-%05d", i))
end
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index c98b9c8b3f..11e57af596 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -1236,15 +1236,10 @@ void dgGroupCollection::registerBindings(binding *b)
cm->setDescription("return the number of dgGroupOfElements");
cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups);
cm->setDescription("return the number of dgGroupOfFaces (interior ones, not the domain boundaries)");
- cm = cb->addMethod("getNbBoundaryGroups", &dgGroupCollection::getNbBoundaryGroups);
- cm->setDescription("return the number of group of boundary faces");
cm = cb->addMethod("getElementGroup", &dgGroupCollection::getElementGroup);
cm->setDescription("get 1 group of elements");
cm->setArgNames("id", NULL);
cm = cb->addMethod("getFaceGroup", &dgGroupCollection::getFaceGroup);
cm->setDescription("get 1 group of faces");
cm->setArgNames("id", NULL);
- cm = cb->addMethod("getBoundaryGroup", &dgGroupCollection::getBoundaryGroup);
- cm->setDescription("get 1 group of boundary faces");
- cm->setArgNames("id", NULL);
}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 0791de9c30..9c3e5a2506 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -236,11 +236,9 @@ class dgGroupCollection {
inline GModel* getModel() {return _model;}
inline int getNbElementGroups() const {return _elementGroups.size();}
inline int getNbFaceGroups() const {return _faceGroups.size();}
- inline int getNbBoundaryGroups() const {return _boundaryGroups.size();}
inline int getNbGhostGroups() const {return _ghostGroups.size();}
inline dgGroupOfElements *getElementGroup(int i) const {return i<getNbElementGroups()?_elementGroups[i]:_ghostGroups[i-getNbElementGroups()];}
inline dgGroupOfFaces *getFaceGroup(int i) const {return _faceGroups[i];}
- inline dgGroupOfFaces *getBoundaryGroup(int i) const {return _boundaryGroups[i];}
inline dgGroupOfElements *getGhostGroup(int i) const {return _ghostGroups[i];}
inline int getNbImageElementsOnPartition(int partId) const {return _elementsToSend[partId].size();}
diff --git a/Solver/dgLimiter.cpp b/Solver/dgLimiter.cpp
index 7e36cae0bf..dd8286ca40 100644
--- a/Solver/dgLimiter.cpp
+++ b/Solver/dgLimiter.cpp
@@ -22,6 +22,8 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution)
fullMatrix<double> TempL, TempR;
for( int iGFace=0; iGFace<groups->getNbFaceGroups(); iGFace++) {
dgGroupOfFaces* group = groups->getFaceGroup(iGFace);
+ if (group->getNbGroupOfConnections()!=2) // skip boundaries and bifurcations
+ continue;
const dgGroupOfConnections &left = group->getGroupOfConnections(0);
const dgGroupOfConnections &right = group->getGroupOfConnections(1);
const dgGroupOfElements *groupLeft = &left.getGroupOfElements();
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index 26d825742b..c872c42ea3 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -292,7 +292,16 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
maximumDiffusivities[i] = _claw.newMaximumDiffusivity(caches[i]);
caches[i].setNbEvaluationPoints(group.getNbIntegrationPoints());
}
- dataCacheDouble *riemannSolver = NULL;
+ dataCacheDouble *riemannSolver = NULL, *neumann = NULL, *dirichlet = NULL, *maximumDiffusivityOut = NULL;
+ if (nbConnections == 1) {
+ const dgBoundaryCondition *boundaryCondition = _claw.getBoundaryCondition(group.getPhysicalTag());
+ riemannSolver = boundaryCondition->newBoundaryTerm(caches[0]);
+ neumann = boundaryCondition->newDiffusiveNeumannBC(caches[0]);
+ dirichlet = boundaryCondition->newDiffusiveDirichletBC(caches[0]);
+ maximumDiffusivityOut = boundaryCondition->newMaximumDiffusivity(caches[0]);
+ if (!maximumDiffusivityOut)
+ maximumDiffusivityOut = maximumDiffusivities[0];
+ }
if (nbConnections == 2)
riemannSolver = _claw.newRiemannSolver(caches[0], caches[1]);
if (nbConnections == 3)
@@ -324,26 +333,39 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
// proxies for the flux
normalFluxQP.setAsProxy(NormalFluxQP, iFace*_nbFields*nbConnections, _nbFields*nbConnections);
// B3 ) compute fluxes
- 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);
+ if (diffusiveFluxes[0]) { // IP only implemented for 2-sides faces (not for bifurcations)
+ if(nbConnections == 2) {
//just for the lisibility :
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 =
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iFace, iPt);
+ for (int k=0;k<_nbFields;k++) {
+ double meanNormalFlux =
((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((*maximumDiffusivities[1])(iPt,0),(*maximumDiffusivities[0])(iPt,0));
- double mu = (p+1)*(p+dim)/dim*nu/minl;
- 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;
+ +(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((*maximumDiffusivities[1])(iPt,0),(*maximumDiffusivities[0])(iPt,0));
+ double mu = (p+1)*(p+dim)/dim*nu/minl;
+ 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;
+ }
+ }
+ } else if(nbConnections == 1) {
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iFace, iPt);
+ for (int k=0;k<_nbFields;k++) {
+ double minl = elementGroups[0]->getElementVolume(connections[0]->getElementId(iFace))/group.getInterfaceSurface(iFace);
+ double nu = std::max((*maximumDiffusivities[0])(iPt,0),(*maximumDiffusivityOut)(iPt,0));
+ double mu = (p+1)*(p+dim)/dim*nu/minl;
+ double solutionJumpPenalty = (caches[0].getSolution(NULL)(iPt,k)-(*dirichlet)(iPt,k))/2*mu;
+ normalFluxQP(iPt,k) -= ((*neumann)(iPt,k)+solutionJumpPenalty)*detJ;
+ }
}
}
}
@@ -354,6 +376,7 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
normalFluxQP(iPt,k) += (*riemannSolver)(iPt,k)*detJ;
}
}
+
}
// C ) redistribute the flux to the residual (at Faces nodes)
if(riemannSolver || diffusiveFluxes[0])
@@ -384,118 +407,6 @@ dgResidualInterface::~dgResidualInterface ()
{
}
-
-dgResidualBoundary::dgResidualBoundary(const dgConservationLaw &claw):_claw(claw){
-}
-
-void dgResidualBoundary::compute1Group(
- dgGroupOfFaces &group,
- const fullMatrix<double> &solution, // solution !! at faces nodes
- const fullMatrix<double> &solutionLeft,
- fullMatrix<double> &residual // residual !! at faces nodes
- )
-{
- //should be splitted like dgResidualInterface and dgResidualVolume
- //but i do not do it know because dgResidualBoundary will probably disapear when we will have list of elements on interfaces
- const dgGroupOfConnections &left = group.getGroupOfConnections(0);
- const dgGroupOfElements &groupLeft = left.getGroupOfElements();
- int nbFields = _claw.getNbFields();
- const dgBoundaryCondition *boundaryCondition = _claw.getBoundaryCondition(group.getPhysicalTag());
- // ----- 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);
-
- dataCacheMap cacheMapLeft;
- cacheMapLeft.setNbEvaluationPoints(group.getNbIntegrationPoints());
- const fullMatrix<double> &DPsiLeftDx = left.getDPsiDx();
- // provided dataCache
- dataCacheDouble &uvw=cacheMapLeft.provideParametricCoordinates();
- dataCacheDouble &solutionQPLeft = cacheMapLeft.provideSolution(nbFields);
- dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideSolutionGradient(nbFields);
- dataCacheDouble &normals = cacheMapLeft.provideNormals();
- dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
- // required data
- // inviscid
- dataCacheDouble *boundaryTerm = boundaryCondition->newBoundaryTerm(cacheMapLeft);
- // viscous
- dataCacheDouble *diffusiveFluxLeft = _claw.newDiffusiveFlux(cacheMapLeft);
- dataCacheDouble *maximumDiffusivityLeft = _claw.newMaximumDiffusivity(cacheMapLeft);
- dataCacheDouble *maximumDiffusivityOut = boundaryCondition->newMaximumDiffusivity(cacheMapLeft);
- dataCacheDouble *neumann = boundaryCondition->newDiffusiveNeumannBC(cacheMapLeft);
- dataCacheDouble *dirichlet = boundaryCondition->newDiffusiveDirichletBC(cacheMapLeft);
-
- fullMatrix<double> normalFluxQP,dPsiLeftDx,dofsLeft;
-
- int p = groupLeft.getOrder();
- int dim = left.getElement(0)->getDim();
-
- for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
- normalFluxQP.setAsProxy(NormalFluxQP, iFace*nbFields, nbFields);
- // ----- 2.3.1 --- provide the data to the cacheMap
- solutionQPLeft.setAsProxy(solutionQP, iFace*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());
- dofsLeft.setAsProxy(solutionLeft, nbFields*left.getElementId(iFace), nbFields);
-
- uvw.setAsProxy(left.getIntegrationPointsOnElement(iFace));
- cacheElementLeft.set(left.getElement(iFace));
-
- // compute the gradient of the solution
- if(gradientSolutionLeft.somethingDependOnMe()){
- dPsiLeftDx.mult(dofsLeft, gradientSolutionLeft.set());
- }
-
- // ----- 2.3.2 --- compute inviscid contribution
- 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) = (*boundaryTerm)(iPt,k)*detJ;
- }
-
- // ----- 2.3.3 --- compute viscous contribution
- if (diffusiveFluxLeft) {
- for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
- const double detJ = group.getDetJ (iFace, iPt);
- //just for the lisibility :
- for (int k=0;k<nbFields;k++) {
- double minl = groupLeft.getElementVolume(left.getElementId(iFace))/group.getInterfaceSurface(iFace);
- double nu = (*maximumDiffusivityLeft)(iPt,0);
- if(maximumDiffusivityOut)
- nu = std::max(nu,(*maximumDiffusivityOut)(iPt,0));
- double mu = (p+1)*(p+dim)/dim*nu/minl;
- double solutionJumpPenalty = (solutionQPLeft(iPt,k)-(*dirichlet)(iPt,k))/2*mu;
- normalFluxQP(iPt,k) -= ((*neumann)(iPt,k)+solutionJumpPenalty)*detJ;
- }
- }
- }
- }
- // ----- 3 ---- do the redistribution at face nodes using BLAS3
- residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
-}
-
-void dgResidualBoundary::computeAndMap1Group(dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual)
-{
- int _nbFields = _claw.getNbFields();
- 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)
{
solution.scatter();
@@ -506,9 +417,6 @@ void dgResidual::compute(dgGroupCollection &groups, dgDofContainer &solution, dg
for(size_t i=0;i<groups.getNbFaceGroups() ; i++){
residualInterface.computeAndMap1Group(*groups.getFaceGroup(i), solution, residual);
}
- dgResidualBoundary residualBoundary(_claw);
- for(size_t i=0;i<groups.getNbBoundaryGroups() ; i++)
- residualBoundary.computeAndMap1Group(*groups.getBoundaryGroup(i), solution, residual);
}
#include "Bindings.h"
@@ -545,13 +453,4 @@ void dgResidual::registerBindings (binding *b)
// 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");
- mb = cb->addMethod("computeAndMap1Group",&dgResidualBoundary::computeAndMap1Group);
- mb->setDescription("compute the residual for one group given a dgDofContainer solution");
- mb->setArgNames("group", "solution", "residual", NULL);
- mb = cb->addMethod("compute1Group", &dgResidualBoundary::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", "residual", NULL);
}
diff --git a/Solver/dgResidual.h b/Solver/dgResidual.h
index 3b59fecfb9..e633056eb8 100644
--- a/Solver/dgResidual.h
+++ b/Solver/dgResidual.h
@@ -42,19 +42,6 @@ class dgResidualInterface {
~dgResidualInterface();
};
-class dgResidualBoundary {
- const dgConservationLaw &_claw;
- public :
- void compute1Group ( //dofManager &dof, // the DOF manager (maybe useless here)
- dgGroupOfFaces &group,
- const fullMatrix<double> &solution, // solution !! at faces nodes
- const fullMatrix<double> &solutionLeft,
- fullMatrix<double> &residual // residual !! at faces nodes
- );
- void computeAndMap1Group (dgGroupOfFaces &faces, dgDofContainer &solution, dgDofContainer &residual);
- dgResidualBoundary (const dgConservationLaw &claw);
-};
-
class dgResidual {
const dgConservationLaw &_claw;
public:
diff --git a/Solver/dgRungeKuttaMultirate.cpp b/Solver/dgRungeKuttaMultirate.cpp
index 23f41da85b..204eb10d70 100644
--- a/Solver/dgRungeKuttaMultirate.cpp
+++ b/Solver/dgRungeKuttaMultirate.cpp
@@ -81,7 +81,6 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
_law=law;
_residualVolume=new dgResidualVolume(*_law);
_residualInterface=new dgResidualInterface(*_law);
- _residualBoundary=new dgResidualBoundary(*_law);
_gc=gc;
_K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){
@@ -113,32 +112,15 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
}
for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){
dgGroupOfFaces *gf=gc->getFaceGroup(iGroup);
- const dgGroupOfElements *ge[2];
- ge[0]=&gf->getGroupOfConnections(0).getGroupOfElements();
- ge[1]=&gf->getGroupOfConnections(1).getGroupOfElements();
- for(int i=0;i<2;i++){
- if(ge[i]->getIsInnerMultirateBuffer()){
- _innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
+ for(int i=0;i<gf->getNbGroupOfConnections();i++){
+ const dgGroupOfElements *ge = &gf->getGroupOfConnections(0).getGroupOfElements();
+ if(ge->getIsInnerMultirateBuffer()){
+ _innerBufferGroupsOfElements[ge->getMultirateExponent()].second.push_back(gf);
}
- else if(ge[i]->getIsOuterMultirateBuffer()){
- _outerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
+ else if(ge->getIsOuterMultirateBuffer()){
+ _outerBufferGroupsOfElements[ge->getMultirateExponent()].second.push_back(gf);
}else{
- _bulkGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
- }
- }
- }
- for(int iGroup=0;iGroup<gc->getNbBoundaryGroups();iGroup++){
- dgGroupOfFaces *gf=gc->getBoundaryGroup(iGroup);
- const dgGroupOfElements *ge[1];
- ge[0]=&gf->getGroupOfConnections(0).getGroupOfElements();
- for(int i=0;i<1;i++){
- if(ge[i]->getIsInnerMultirateBuffer()){
- _innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
- }
- else if(ge[i]->getIsOuterMultirateBuffer()){
- _outerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
- }else{
- _bulkGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
+ _bulkGroupsOfElements[ge->getMultirateExponent()].second.push_back(gf);
}
}
}
@@ -163,7 +145,6 @@ dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
delete _currentInput;
delete _residualVolume;
delete _residualInterface;
- delete _residualBoundary;
}
void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
@@ -191,7 +172,7 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
for(std::set<dgGroupOfFaces *>::iterator it=gOFSet.begin();it!=gOFSet.end();it++){
dgGroupOfFaces *faces=*it;
if(faces->getNbGroupOfConnections()==1){
- _residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual);
+ _residualInterface->computeAndMap1Group(*faces,*_currentInput,*_residual);
}
else{
const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements();
@@ -243,7 +224,7 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
for(int i=0;i<vf.size();i++){
dgGroupOfFaces *faces=vf[i];
if(faces->getNbGroupOfConnections()==1){
- _residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual);
+ _residualInterface->computeAndMap1Group(*faces,*_currentInput,*_residual);
}
else{
const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements();
diff --git a/Solver/dgRungeKuttaMultirate.h b/Solver/dgRungeKuttaMultirate.h
index 79004a329a..91254fbe57 100644
--- a/Solver/dgRungeKuttaMultirate.h
+++ b/Solver/dgRungeKuttaMultirate.h
@@ -12,7 +12,6 @@ class dgRungeKuttaMultirate: public dgRungeKutta{
dgDofContainer *_residual;
dgResidualVolume *_residualVolume;
dgResidualInterface *_residualInterface;
- dgResidualBoundary *_residualBoundary;
dgGroupCollection *_gc;
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_bulkGroupsOfElements;// int is the multirateExponent
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_innerBufferGroupsOfElements;// int is the multirateExponent
--
GitLab