No commit message

No commit message

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]);
 }

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;}