From 174749368663d8f1fd69c6017ec664aed3a1ccec Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 4 Nov 2009 08:54:51 +0000
Subject: [PATCH]
---
Solver/dgAlgorithm.cpp | 6 +++---
Solver/dgGroupOfElements.h | 28 ++++++++++++++++++++--------
2 files changed, 23 insertions(+), 11 deletions(-)
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index faac9de1fb..447be98377 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -35,8 +35,8 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
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);
+ // proxy to the elementary values
+ fullMatrix<double> fuvwe(Fuvw[iUVW], iElement*claw.nbFields(),claw.nbFields());
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
@@ -46,7 +46,7 @@ void dgAlgorithm::residualVolume ( dofManager &dof, // the DOF manager (maybe us
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;
+ fuvwe(iPt,k) += ( fConv[iXYZ](iPt,k) + fDiff[iXYZ](iPt,k)) * factor;
}
}
}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index b211fc110f..6f1922ad9a 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -20,9 +20,14 @@ class functionSpace;
class dgElement {
MElement *_element;
- double *_detJ;
+ const fullMatrix<double> &_solution, &_detJ, &_integration;
public:
- dgElement (MElement *e, double *detJ);
+ dgElement (MElement *e,
+ const fullMatrix<double> &sol,
+ const fullMatrix<double> &detJ,
+ const fullMatrix<double> &integ)
+ : _element(e), _solution(sol), _detJ(detJ), _integration(integ)
+ {}
};
class dgGroupOfElements {
@@ -50,24 +55,31 @@ class dgGroupOfElements {
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;
+ int _dimUVW, _dimXYZ, _nbFields;
// forbid the copy
// dgGroupOfElements (const dgGroupOfElements &e, int order) {}
// dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
public:
- dgGroupOfElements (const std::vector<MElement*> &e, int);
+ dgGroupOfElements (const std::vector<MElement*> &e, int pOrder);
virtual ~dgGroupOfElements ();
- inline dgElement getDGElement(int i) const;
+ 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 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;}
+ // 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;}
- inline double getDetJ (int iElement, int iGaussPoint)
+ // get a proxy on the residual for element iElement
+ inline fullMatrix<double> & getResidual (int iElement) const {return fullMatrix<double>(*_residual, iElement*_nbFields, _nbFields);}
+ 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> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
};
/*class dgFace {
--
GitLab