diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index bc829bd231b0accfa99f840253acb5a7d0e2e94d..b4561cb25cd847ba7c556af82d7c3de6627e57c2 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -103,6 +103,14 @@ class fullMatrix
}
fullMatrix() : _own_data(false),_r(0), _c(0), _data(0) {}
~fullMatrix() { if(_data && _own_data) delete [] _data; }
+ void setAsProxy(fullMatrix<scalar> &original, int c_start, int c) {
+ if(_data && _own_data)
+ delete [] _data;
+ _c = c;
+ _r = original._r;
+ _own_data = false;
+ _data = original._data + c_start * _r;
+ }
inline int size1() const { return _r; }
inline int size2() const { return _c; }
fullMatrix<scalar> & operator = (const fullMatrix<scalar> &other)
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 5aaaf3e1de364f275a56a40c964474a995dc20aa..abe4b8c5fb33045e246534bab257df4a99046d2d 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -60,9 +60,9 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
// ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
- solutionQPe.set(fullMatrix<double>(solutionQP, iElement*nbFields, nbFields ));
+ solutionQPe.setAsProxy(solutionQP, iElement*nbFields, nbFields );
if (gradSolutionQPe.somethingDependOnMe())
- gradSolutionQPe.set(fullMatrix<double>(*gradientSolutionQP, 3*iElement*nbFields,3*nbFields ));
+ gradSolutionQPe.setAsProxy(*gradientSolutionQP, 3*iElement*nbFields,3*nbFields );
cacheElement.set(group.getElement(iElement));
if(convectiveFlux || diffusiveFlux) {
// ----- 2.3.3 --- compute fluxes in UVW coordinates
@@ -120,43 +120,97 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
const dgConservationLaw &claw, // the conservation law
const dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &solutionRight,
fullMatrix<double> &residual // residual !! at faces nodes
)
{
+ // A) global operations before entering the loop over the faces
+ // A1 ) copy locally some references from the group for the sake of lisibility
int nbFields = claw.nbFields();
+ int nbNodesLeft = group.getGroupLeft().getNbNodes();
+ int nbNodesRight = group.getGroupRight().getNbNodes();
+ int nbFaces = group.getNbElements();
+ //get matrices needed to compute the gradient on both sides
+ fullMatrix<double> &DPsiLeftDx = group.getDPsiLeftDxMatrix();
+ fullMatrix<double> &DPsiRightDx = group.getDPsiRightDxMatrix();
+
// ----- 1 ---- get the solution at quadrature points
- fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields*2);
+ fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),nbFaces * nbFields*2);
group.getCollocationMatrix().mult(solution, solutionQP);
- // ----- 2 ---- compute normal fluxes at integration points
- fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields*2);
+ // needed tocompute normal fluxes at integration points
+ fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), nbFaces*nbFields*2);
+
+ // create one dataCache for each side
dataCacheMap cacheMapLeft, cacheMapRight;
- // provided dataCache
- cacheMapLeft.provideData("UVW").set(group.getIntegrationPointsMatrix());
- dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
- dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution");
+
+ // data this algorithm provide to the cache map (we can maybe move each data to a separate function but
+ // I think It's easier like this)
dataCacheDouble &normals = cacheMapLeft.provideData("Normals");
dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
dataCacheElement &cacheElementRight = cacheMapRight.getElement();
- // required data
+ cacheMapLeft.provideData("UVW").set(group.getIntegrationPointsMatrix());
+ dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
+ dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution");
+ dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("GradientSolution");
+ dataCacheDouble &gradientSolutionRight = cacheMapRight.provideData("GradientSolution");
+ //resize gradientSolutionLeft and Right
+ gradientSolutionLeft.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
+ gradientSolutionRight.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
+
+ // A2 ) ask the law to provide the fluxes (possibly NULL)
dataCacheDouble *riemannSolver = claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
- for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
- // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
- solutionQPLeft.set(fullMatrix<double>(solutionQP, iFace*2*nbFields, nbFields));
- solutionQPRight.set(fullMatrix<double>(solutionQP, (iFace*2+1)*nbFields, nbFields));
- normals.set(group.getNormals(iFace));
+ dataCacheDouble *diffusiveFluxLeft = claw.newDiffusiveFlux(cacheMapLeft);
+ dataCacheDouble *diffusiveFluxRight = claw.newDiffusiveFlux(cacheMapRight);
+
+ for (int iFace=0 ; iFace < nbFaces ; ++iFace) {
+ // B1 ) adjust the proxies for this element
+ // NB : the B1 section will almost completely disapear
+ // if we write a function (from the function class) for moving proxies across big matrices
+ // give the current elements to the dataCacheMap
cacheElementLeft.set(group.getElementLeft(iFace));
cacheElementRight.set(group.getElementRight(iFace));
+ // proxies for the solution
+ solutionQPLeft.setAsProxy(solutionQP, iFace*2*nbFields, nbFields);
+ solutionQPRight.setAsProxy(solutionQP, (iFace*2+1)*nbFields, nbFields);
+ // this should be changed to use a proxy instead of copying the normals values
+ normals.set(group.getNormals(iFace));
+ // proxies needed to compute the gradient of the solution and the IP term
+ fullMatrix<double> dPsiLeftDx (DPsiLeftDx,iFace*nbNodesLeft,nbNodesLeft);
+ fullMatrix<double> dPsiRightDx (DPsiLeftDx,iFace*nbNodesRight,nbNodesRight);
+ fullMatrix<double> dofsLeft (solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
+ fullMatrix<double> dofsRight (solutionRight, nbFields*group.getElementRightId(iFace), nbFields);
+ // proxies for the flux
fullMatrix<double> normalFluxQP (NormalFluxQP, iFace*nbFields*2, nbFields*2);
- // ----- 2.3.2 --- compute fluxes
+
+ // B2 ) compute the gradient of the solution
+ if(gradientSolutionLeft.somethingDependOnMe()){
+ dPsiLeftDx.mult(dofsLeft, gradientSolutionLeft.set());
+ dPsiRightDx.mult(dofsRight, gradientSolutionRight.set());
+ }
+
+ // B3 ) compute fluxes
+
+ //JF : here you can test (diffusiveFluxLeft!=NULL) to know if the law is diffusive
+ //you can access the values by (*diffusiveFluxLeft)(iQP*3+iXYZ, iField)
+ //use gradientSolutionLef(IQP*3+IXYZ, iField) and dPsiLeftDx(iQP*3+iXYZ, iPsi)
+
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
for (int k=0;k<nbFields*2;k++)
normalFluxQP(iPt,k) = (*riemannSolver)(iPt,k)*detJ;
}
}
- // ----- 3 ---- do the redistribution at face nodes using BLAS3
+
+ // C ) redistribute the flux to the residual (at Faces nodes)
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
+
+ // D ) delete the dataCacheDouble provided by the law
delete riemannSolver;
+ if(diffusiveFluxLeft) {
+ delete diffusiveFluxLeft;
+ delete diffusiveFluxRight;
+ }
}
void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
@@ -217,6 +271,7 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
const dgConservationLaw &claw, // the conservation law
const dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
)
{
@@ -304,7 +359,7 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
std::vector<dgGroupOfElements*> &eGroups, //group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- const fullMatrix<double> &solution, // solution
+ fullMatrix<double> &solution, // solution
fullMatrix<double> &residu) // residual
{
dgAlgorithm algo;
@@ -319,7 +374,7 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2);
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2);
faces.mapToInterface(1, solution, solution, solInterface);
- algo.residualInterface(claw,faces,solInterface,residuInterface);
+ algo.residualInterface(claw,faces,solInterface,solution,solution,residuInterface);
faces.mapFromInterface(1, residuInterface, residu, residu);
}
//boundaries
@@ -328,7 +383,7 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements());
fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements());
faces.mapToInterface(1, solution, solution, solInterface);
- algo.residualBoundary(claw,faces,solInterface,residuInterface);
+ algo.residualBoundary(claw,faces,solInterface,solution,residuInterface);
faces.mapFromInterface(1, residuInterface, residu, residu);
}
}
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 4d1c81bb080f3a260745d19bc1380f400def381e..f51f050e423634c154b1258b472f86b835b0a0ef 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -11,34 +11,33 @@ class dgTerm;
class dgAlgorithm {
public :
- void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
+ static void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group,
const fullMatrix<double> &solution,
fullMatrix<double> &residual);
- void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
+ static void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfFaces & group,
const fullMatrix<double> &solution, // ! at face nodes
+ fullMatrix<double> &solutionRight,
+ fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual); // ! at face nodes
- void residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
+ static void residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfFaces &group,
const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
);
// works only if there is only 1 group of element
- void residual( const dgConservationLaw &claw,
+ static void residual( const dgConservationLaw &claw,
std::vector<dgGroupOfElements*> &eGroups, //group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solution, // solution !! at faces nodes
fullMatrix<double> &residual // residual !! at faces nodes
);
- void multAddInverseMassMatrix ( /*dofManager &dof,*/
- const dgGroupOfElements & group,
- fullMatrix<double> &residu,
- fullMatrix<double> &sol);
void rungeKutta (
const dgConservationLaw &claw,
std::vector<dgGroupOfElements*> &eGroups, //group of elements
@@ -48,7 +47,11 @@ class dgAlgorithm {
fullMatrix<double> &residu,
fullMatrix<double> &sol,
int orderRK=4);
- void buildGroups(GModel *model, int dim, int order,
+ static void multAddInverseMassMatrix ( /*dofManager &dof,*/
+ const dgGroupOfElements & group,
+ fullMatrix<double> &residu,
+ fullMatrix<double> &sol);
+ static void buildGroups(GModel *model, int dim, int order,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroups,
std::vector<dgGroupOfFaces*> &bGroups);
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index a886ef5bf8aeb7c0eb2c3e1088eb1485b8de8d7d..ce654004be7002c474f720a3a34b58ca59ccdaf3 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -16,7 +16,7 @@ class dgConservationLawAdvection : public dgConservationLaw {
v(cacheMap.get(vFunctionName,this))
{};
void _eval () {
- if(_value.size1() !=sol().size1())
+ if(_value.size1() != sol().size1())
_value=fullMatrix<double>(sol().size1(),3);
for(int i=0; i< sol().size1(); i++) {
_value(i,0) = sol(i,0)*v(i,0);
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 5b28f8acbea95df7043dcc448f284768369edf70..7694e1c065ef8709f0f01fb414203c7183074be5 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -119,40 +119,6 @@ dgGroupOfElements::~dgGroupOfElements(){
}
-void dgGroupOfFaces::computeFaceNormals () {
- double g[256][3];
- _normals = new fullMatrix<double> (3,_fsFace->points.size1()*_faces.size());
- int index = 0;
- for (size_t i=0; i<_faces.size();i++){
- const std::vector<int> &closure=*_closuresLeft[i];
- fullMatrix<double> *intLeft=dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,&closure);
- double jac[3][3],ijac[3][3];
- for (int j=0; j<intLeft->size1(); j++){
- _fsLeft->df((*intLeft)(j,0),(*intLeft)(j,1),(*intLeft)(j,2),g);
- getElementLeft(i)->getJacobian ((*intLeft)(j,0), (*intLeft)(j,1), (*intLeft)(j,2), jac);
- inv3x3(jac,ijac);
- double &nx=(*_normals)(0,index);
- double &ny=(*_normals)(1,index);
- double &nz=(*_normals)(2,index);
- double nu=0,nv=0,nw=0;
- for (size_t k=0; k<closure.size(); k++){
- nu += g[closure[k]][0];
- nv += g[closure[k]][1];
- nw += g[closure[k]][2];
- }
- nx = nu*ijac[0][0]+nv*ijac[0][1]+nw*ijac[0][2];
- ny = nu*ijac[1][0]+nv*ijac[1][1]+nw*ijac[1][2];
- nz = nu*ijac[2][0]+nv*ijac[2][1]+nw*ijac[2][2];
- double norm = sqrt(nx*nx+ny*ny+nz*nz);
- nx/=norm;
- ny/=norm;
- nz/=norm;
- index++;
- }
- delete intLeft;
- }
-}
-
void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
// compute all closures
// compute closures for the interpolation
@@ -221,7 +187,7 @@ void dgGroupOfFaces::addEdge(const MEdge &topoEdge, int iElLeft, int iElRight){
void dgGroupOfFaces::init(int pOrder) {
_fsFace = _faces[0]->getFunctionSpace (pOrder);
- _integration=dgGetIntegrationRule (_faces[0],pOrder);
+ _integration = dgGetIntegrationRule (_faces[0],pOrder);
_redistribution = new fullMatrix<double> (_fsFace->coefficients.size1(),_integration->size1());
_collocation = new fullMatrix<double> (_integration->size1(), _fsFace->coefficients.size1());
_detJac = new fullMatrix<double> (_integration->size1(), _faces.size());
@@ -241,7 +207,71 @@ void dgGroupOfFaces::init(int pOrder) {
(*_detJac)(j,i)=f->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
}
}
- computeFaceNormals();
+
+ // compute data on quadrature points : normals and dPsidX
+ double g[256][3];
+ _normals = new fullMatrix<double> (3,_fsFace->points.size1()*_faces.size());
+ _dPsiLeftDxOnQP = new fullMatrix<double> ( _integration->size1()*3,_fsLeft->points.size1()*_faces.size());
+ if(_fsRight){
+ _dPsiRightDxOnQP = new fullMatrix<double> ( _integration->size1()*3,_fsRight->points.size1()*_faces.size());
+ }
+ int index = 0;
+ for (size_t i=0; i<_faces.size();i++){
+ const std::vector<int> &closureLeft=*_closuresLeft[i];
+ fullMatrix<double> *intLeft=dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,&closureLeft);
+ double jac[3][3],ijac[3][3];
+ for (int j=0; j<intLeft->size1(); j++){
+ _fsLeft->df((*intLeft)(j,0),(*intLeft)(j,1),(*intLeft)(j,2),g);
+ getElementLeft(i)->getJacobian ((*intLeft)(j,0), (*intLeft)(j,1), (*intLeft)(j,2), jac);
+ inv3x3(jac,ijac);
+ //compute dPsiLeftDxOnQP
+ //(iPsi*3+iXYZ,iQP+iFace*NQP);
+ int nPsi = _fsLeft->coefficients.size1();
+ for (int iPsi=0; iPsi< nPsi; iPsi++) {
+ (*_dPsiLeftDxOnQP)(j*3 ,i*nPsi+iPsi) = g[iPsi][0]*ijac[0][0]+g[iPsi][1]*ijac[0][1]+g[iPsi][2]*ijac[0][2];
+ (*_dPsiLeftDxOnQP)(j*3+1,i*nPsi+iPsi) = g[iPsi][0]*ijac[1][0]+g[iPsi][1]*ijac[1][1]+g[iPsi][2]*ijac[1][2];
+ (*_dPsiLeftDxOnQP)(j*3+2,i*nPsi+iPsi) = g[iPsi][0]*ijac[2][0]+g[iPsi][1]*ijac[2][1]+g[iPsi][2]*ijac[2][2];
+ }
+ //compute face normals
+ double &nx=(*_normals)(0,index);
+ double &ny=(*_normals)(1,index);
+ double &nz=(*_normals)(2,index);
+ double nu=0,nv=0,nw=0;
+ for (size_t k=0; k<closureLeft.size(); k++){
+ nu += g[closureLeft[k]][0];
+ nv += g[closureLeft[k]][1];
+ nw += g[closureLeft[k]][2];
+ }
+ nx = nu*ijac[0][0]+nv*ijac[0][1]+nw*ijac[0][2];
+ ny = nu*ijac[1][0]+nv*ijac[1][1]+nw*ijac[1][2];
+ nz = nu*ijac[2][0]+nv*ijac[2][1]+nw*ijac[2][2];
+ double norm = sqrt(nx*nx+ny*ny+nz*nz);
+ nx/=norm;
+ ny/=norm;
+ nz/=norm;
+ index++;
+ }
+ delete intLeft;
+ // there is nothing on the right for boundary groups
+ if(_fsRight){
+ fullMatrix<double> *intRight=dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsRight,_closuresRight[i]);
+ for (int j=0; j<intRight->size1(); j++){
+ _fsRight->df((*intRight)(j,0),(*intRight)(j,1),(*intRight)(j,2),g);
+ getElementRight(i)->getJacobian ((*intRight)(j,0), (*intRight)(j,1), (*intRight)(j,2), jac);
+ inv3x3(jac,ijac);
+ //compute dPsiRightDxOnQP
+ // (iQP*3+iXYZ , iFace*NPsi+iPsi)
+ int nPsi = _fsRight->coefficients.size1();
+ for (int iPsi=0; iPsi< nPsi; iPsi++) {
+ (*_dPsiRightDxOnQP)(j*3 ,i*nPsi+iPsi) = g[iPsi][0]*ijac[0][0]+g[iPsi][1]*ijac[0][1]+g[iPsi][2]*ijac[0][2];
+ (*_dPsiRightDxOnQP)(j*3+1,i*nPsi+iPsi) = g[iPsi][0]*ijac[1][0]+g[iPsi][1]*ijac[1][1]+g[iPsi][2]*ijac[1][2];
+ (*_dPsiRightDxOnQP)(j*3+2,i*nPsi+iPsi) = g[iPsi][0]*ijac[2][0]+g[iPsi][1]*ijac[2][1]+g[iPsi][2]*ijac[2][2];
+ }
+ }
+ delete intRight;
+ }
+ }
+
}
dgGroupOfFaces::~dgGroupOfFaces()
@@ -249,6 +279,9 @@ dgGroupOfFaces::~dgGroupOfFaces()
delete _redistribution;
delete _collocation;
delete _detJac;
+ delete _normals;
+ delete _dPsiLeftDxOnQP;
+ delete _dPsiRightDxOnQP;
}
dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges):
_groupLeft(elGroup),_groupRight(elGroup)
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index da35c02100acbff7c886136eeb5ada8f5d07b68e..45d51672836e7f6b01228d8bd7a143eee79334c9 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -83,35 +83,12 @@ public:
inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
};
-class dgFace {
- int nbGaussPoints;
- MElement *_left, *_right;
- MElement *_face;
- const fullMatrix<double> *_solutionRight, *_solutionLeft, *_integration,*_normals;
-public:
- dgFace (MElement *face,MElement *left, MElement *right,
- const fullMatrix<double> &solLeft,
- const fullMatrix<double> &solRight,
- const fullMatrix<double> &integration,
- const fullMatrix<double> &normals
- ) : _left(left), _right(right), _face(face),_solutionRight(&solRight),_solutionLeft(&solLeft),_integration(&integration),_normals(&normals)
- {}
- inline const fullMatrix<double> &solutionRight() const { return *_solutionRight; }
- inline const fullMatrix<double> &solutionLeft() const { return *_solutionLeft; }
- inline const fullMatrix<double> &integration() const { return *_integration; }
- inline const fullMatrix<double> &normals() const { return *_normals; }
- inline MElement *left() const { return _left;}
- inline MElement *right() const { return _right;}
- inline MElement *face() const { return _face;}
-};
-
class dgGroupOfFaces {
// only used if this is a group of boundary faces
std::string _boundaryTag;
const dgGroupOfElements &_groupLeft,&_groupRight;
void addFace(const MFace &topoFace, int iElLeft, int iElRight);
void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
- void computeFaceNormals();
// Two polynomialBases for left and right elements
// the group has always the same types for left and right
const polynomialBasis *_fsLeft,*_fsRight, *_fsFace;
@@ -127,6 +104,10 @@ class dgGroupOfFaces {
// GEOMETRICAL CLOSURE !!!
std::vector<const std::vector<int> * > _closuresLeft;
std::vector<const std::vector<int> * > _closuresRight;
+ // XYZ gradient of the shape functions of both elements on the integrations points of the face
+ // (iQP*3+iXYZ , iFace*NPsi+iPsi)
+ fullMatrix<double> *_dPsiLeftDxOnQP;
+ fullMatrix<double> *_dPsiRightDxOnQP;
// normals at integration points (N*Ni) x 3
fullMatrix<double> *_normals;
// detJac at integration points (N*Ni) x 1
@@ -139,6 +120,10 @@ class dgGroupOfFaces {
//common part of the 3 constructors
void init(int pOrder);
public:
+ inline const dgGroupOfElements &getGroupLeft()const {return _groupLeft; }
+ inline const dgGroupOfElements &getGroupRight()const {return _groupRight; }
+ inline int getElementLeftId (int i) const {return _left[i];};
+ inline int getElementRightId (int i) const {return _right[i];};
inline MElement* getElementLeft (int i) const {return _groupLeft.getElement(_left[i]);}
inline MElement* getElementRight (int i) const {return _groupRight.getElement(_right[i]);}
inline MElement* getFace (int iElement) const {return _faces[iElement];}
@@ -151,6 +136,8 @@ public:
virtual ~dgGroupOfFaces ();
inline bool isBoundary() const {return !_boundaryTag.empty();}
inline const std::string getBoundaryTag() const {return _boundaryTag;}
+ inline fullMatrix<double> & getDPsiLeftDxMatrix() const { return *_dPsiLeftDxOnQP;}
+ inline fullMatrix<double> & getDPsiRightDxMatrix() const { return *_dPsiRightDxOnQP;}
//this part is common with dgGroupOfElements, we should try polymorphism
inline int getNbElements() const {return _faces.size();}
inline int getNbNodes() const {return _collocation->size1();}
diff --git a/Solver/function.h b/Solver/function.h
index 4f64dfe1fcc1f5ca5f124c505aa8ba03d2685d94..28c725378a24a42e8ed9a2fbe67b1bddb761b256 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -48,6 +48,21 @@ class dataCacheDouble : public dataCache {
_value=mat;
_valid=true;
}
+ // take care if you use this you must ensure that the value pointed to are not modified
+ // without further call to setAsProxy because the dependencies won't be invalidate
+ inline void setAsProxy(fullMatrix<double> &mat, int cShift, int c) {
+ _invalidateDependencies();
+ _value.setAsProxy(mat,cShift,c);
+ _valid=true;
+ }
+ // this function is needed to be able to pass the _value to functions like gemm or mult
+ // but you cannot keep the reference to the _value, you should always use the set function
+ // to modify the _value
+ inline fullMatrix<double> &set() {
+ _invalidateDependencies();
+ _valid=true;
+ return _value;
+ }
inline const double &operator () (int i, int j)
{
if(!_valid) {