From 240d806418f5a5c32b55f4d0afa24cac43c0bee1 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 3 Nov 2009 21:29:09 +0000
Subject: [PATCH]
---
Solver/dgAlgorithm.cpp | 59 ++++++++++++++++++++++++++++++++++++++
Solver/dgConservationLaw.h | 4 +--
Solver/dgGroupOfElements.h | 25 ++++++++++++++++
3 files changed, 86 insertions(+), 2 deletions(-)
create mode 100644 Solver/dgAlgorithm.cpp
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
new file mode 100644
index 0000000000..faac9de1fb
--- /dev/null
+++ b/Solver/dgAlgorithm.cpp
@@ -0,0 +1,59 @@
+#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]);
+}
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 805806f3ab..8c5bc9ff33 100644
--- a/Solver/dgConservationLaw.h
+++ b/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 {
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 2b517b6773..b211fc110f 100644
--- a/Solver/dgGroupOfElements.h
+++ b/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 {
--
GitLab