No commit message

No commit message

Solver/dgAlgorithm.cpp

+#include "dgAlgorithm.h"
+
+/*
+  compute 
+    \int \vec{f} \cdot \grad \phi dv   
+*/
+
+void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe useless here)
+				   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); 
+  
+  // for all elements of the group,
+  // compute all fluxes, both diffusive and convective, in reference 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() )};
+  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
+  for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {        
+    dgElement DGE = group.getDGElement(iElement);
+    // compute convective and diffusive fluxes in XYZ coordinates
+    if (claw.convectiveFlux()) (*claw.convectiveFlux())(DGE,fConv);
+    if (claw.diffusiveFlux()) (*claw.diffusiveFlux())(DGE,fDiff);
+    
+    for (int iUVW=0;iUVW<group.getDimUVW();iUVW++) {
+      // proxies should be done...
+      fullMatrix<double> fuvwe  = Fuvw[iUVW].linesProxy(iElement);
+      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
+	  const double invJ = group.getInvJ (iElement, iPt, iXYZ, iUVW);
+	  // get the detJ at this point
+	  const double detJ = group.getDetJ (iElement, iPt);
+	  const double factor = invJ * detJ;
+	  // compute fluxes in the reference coordinate system at this integration point
+	  for (int k=0;k<group.nbFields();k++) {
+	    fuvwe(i,k) += ( fConv[iXYZ](i,k) + fDiff[iXYZ](i,k)) * factor;
+	  }
+	}
+      } 
+    }
+  }
+  // finally, premultiply with the redistribution matrices
+  // this is the "heavy load" part
+  for (int iUVW=0;iUVW<group.getDimUVW();iUVW++)
+    group.getResidual().dgemm(group.getCollocationMatrix(iUVW), Fuvw[iUVW]);
+}

Solver/dgConservationLaw.h

@@ -13,13 +13,13 @@ class dgFace;
 class dgTerm{
  public:
   virtual ~dgTerm () {}
-  virtual void operator () (const dgElement &, fullMatrix<double> &fcx) const = 0;
+  virtual void operator () (const dgElement &, fullMatrix<double> fcx[]) const = 0;
 };
 
 class dgFaceTerm{
  public:
   virtual ~dgFaceTerm () {}
-  virtual void operator () (const dgFace &, fullMatrix<double> &fcx) const = 0;
+  virtual void operator () (const dgFace &, fullMatrix<double> fcx[]) const = 0;
 };
 
 class dgConservationLaw {

Solver/dgGroupOfElements.h

@@ -18,6 +18,13 @@ class MFace;
 class MEdge;
 class functionSpace;
 
+class dgElement {
+  MElement *_element;
+  double *_detJ;
+public:
+  dgElement (MElement *e, double *detJ);
+};
+
 class dgGroupOfElements {
   // N elements in the group
   std::vector<MElement*> _elements;
@@ -37,12 +44,30 @@ class dgGroupOfElements {
   fullMatrix<double> *_redistributionFluxes[3];
   // redistribution for the source term
   fullMatrix<double> *_redistributionSource;
+  // the solution for this group (not owned)
+  fullMatrix<double> *_solution;
+  // the residual for this group (not owned)
+  fullMatrix<double> *_residual;
+  // dimension of the parametric space and of the real space
+  // may be different if the domain is a surface in 3D (manifold)
+  int dimUVW, dimXYZ;
   // forbid the copy 
   //  dgGroupOfElements (const dgGroupOfElements &e, int order) {}
   //  dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
 public:
   dgGroupOfElements (const std::vector<MElement*> &e, int);
   virtual ~dgGroupOfElements ();
+  inline dgElement getDGElement(int i) const;
+  inline int getNbElements() const {return _elements.size();}
+  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 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> & getResidual () const {return *_solution;}
+  inline double getDetJ (int iElement, int iGaussPoint) 
 };
 
 /*class dgFace {