Skip to content
Snippets Groups Projects
Commit 2b70dffc authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dg : compute dgPsiLeftDx and Solution gradient on faces qp

parent 0b82c444
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Numeric/fullMatrix.h
+ 8
− 0
View file @ 2b70dffc
......@@ -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)
......
Solver/dgAlgorithm.cpp
+ 75
− 20
View file @ 2b70dffc
......@@ -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);
}
}
Solver/dgAlgorithm.h
+ 13
− 10
View file @ 2b70dffc
......@@ -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);
......
Solver/dgConservationLawAdvection.cpp
+ 1
− 1
View file @ 2b70dffc
Solver/dgGroupOfElements.cpp
+ 69
− 36
View file @ 2b70dffc
......@@ -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
......@@ -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)
......
Solver/dgGroupOfElements.h
+ 10
− 23
View file @ 2b70dffc
......@@ -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();}
......
Solver/function.h
+ 15
− 0
View file @ 2b70dffc
......@@ -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) {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment