diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 447be9837702b26c603a91fbb924cae8d4d865d7..6fe53fa02be80617b4987fddb6ae3b20a80e8ea6 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -9,34 +9,48 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
const dgConservationLaw &claw, // the conservation law
const groupOfElements & group ) // the residual
{
- // get the solution at quadrature points
- fullMatrix<double> solutionQp (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints());
- // solutionQp ( NbE x NbFields , NbQuad)
- // collocation ( NbNod , NbQuad)
- // solution ( NbE x NbFields, nbNodes)
- group.getSolution().mult ( group.getCollocationMatrix() , solutionQp);
+ // ----- 1 ---- get the solution at quadrature points
+ // ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints)
+ fullMatrix<double> solutionQP (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints());
+ // ----- 1.2 --- multiply the solution by the collocation matrix
+ group.getSolution().mult ( group.getCollocationMatrix() , 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.getNbElements() * group.getNbFields() * 3, group.getNbIntegrationPoints());
+ group.getGradientOfSolution().mult ( group.getCollocationMatrix() , solutionQP);
+ }
- // for all elements of the group,
- // compute all fluxes, both diffusive and convective, in reference coordinates
+
+ // ----- 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(), group.nbFields() ),
fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() )};
fullMatrix<double> fDiff[3] = {fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() ),
fullMatrix<double>( group.getNbIntegrationPoints(), group.nbFields() )};
+ // ----- 2.2 --- allocate parametric fluxes (compued in UVW coordinates) for all elements at all integration points
fullMatrix<double> Fuvw[3] = {fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints()),
fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints()),
fullMatrix<double> (group.getNbElements() * group.getNbFields(), group.getNbIntegrationPoints())};
- // for all elements of the group
+ // ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
- dgElement DGE = group.getDGElement(iElement);
- // compute convective and diffusive fluxes in XYZ coordinates
+ // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
+ fullMatrix<double> solutionQPe (solutionQP, iElement*claw.nbFields(),claw.nbFields() );
+ fullMatrix<double> *gradSolutionQPe;
+ if (claw.diffusiveFlux()) gradSolutionQPe = new fullMatrix<double>(*gradSolutionQP, 3*iElement*claw.nbFields(),3*claw.nbFields() );
+ else gradSolutionQPe = new fullMatrix<double>
+ dgElement DGE( group.getElement(iElement), solutionQPe, *gradSolutionQPe, group.getIntegrationPointsMatrix());
+ // ----- 2.3.2 --- compute fluxes in XYZ coordinates
if (claw.convectiveFlux()) (*claw.convectiveFlux())(DGE,fConv);
if (claw.diffusiveFlux()) (*claw.diffusiveFlux())(DGE,fDiff);
+ // ----- 2.3.3 --- compute fluxes in UVW coordinates
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++) {
- // proxy to the elementary values
+ // ----- 2.3.3.1 --- get a proxy on the big local flux matrix
fullMatrix<double> fuvwe(Fuvw[iUVW], iElement*claw.nbFields(),claw.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++) {
// get the right inv jacobian matrix dU/dX element
@@ -52,8 +66,7 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
}
}
}
- // finally, premultiply with the redistribution matrices
- // this is the "heavy load" part
+ // ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++)
- group.getResidual().dgemm(group.getCollocationMatrix(iUVW), Fuvw[iUVW]);
+ group.getResidual().dgemm(group.getRedistributionMatrix(iUVW), Fuvw[iUVW]);
}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 6f1922ad9ac7f8d613cb3ce5a2fa9de935f7eacd..3f793ea486915ac5b544a2ddd1bf95946bee76cf 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -20,13 +20,14 @@ class functionSpace;
class dgElement {
MElement *_element;
- const fullMatrix<double> &_solution, &_detJ, &_integration;
+ // solution at points
+ const fullMatrix<double> &_solution, &_integration, &_gradients;
public:
- dgElement (MElement *e,
- const fullMatrix<double> &sol,
- const fullMatrix<double> &detJ,
- const fullMatrix<double> &integ)
- : _element(e), _solution(sol), _detJ(detJ), _integration(integ)
+ dgElement (MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &integ)
+ : _element(e), _solution(sol), _integration(integ), _gradients(sol)
+ {}
+ dgElement (MElement *e, const fullMatrix<double> &sol, const fullMatrix<double> &grads, const fullMatrix<double> &integ)
+ : _element(e), _solution(sol), _integration(integ), _gradients(grad)
{}
};
@@ -49,6 +50,8 @@ class dgGroupOfElements {
fullMatrix<double> *_redistributionFluxes[3];
// redistribution for the source term
fullMatrix<double> *_redistributionSource;
+ // the "finite element" gradient of the solution for this group (not owned)
+ fullMatrix<double> *_gradSolution;
// the solution for this group (not owned)
fullMatrix<double> *_solution;
// the residual for this group (not owned)
@@ -62,16 +65,18 @@ class dgGroupOfElements {
public:
dgGroupOfElements (const std::vector<MElement*> &e, int pOrder);
virtual ~dgGroupOfElements ();
- inline dgElement getDGElement(int i) const {return dgElement(_elements[i], getSolution(i), getMapping(i), *_integration);}
inline int getNbElements() const {return _elements.size();}
inline int getNbFields() const {return _nbFields;}
inline int getNbNodes() const {return _collocation.size1();}
inline int getNbIntegrationPoints() const {return _collocation.size2();}
inline int getDimUVW () const {return _dimUVW;}
inline int getDimXYZ () const {return _dimXYZ;}
+ inline const MElement* getElement (int iElement) const {return _elements[iElement];}
+ inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
inline const fullMatrix<double> & getRedistributionMatrix (int i) const {return *_redistribution[i];}
inline const fullMatrix<double> & getSolution () const {return *_solution;}
+ inline const fullMatrix<double> & getGradientOfSolution () const {return *_gradSolution;}
// get a proxy on the solution for element iElement
inline fullMatrix<double> & getSolution (int iElement) const {return fullMatrix<double>(*_solution, iElement*_nbFields, _nbFields);}
inline const fullMatrix<double> & getResidual () const {return *_solution;}