diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 0654da1299b4bb3380d9202ce4428c825dfe55c8..be697b95cd7be8f09b2e70f4a4004ca5132fcfc8 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -84,7 +84,7 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
dataCacheDouble &sol = cacheMap.provideSolution(nbFields);
dataCacheDouble &UVW = cacheMap.provideParametricCoordinates();
- UVW.set(group.getIntegrationPointsMatrix());
+ UVW.set()=group.getIntegrationPointsMatrix();
// provided dataCache
dataCacheDouble *maxConvectiveSpeed = claw.newMaxConvectiveSpeed(cacheMap);
dataCacheDouble *maximumDiffusivity = claw.newMaximumDiffusivity(cacheMap);
@@ -98,7 +98,7 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
double l_red_sq = l_red*l_red;
DT.resize(group.getNbElements());
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
- sol.setAsProxy(solution, iElement*nbFields, nbFields);
+ sol.set().setAsProxy(solution, iElement*nbFields, nbFields);
cacheMap.setElement(group.getElement(iElement));
const double L = group.getInnerRadius(iElement);
double spectralRadius = 0.0;
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index 2b66374db7030e0adf49e3a39325566e09840d80..6d0837ed314f74a3a32dfb3f36832af9278e7ca7 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/Solver/dgConservationLaw.cpp
@@ -22,7 +22,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
}
void _eval() {
- solutionRight.set(outsideValue());
+ solutionRight.set()=outsideValue();
if(riemannSolver){
for(int i=0;i<_value.size1(); i++)
for(int j=0;j<_value.size2(); j++)
@@ -62,7 +62,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
}
void _eval() {
if(maxDif){
- solutionRight.set(outsideValue());
+ solutionRight.set()=outsideValue();
for(int i=0;i<_value.size1(); i++)
for(int j=0;j<_value.size2(); j++)
_value(i,j) = (*maxDif)(i,j);
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 2e599a6ec2d879ec3b5bdf0fbf69f2303dcf8078..61b0d671ba23558bc879fa0dcb061fdb24261900 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -229,7 +229,7 @@ void dgDofContainer::L2Projection(const function *f){
fullMatrix<double> Source = fullMatrix<double> (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
dataCacheMap cacheMap;
cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
- cacheMap.provideParametricCoordinates().set(group.getIntegrationPointsMatrix());
+ cacheMap.provideParametricCoordinates().set()=group.getIntegrationPointsMatrix();
dataCacheDouble &sourceTerm = cacheMap.get(f);
fullMatrix<double> source;
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
diff --git a/Solver/dgFunctionIntegrator.cpp b/Solver/dgFunctionIntegrator.cpp
index 130a7fab5a342673c50b2954a00a9366b9f2f776..71b7559b91c1a66f27860a78ad5aa42122ad5d23 100644
--- a/Solver/dgFunctionIntegrator.cpp
+++ b/Solver/dgFunctionIntegrator.cpp
@@ -20,13 +20,13 @@ void dgFunctionIntegrator::compute(dgDofContainer *sol,fullMatrix<double> &resul
dgGroupOfElements &group=*sol->getGroups()->getElementGroup(iGroup);
fullMatrix<double> &solProxy=sol->getGroupProxy(&group);
cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
- UVW.set(group.getIntegrationPointsMatrix());
+ UVW.set()=group.getIntegrationPointsMatrix();
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
group.getCollocationMatrix().mult(solProxy , solutionQP);
fullMatrix<double> IPMatrix = group.getIntegrationPointsMatrix();
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
cacheMap.setElement(group.getElement(iElement));
- solutionQPe.setAsProxy(solutionQP, iElement*nbFields, nbFields );
+ solutionQPe.set().setAsProxy(solutionQP, iElement*nbFields, nbFields );
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iElement, iPt);
for (int k=0;k<nbRowResult;k++){
diff --git a/Solver/dgLimiter.cpp b/Solver/dgLimiter.cpp
index fab0d9d96e965257c212b3fb5d598ea7c7b489c0..c8cb8461fe0b69948458f204dd9a0dfff787ced7 100644
--- a/Solver/dgLimiter.cpp
+++ b/Solver/dgLimiter.cpp
@@ -118,8 +118,8 @@ int dgSlopeLimiter::apply ( dgDofContainer *solution)
dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
if (solutionEClipped){
for (int iElement=0 ; iElement<egroup->getNbElements() ;++iElement) {
- solutionE.setAsProxy(solGroup, iElement*nbFields, nbFields );
- solutionE.set((*solutionEClipped)());
+ solutionE.set().setAsProxy(solGroup, iElement*nbFields, nbFields );
+ solutionE.set()=(*solutionEClipped)();
}
}
}
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index 3b147ec38a01aadf4c5063333610c20b252022f8..71d15e8ae228329ace8d7f9f9303fb12a9b76d7a 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -28,7 +28,7 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
{
residual.scale(0);
_cacheMap->setNbEvaluationPoints(group.getNbIntegrationPoints());
- _UVW.set(group.getIntegrationPointsMatrix());
+ _UVW.set()=group.getIntegrationPointsMatrix();
// ----- 1 ---- get the solution at quadrature points
// ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints)
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
@@ -53,7 +53,7 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
// ----- 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.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
+ _solutionQPe.set().setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
if(_gradientSolutionQPe.somethingDependOnMe()){
dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
@@ -112,7 +112,7 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
{
residual.scale(0);
_cacheMap->setNbEvaluationPoints(group.getNbIntegrationPoints());
- _UVW.set(group.getIntegrationPointsMatrix());
+ _UVW.set()=group.getIntegrationPointsMatrix();
// ----- 1 ---- get the solution at quadrature points
// ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints)
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
@@ -145,7 +145,7 @@ void dgResidualVolume::compute1GroupWithJacobian(dgGroupOfElements &group, fullM
// ----- 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.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
+ _solutionQPe.set().setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
a.setAsProxy(A, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
b.setAsProxy(B, iElement*_nbFields*_nbFields, _nbFields*_nbFields);
@@ -321,9 +321,9 @@ void dgResidualInterface::compute1Group ( //dofManager &dof, // the DOF manager
for (int i=0; i<nbConnections; i++) {
// B1 ) adjust the proxies for this element
caches[i].setElement(connections[i]->getElement(iFace));
- caches[i].getParametricCoordinates(NULL).setAsProxy(connections[i]->getIntegrationPointsOnElement(iFace));
- caches[i].getSolution(NULL).setAsProxy(solutionQP, (iFace*nbConnections+i)*_nbFields, _nbFields);
- caches[i].getNormals(NULL).setAsProxy(connections[i]->getNormals(), iFace*group.getNbIntegrationPoints(), group.getNbIntegrationPoints());
+ caches[i].getParametricCoordinates(NULL).set() = connections[i]->getIntegrationPointsOnElement(iFace);
+ caches[i].getSolution(NULL).set().setAsProxy(solutionQP, (iFace*nbConnections+i)*_nbFields, _nbFields);
+ caches[i].getNormals(NULL).set().setAsProxy(connections[i]->getNormals(), iFace*group.getNbIntegrationPoints(), group.getNbIntegrationPoints());
// B2 ) compute the gradient of the solution
if(caches[i].getSolutionGradient(NULL).somethingDependOnMe()) {
dofs.setAsProxy(*solutionOnElements[i], _nbFields*connections[i]->getElementId(iFace), _nbFields);
diff --git a/Solver/function.h b/Solver/function.h
index 37ff1c7553d2b5c9f407a45ffd7176336ba07674..ff25b7bf09bcf01c32f8c9aa9aa0c22281617d00 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -76,13 +76,6 @@ public :
inline bool somethingDependOnMe() {
return !_dependOnMe.empty();
}
- inline int howManyDependOnMe() {
- return _dependOnMe.size();
- }
- inline int howManyDoIDependOn() {
- return _iDependOn.size();
- }
-
std::vector<dataCacheDouble*> _dependencies;
std::vector<const fullMatrix<double>*> _depM;
@@ -92,6 +85,7 @@ public :
protected:
fullMatrix<double> _value;
// do the actual computation and put the result into _value
+ // still virtual because it is overrided by conservation law terms, as soon as conservation law terms will be regular functions, we will remove this
virtual void _eval()
{
for(int i=0;i<_dependencies.size(); i++)
@@ -100,42 +94,22 @@ public :
}
public:
//set the value (without computing it by _eval) and invalidate the dependencies
- inline void set(const fullMatrix<double> &mat) {
- if(_valid)
- _invalidateDependencies();
- _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(const fullMatrix<double> &mat, int cShift, int c) {
- if(_valid)
- _invalidateDependencies();
- _value.setAsProxy(mat,cShift,c);
- _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(const fullMatrix<double> &mat) {
- if(_valid)
- _invalidateDependencies();
- _value.setAsProxy(mat,0,mat.size2());
- _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
+ // take care if you use this to set a proxy you must ensure that the value pointed to are not modified
+ // without further call to set because the dependencies won't be invalidate
inline fullMatrix<double> &set() {
if(_valid)
_invalidateDependencies();
- _valid=true;
+ _valid = true;
return _value;
}
inline const double &operator () (int i, int j)
{
if(!_valid) {
_eval();
- _valid=true;
+ _valid = true;
}
return _value(i,j);
}
diff --git a/Solver/functionDerivator.cpp b/Solver/functionDerivator.cpp
index 88d0e93541396f3f732d21b4558198285992172a..ea9f7d6150855d8bf959768067cd50a7d7aa3433 100644
--- a/Solver/functionDerivator.cpp
+++ b/Solver/functionDerivator.cpp
@@ -9,10 +9,10 @@ void functionDerivator::compute() {
_xDx = _xRef;
for (int i=0;i<_fRef.size1();i++)
_xDx(i,j) += _epsilon;
- _x.set(_xDx);
+ _x.set()=_xDx;
for (int i=0; i<_fRef.size1(); i++)
for (int k=0; k<_fRef.size2(); k++)
_dfdx(i, k*_xRef.size2()+j) = (_f(i,k)-_fRef(i,k))/_epsilon;
}
- _x.set(_xRef);
+ _x.set()=_xRef;
};