diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index abe4b8c5fb33045e246534bab257df4a99046d2d..551e0e9e8db61481da2e41ddffd7a23e1957569e 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -57,6 +57,8 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
gradientSolutionQP = new fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields * 3);
//group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP);
}
+ fullMatrix<double> fuvwe;
+ fullMatrix<double> source;
// ----- 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
@@ -68,7 +70,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
// ----- 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
- fullMatrix<double> fuvwe(Fuvw[iUVW], iElement*nbFields,nbFields);
+ fuvwe.setAsProxy(Fuvw[iUVW], iElement*nbFields,nbFields);
// ----- 2.3.3.1 --- compute \vec{f}^{UVW} =\vec{f}^{XYZ} J^{-1} and multiply bt detJ
for (int iXYZ=0;iXYZ<group.getDimXYZ();iXYZ++) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
@@ -89,7 +91,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
}
}
if (sourceTerm){
- fullMatrix<double> source(Source, iElement*claw.nbFields(),claw.nbFields());
+ source.setAsProxy(Source, iElement*claw.nbFields(),claw.nbFields());
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iElement, iPt);
for (int k=0;k<nbFields;k++){
@@ -162,6 +164,7 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
dataCacheDouble *riemannSolver = claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
dataCacheDouble *diffusiveFluxLeft = claw.newDiffusiveFlux(cacheMapLeft);
dataCacheDouble *diffusiveFluxRight = claw.newDiffusiveFlux(cacheMapRight);
+ fullMatrix<double> dPsiLeftDx,dPsiRightDx,dofsLeft,dofsRight,normalFluxQP;
for (int iFace=0 ; iFace < nbFaces ; ++iFace) {
// B1 ) adjust the proxies for this element
@@ -173,15 +176,14 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
// proxies for the solution
solutionQPLeft.setAsProxy(solutionQP, iFace*2*nbFields, nbFields);
solutionQPRight.setAsProxy(solutionQP, (iFace*2+1)*nbFields, nbFields);
- // this should be changed to use a proxy instead of copying the normals values
- normals.set(group.getNormals(iFace));
+ normals.setAsProxy(group.getNormals(), iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
// proxies needed to compute the gradient of the solution and the IP term
- fullMatrix<double> dPsiLeftDx (DPsiLeftDx,iFace*nbNodesLeft,nbNodesLeft);
- fullMatrix<double> dPsiRightDx (DPsiLeftDx,iFace*nbNodesRight,nbNodesRight);
- fullMatrix<double> dofsLeft (solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
- fullMatrix<double> dofsRight (solutionRight, nbFields*group.getElementRightId(iFace), nbFields);
+ dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*nbNodesLeft,nbNodesLeft);
+ dPsiRightDx.setAsProxy(DPsiLeftDx,iFace*nbNodesRight,nbNodesRight);
+ dofsLeft.setAsProxy(solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
+ dofsRight.setAsProxy(solutionRight, nbFields*group.getElementRightId(iFace), nbFields);
// proxies for the flux
- fullMatrix<double> normalFluxQP (NormalFluxQP, iFace*nbFields*2, nbFields*2);
+ normalFluxQP.setAsProxy(NormalFluxQP, iFace*nbFields*2, nbFields*2);
// B2 ) compute the gradient of the solution
if(gradientSolutionLeft.somethingDependOnMe()){
@@ -218,10 +220,14 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
fullMatrix<double> &residu,
fullMatrix<double> &sol)
{
+ fullMatrix<double> residuEl;
+ fullMatrix<double> solEl;
+ fullMatrix<double> iMassEl;
for(int i=0;i<group.getNbElements();i++) {
- fullMatrix<double> residuEl(residu,i,1);
- fullMatrix<double> solEl(sol,i,1);
- solEl.gemm(group.getInverseMassMatrix(i),residuEl);
+ residuEl.setAsProxy(residu,i,1);
+ solEl.setAsProxy(sol,i,1);
+ iMassEl.setAsProxy(group.getInverseMassMatrix(),i*group.getNbNodes(),group.getNbNodes());
+ solEl.gemm(iMassEl,residuEl);
}
}
@@ -246,23 +252,26 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
fullMatrix<double> K(sol);
// Current updated solution
fullMatrix<double> Unp(sol);
+ fullMatrix<double> residuEl, KEl;
+ fullMatrix<double> iMassEl;
+
+ int nbNodes=eGroups[0]->getNbNodes();
for(int j=0; j<orderRK;j++){
- if(j!=0){
- K.scale(b[j]);
- K.add(sol);
- }
- this->residual(claw,eGroups,fGroups,bGroups,K,residu);
- K.scale(0.);
- for(int i=0;i<eGroups[0]->getNbElements();i++) {
- fullMatrix<double> residuEl(residu,i,1);
- fullMatrix<double> KEl(K,i,1);
- (eGroups[0]->getInverseMassMatrix(i)).mult(residuEl,KEl);
+ if(j!=0){
+ K.scale(b[j]);
+ K.add(sol);
+
+ this->residual(claw,eGroups,fGroups,bGroups,K,residu);
+ K.scale(0.);
+ for(int i=0;i<eGroups[0]->getNbElements();i++) {
+ residuEl.setAsProxy(residu,i,1);
+ KEl.setAsProxy(K,i,1);
+ iMassEl.setAsProxy(eGroups[0]->getInverseMassMatrix(),i*nbNodes,nbNodes);
+ iMassEl.mult(residuEl,KEl);
}
-
- Unp.add(K,a[j]);
+ Unp.add(K,a[j]);
}
-
sol=Unp;
}
@@ -291,15 +300,16 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
// required data
dataCacheDouble *boundaryTerm = boundaryCondition->newBoundaryTerm(cacheMapLeft);
+ fullMatrix<double> normalFluxQP;
for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
+ normalFluxQP.setAsProxy(NormalFluxQP, iFace*nbFields, nbFields);
// ----- 2.3.1 --- provide the data to the cacheMap
- solutionQPLeft.set(fullMatrix<double>(solutionQP, iFace*nbFields, nbFields));
- normals.set(group.getNormals(iFace));
+ solutionQPLeft.setAsProxy(solutionQP, iFace*nbFields, nbFields);
+ normals.setAsProxy(group.getNormals(),iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
cacheElementLeft.set(group.getElementLeft(iFace));
// ----- 2.3.2 --- compute fluxes
- 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++)
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 45d51672836e7f6b01228d8bd7a143eee79334c9..7e23dfe336bcd01a5c2e0328725491b8ebb31800 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -79,7 +79,7 @@ public:
inline const fullMatrix<double> & getSourceRedistributionMatrix () const {return *_redistributionSource;}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_mapping)(iElement, 10*iGaussPoint + 9);}
inline double getInvJ (int iElement, int iGaussPoint, int i, int j) const {return (*_mapping)(iElement, 10*iGaussPoint + i + 3*j);}
- inline fullMatrix<double> getInverseMassMatrix (int iElement) const {return fullMatrix<double>(*_imass,iElement*getNbNodes(),getNbNodes());}
+ inline fullMatrix<double> &getInverseMassMatrix () const {return *_imass;}
inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
};
@@ -129,7 +129,7 @@ public:
inline MElement* getFace (int iElement) const {return _faces[iElement];}
const std::vector<int> * getClosureLeft(int iFace) const{ return _closuresLeft[iFace];}
const std::vector<int> * getClosureRight(int iFace) const{ return _closuresRight[iFace];}
- inline fullMatrix<double> getNormals (int iFace) const {return fullMatrix<double>(*_normals,iFace*getNbIntegrationPoints(),getNbIntegrationPoints());}
+ inline fullMatrix<double> &getNormals () const {return *_normals;}
dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces);