diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 0d62971981b26f5f3c6611c84df203ce45721bbb..a7701284f5148266b0e3e1c1a97360c381d8bda7 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -19,6 +19,7 @@ set(SRC
dgConservationLawShallowWater2d.cpp
dgConservationLawWaveEquation.cpp
dgConservationLawPerfectGas.cpp
+ SlopeLimiter.cpp
function.cpp
functionDerivator.cpp
luaFunction.cpp
diff --git a/Solver/SlopeLimiter.cc b/Solver/SlopeLimiter.cc
deleted file mode 100644
index 81800c32cba824ada1fdd9035043c7b80f9d57d9..0000000000000000000000000000000000000000
--- a/Solver/SlopeLimiter.cc
+++ /dev/null
@@ -1,143 +0,0 @@
-#include "dgLimiter.h"
-#include "dgGroupOfElements.h"
-#include "dgSystemOfEquations.h"
-#include "dgConservationLaw.h"
-
-
-// //----------------------------------------------------------------------------------
-// bool SlopeLimiter::apply ( dgDofContainer &solution ,
-// std::vector<dgGroupOfElements*> &eGroups,
-// std::vector<dgGroupOfFaces*> &fGroups)
-// {
-
-// const DGIntMatrix &topo = domain->getFaceToRegion();
-// const int fSize = solution->FunctionSpaceSize();
-// const ClosedFunctionSpace* fsC = DGModelDataManager::instance().getFunctionSpace(problemName,0);
-
-// // //--- Correct unphysical values
-// // double Amin = 1.e-5;
-// // for (int iE=0;iE<solution->getNbElements();iE++){
-// // fullMatrix soljElem = solution->touchElementCoefs(iE);
-// // for (int k=0;k< soljElem.size1() ;k++){
-// // if (soljElem(k,0) < Amin) soljElem(k,0) = Amin;
-// // }
-// // }
-
-// //WARNING FOR ONLY 1 GROUP OF FACES
-
-// // first compute max and min of all fields for all stencils
-// //----------------------------------------------------------
-// fullMatrix MIN (solution->getNbElements(), solution->getNbFields());
-// fullMatrix MAX (solution->getNbElements(), solution->getNbFields());
-// MIN.set_all ( 1.e22);
-// MAX.set_all (-1.e22);
-
-// fullMatrix MEANL (solution->getNbElements(), solution->getNbFields());
-// fullMatrix MEANR (solution->getNbElements(), solution->getNbFields());
-// MEANL.set_all ( 0.01);
-// MEANR.set_all ( 0.01);
-
-// dgGroupOfElements* group = fGroups[0];
-// for(int iFace=0 ; iFace<group->getNbElements();iFace++)
-// {
-// const int iElementL = group->getElementLeft(iFace);
-// const int iElementR = group->getElementRight(iFace);
-// if (iElementR >= 0)
-// {
-// const fullMatrix TempR = fullMatrix(solution->seeElementCoefs (iElementR));
-// const fullMatrix TempL = fullMatrix(solution->seeElementCoefs (iElementL));
-// for (int k=0; k<nbFields; ++k)
-// {
-// double AVGL = 0;
-// double AVGR = 0;
-// for (int i=0; i<fSize; ++i)
-// {
-// AVGL += TempL(i,k);
-// AVGR += TempR(i,k);
-// }
-// AVGL /= (double) fSize;
-// AVGR /= (double) fSize;
-// MIN ( iElementL , k ) = std::min ( AVGR , MIN ( iElementL , k ) );
-// MAX ( iElementL , k ) = std::max ( AVGR , MAX ( iElementL , k ) );
-// MIN ( iElementR , k ) = std::min ( AVGL , MIN ( iElementR , k ) );
-// MAX ( iElementR , k ) = std::max ( AVGL , MAX ( iElementR , k ) );
-
-// MEANR( iElementL,k ) = AVGR;
-// MEANL( iElementR,k ) = AVGL;
-
-// }
-// }
-// else{
-// const fullMatrix TempL = fullMatrix(solution->seeElementCoefs (iElementL));
-// for (int k=0; k<nbFields; ++k){
-// for (int i=0; i<fSize; ++i){
-// MIN ( iElementL , k ) = std::min ( TempL(i,k) , MIN ( iElementL , k ) );
-// MAX ( iElementL , k ) = std::max ( TempL(i,k) , MAX ( iElementL , k ) );
-
-
-// double AVG = 0;
-// for (int i=0; i<fSize; ++i) AVG += TempL(i,k);
-
-// }
-// }
-// }
-// }
-// // some parallel should be done here in order to compute averages on other processors
-// //----------------------------------------------------------
-
-// //...
-
-// // then limit the solution
-// //----------------------------------------------------------
-
-// for (int iElement=0 ; iElement<solution->getNbElements(); ++iElement)
-// {
-// fullMatrix Temp = fullMatrix(solution->touchElementCoefs (iElement));
-// for (int k=0; k<nbFields; ++k)
-// {
-// double AVG = 0.;
-// double locMax = -1.e22;
-// double locMin = 1.e22;
-// double neighMax = MAX (iElement,k);
-// double neighMin = MIN (iElement,k);
-// for (int i=0; i<fSize; ++i)
-// {
-// AVG += Temp(i,k);
-// locMax = std::max (locMax, Temp (i,k));
-// locMin = std::min (locMin, Temp (i,k));
-// }
-// AVG /= (double) fSize;
-
-// //SLOPE LIMITING DG
-// //-------------------
-// for (int i=0; i<fSize; ++i)Temp(i,k) -= AVG;
-
-// double slopeLimiterValue = 1.0;
-// if (locMax != AVG && locMax > neighMax) slopeLimiterValue = (neighMax-AVG) / (locMax-AVG);
-// if (locMin != AVG && locMin < neighMin) slopeLimiterValue = std::min ( slopeLimiterValue , (AVG-neighMin) / (AVG-locMin) );
-// if (AVG < neighMin) slopeLimiterValue = 0;
-// if (AVG > neighMax) slopeLimiterValue = 0;
-
-// if (detectDISC(iElement) == 0.0) slopeLimiterValue = 1.0; // do not limit
-// //if (slopeLimiterValue != 1.0) printf("limit (%g) elem =%d \n", slopeLimiterValue, iElement);
-
-// for (int i=0; i<fSize; ++i) Temp(i,k) *= slopeLimiterValue;
-
-// for (int i=0; i<fSize; ++i) Temp(i,k) += AVG;
-
-// //MINMOD DG (IF CHANGE TO THIS ADD LIMITER LINE IN RUNGEKUTTANNOPERATOR.cc)
-// //--------------------------------
-// // if (detectDISC(iElement) != 0.0){
-// // const size_t nn = Temp.size1();
-// // double alpha = 0.5; //0.5 1.0 1.3 (alpha=0.5 => MUSC Schem VAN LEER)
-// // Temp(0,k) = AVG-minmod(AVG-Temp(0,k), alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
-// // Temp(nn-1,k) = AVG+minmod(Temp(nn-1,k)-AVG, alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
-// // }
-
-
-// }
-// }
-
-// return true;
-
-// }
diff --git a/Solver/SlopeLimiter.cpp b/Solver/SlopeLimiter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f8f3f56cb7be073b06fd8b2db4a6c6a33b63607c
--- /dev/null
+++ b/Solver/SlopeLimiter.cpp
@@ -0,0 +1,139 @@
+#include "dgLimiter.h"
+#include "dgGroupOfElements.h"
+#include "dgSystemOfEquations.h"
+
+
+// //----------------------------------------------------------------------------------
+// bool SlopeLimiter::apply ( dgDofContainer &solution,
+// std::vector<dgGroupOfElements*> &eGroups,
+// std::vector<dgGroupOfFaces*> &fGroups)
+// {
+
+// //WARNING FOR ONLY 1 GROUP OF FACES
+// dgGroupOfFaces* group = fGroups[0];
+// fullMatrix<double> SolLeft = solution._dataProxys[0];
+// fullMatrix<double> SolRight = solution._dataProxys[0];
+// int nbFields = solution->getNbFields();
+// int totNbElems = solution->getNbElements();
+
+// // first compute max and min of all fields for all stencils
+// //----------------------------------------------------------
+// fullMatrix<double> MIN (totNbElems,nbFields );
+// fullMatrix<double> MAX (totNbElems, nbFields);
+// MIN.set_all ( 1.e22);
+// MAX.set_all (-1.e22);
+
+// fullMatrix<double> MEANL (totNbElems, nbFields);
+// fullMatrix<double> MEANR (totNbElems, nbFields);
+// MEANL.set_all ( 0.01);
+// MEANR.set_all ( 0.01);
+
+// for(int iFace=0 ; iFace<group->getNbElements();iFace++)
+// {
+// int iElementL = group->getElementLeftId(iFace);
+// int iElementR = group->getElementRightId(iFace);
+// if (iElementR >= 0)
+// {
+// fullMatrix<double> TempL, TempR;
+// TempL.setAsProxy(SolLeft, nbFields*iELementL, nbFields );
+// TempR.setAsProxy(SolRight, nbFields*iELementR, nbFields );
+// printf("TempL size 1 = %d %d ", TempL.size1(), TempL.size2());
+// exit(1);
+
+// int fSize = TempL.size2();
+// for (int k=0; k< fSize; ++k)
+// {
+// double AVGL = 0;
+// double AVGR = 0;
+// for (int i=0; i<fSize; ++i)
+// {
+// AVGL += TempL(i,k);
+// AVGR += TempR(i,k);
+// }
+// AVGL /= (double) fSize;
+// AVGR /= (double) fSize;
+// MIN ( iElementL , k ) = std::min ( AVGR , MIN ( iElementL , k ) );
+// MAX ( iElementL , k ) = std::max ( AVGR , MAX ( iElementL , k ) );
+// MIN ( iElementR , k ) = std::min ( AVGL , MIN ( iElementR , k ) );
+// MAX ( iElementR , k ) = std::max ( AVGL , MAX ( iElementR , k ) );
+
+// MEANR( iElementL,k ) = AVGR;
+// MEANL( iElementR,k ) = AVGL;
+
+// }
+// }
+// else{
+// fullMatrix<double> TempL, ;
+// TempL.setAsProxy(SolLeft, nbFields*iELementL, nbFields );
+// for (int k=0; k<nbFields; ++k){
+// for (int i=0; i<fSize; ++i){
+// MIN ( iElementL , k ) = std::min ( TempL(i,k) , MIN ( iElementL , k ) );
+// MAX ( iElementL , k ) = std::max ( TempL(i,k) , MAX ( iElementL , k ) );
+
+
+// double AVG = 0;
+// for (int i=0; i<fSize; ++i) AVG += TempL(i,k);
+
+// }
+// }
+// }
+// }
+// // some parallel should be done here in order to compute averages on other processors
+// //----------------------------------------------------------
+
+// //...
+
+// // then limit the solution
+// //----------------------------------------------------------
+
+// // for (int iElement=0 ; iElement<totNbElems; ++iElement)
+// // {
+// // fullMatrix Temp = fullMatrix(solution->touchElementCoefs (iElement));
+// // for (int k=0; k<nbFields; ++k)
+// // {
+// // double AVG = 0.;
+// // double locMax = -1.e22;
+// // double locMin = 1.e22;
+// // double neighMax = MAX (iElement,k);
+// // double neighMin = MIN (iElement,k);
+// // for (int i=0; i<fSize; ++i)
+// // {
+// // AVG += Temp(i,k);
+// // locMax = std::max (locMax, Temp (i,k));
+// // locMin = std::min (locMin, Temp (i,k));
+// // }
+// // AVG /= (double) fSize;
+
+// // //SLOPE LIMITING DG
+// // //-------------------
+// // for (int i=0; i<fSize; ++i)Temp(i,k) -= AVG;
+
+// // double slopeLimiterValue = 1.0;
+// // if (locMax != AVG && locMax > neighMax) slopeLimiterValue = (neighMax-AVG) / (locMax-AVG);
+// // if (locMin != AVG && locMin < neighMin) slopeLimiterValue = std::min ( slopeLimiterValue , (AVG-neighMin) / (AVG-locMin) );
+// // if (AVG < neighMin) slopeLimiterValue = 0;
+// // if (AVG > neighMax) slopeLimiterValue = 0;
+
+// // if (detectDISC(iElement) == 0.0) slopeLimiterValue = 1.0; // do not limit
+// // //if (slopeLimiterValue != 1.0) printf("limit (%g) elem =%d \n", slopeLimiterValue, iElement);
+
+// // for (int i=0; i<fSize; ++i) Temp(i,k) *= slopeLimiterValue;
+
+// // for (int i=0; i<fSize; ++i) Temp(i,k) += AVG;
+
+// // //MINMOD DG (IF CHANGE TO THIS ADD LIMITER LINE IN RUNGEKUTTANNOPERATOR.cc)
+// // //--------------------------------
+// // // if (detectDISC(iElement) != 0.0){
+// // // const size_t nn = Temp.size1();
+// // // double alpha = 0.5; //0.5 1.0 1.3 (alpha=0.5 => MUSC Schem VAN LEER)
+// // // Temp(0,k) = AVG-minmod(AVG-Temp(0,k), alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
+// // // Temp(nn-1,k) = AVG+minmod(Temp(nn-1,k)-AVG, alpha*(AVG-MEANL(iElement,k)), alpha*(MEANR(iElement,k)-AVG));
+// // // }
+
+
+// // }
+// // }
+
+// return true;
+
+// }
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 4517596aafbef51ba83518c305cbd05bd4e473e7..e93d4164964467e15076f78e9f8ee3a1b68dbd9e 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -313,14 +313,14 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
}
}
Unp._data->axpy(*(K._data),a[j]);
- if (limiter) limiter->apply(Unp);
+ if (limiter) limiter->apply(Unp, eGroups, fGroups);
}
for (int i=0;i<sol._dataSize;i++){
// printf("tempSol[%d] = %g\n",i,(*tempSol._data)(i));
// memcp
(*sol._data)(i)=(*Unp._data)(i);
- //if (limiter) limiter->apply(sol);
+ //if (limiter) limiter->apply(sol, eGroups, fGroups);
}
}
diff --git a/Solver/dgLimiter.h b/Solver/dgLimiter.h
index 576f2ad7804cc833ea2b90d98caee0d2eca17cc6..d26aeccc183131f668b841a122d2a31857d3cad9 100644
--- a/Solver/dgLimiter.h
+++ b/Solver/dgLimiter.h
@@ -3,24 +3,23 @@
#include "fullMatrix.h"
#include <vector>
+struct dgDofContainer;
+class dgGroupOfElements;
+class dgGroupOfFaces;
-
-class dgLimiter
-{
+class dgLimiter{
public:
dgLimiter () {}
- virtual ~dgLimiter() {}
- virtual bool apply ( dgDofContainer &sol ) =0;
+ virtual bool apply ( dgDofContainer &sol, std::vector<dgGroupOfElements*> &eGroups,
+ std::vector<dgGroupOfFaces*> &fGroup ) = 0;
};
-class SlopeLimiter : public dgLimiter
-{
-
- public :
+class SlopeLimiter : public dgLimiter{
+public :
SlopeLimiter () : dgLimiter () {}
- ~SlopeLimiter() {}
- bool apply ( dgDofContainer &sol );
-
+ virtual bool apply ( dgDofContainer &solution,
+ std::vector<dgGroupOfElements*> &eGroups,
+ std::vector<dgGroupOfFaces*> &fGroup);
};
#endif
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 31fcb63c559cc55e375309a2556f8205cf7e2fa1..0b2776516d4326d0926afccb0c5369619127e0e6 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -5,7 +5,7 @@
#include "MElement.h"
#include "PView.h"
#include "PViewData.h"
-
+#include "dgLimiter.h"
class dgConservationLawL2Projection : public dgConservationLaw {
std::string _functionName;
@@ -77,7 +77,8 @@ void dgSystemOfEquations::setup(){
double dgSystemOfEquations::RK44(double dt){
- _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide);
+ //dgLimiter *sl = new SlopeLimiter();
+ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, NULL);
return _solution->_data->norm();
}