No commit message

No commit message

Common/GmshDefines.h

@@ -97,7 +97,6 @@
 #define MSH_QUA_16I 40
 #define MSH_QUA_20 41
 
-
 #define MSH_NUM_TYPE 35
 
 // Geometric entities

Geo/MLine.cpp

@@ -2,7 +2,7 @@
 //
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
-
+#include "GmshDefines.h"
 #include "MLine.h"
 #include "GaussLegendre1D.h"
 #include "Context.h"

Geo/MTriangle.cpp

@@ -6,6 +6,7 @@
 #include "MTriangle.h"
 #include "Numeric.h"
 #include "Context.h"
+#include "GmshDefines.h"
 
 #if defined(HAVE_MESH)
 #include "qualityMeasures.h"

Numeric/polynomialBasis.cpp

@@ -741,11 +741,12 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
 static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
 {
 
+  closure.clear();
   for (int iRotate = 0; iRotate < 3; iRotate++){
     for (int iSign = 1; iSign >= -1; iSign -= 2){
       for (int iFace = 0; iFace < 4; iFace++){
 	std::vector<int> closure_face;
-	getFaceClosure(iFace, iSign, iRotate, closure_face, order);
+	getFaceClosure(iFace, iSign, iRotate, closure_face, order);	
 	closure.push_back(closure_face);
       }
     }
@@ -768,24 +769,23 @@ static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, i
   }
 }
 
-
-
 static void generate1dVertexClosure(polynomialBasis::clCont &closure)
 {
+
   closure.clear();
   closure.resize(2);
   closure[0].push_back(0);
   closure[1].push_back(1);
-  }
+  
+}
 std::map<int, polynomialBasis> polynomialBases::fs;
 
 const polynomialBasis &polynomialBases::find(int tag) 
 {
   std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
-  if (it != fs.end()) return it->second;
-  
+  if (it != fs.end())     return it->second;
   polynomialBasis F;
-  
+
   switch (tag){
   case MSH_PNT:
     F.monomials = generate1DMonomials(0);
@@ -948,6 +948,7 @@ const polynomialBasis &polynomialBases::find(int tag)
   return fs[tag];
 }
 
+
 std::map<std::pair<int, int>, fullMatrix<double> > polynomialBases::injector;
 
 const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)

Solver/TESTCASES/Advection1D.lua

@@ -4,8 +4,7 @@ model:load ('edge.geo')
 -- model:save ('edge.msh')
 model:load ('edge.msh')
 dg = dgSystemOfEquations (model)
-dg:setOrder(1)
-
+dg:setOrder(0)
 
 -- conservation law
 -- advection speed
@@ -45,16 +44,17 @@ function initial_condition( xyz , f )
   end
 end
 dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
+print'***exporting init solution ***'
 
 dg:exportSolution('output/Adv1D_unlimited')
 dg:limitSolution()
-dg:exportSolution('output/Adv1D_00000')
+dg:exportSolution('output/Adv1D-00000')
 
 
 -- main loop
 n = 5
 for i=1,100*n do
-  norm = dg:RK44_limiter(0.03)
+  norm = dg:RK44(0.03)
   if (i % n == 0) then 
     print('iter',i,norm)
     dg:exportSolution(string.format("output/Adv1D-%05d", i)) 

Solver/TESTCASES/Advection3D.lua

@@ -21,9 +21,9 @@ end
 -- conservation law
 -- advection speed
 v=fullMatrix(3,1);
-v:set(0,0,0.15)
+v:set(0,0,0)
 v:set(1,0,0)
-v:set(2,0,0)
+v:set(2,0,0.15)
 
 law = dgConservationLawAdvection(functionConstant(v):getName(),'')
 dg:setConservationLaw(law)
@@ -49,7 +49,7 @@ print'***exporting init solution ***'
 -- main loop
 n = 5
 for i=1,100*n do
-  norm = dg:RK44(0.0003)
+  norm = dg:RK44(0.03)
   if (i % n == 0) then 
     print('iter',i,norm)
     dg:exportSolution(string.format("output/Adv3D-%05d", i)) 

Solver/dgAlgorithm.cpp

@@ -617,35 +617,35 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
   //  residu[0]->print("Volume");
   //interface term
   for(size_t i=0;i<fGroups.size() ; i++) {
-    dgGroupOfFaces &faces = *fGroups[i];
-    int iGroupLeft = -1, iGroupRight = -1;
-    for(size_t j=0;j<eGroups.size() ; j++) {
-      if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
-      if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
+      dgGroupOfFaces &faces = *fGroups[i];
+      int iGroupLeft = -1, iGroupRight = -1;
+      for(size_t j=0;j<eGroups.size() ; j++) {
+	if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
+	if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
+      }
+      fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
+      fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
+      faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
+      residualInterface(claw,faces,solInterface,*solution[iGroupLeft], *solution[iGroupRight],residuInterface);
+      faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
     }
-    
-    fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
-    fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
-    faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
-    residualInterface(claw,faces,solInterface,*solution[iGroupLeft], *solution[iGroupRight],residuInterface);
-    faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
-  }
-  //  residu[0]->print("Interfaces");
-  //boundaries
-  for(size_t i=0;i<bGroups.size() ; i++) {
-    dgGroupOfFaces &faces = *bGroups[i];
-    int iGroupLeft = -1, iGroupRight = -1;
-    for(size_t j=0;j<eGroups.size() ; j++) {
-      if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
-      if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
+    //  residu[0]->print("Interfaces");
+    //boundaries
+    for(size_t i=0;i<bGroups.size() ; i++) {
+      dgGroupOfFaces &faces = *bGroups[i];
+      int iGroupLeft = -1, iGroupRight = -1;
+      for(size_t j=0;j<eGroups.size() ; j++) {
+	if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
+	if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
+      }
+      fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
+      fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
+      faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
+      residualBoundary(claw,faces,solInterface,*solution[iGroupLeft],residuInterface);
+      faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
     }
-    fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
-    fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
-    faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
-    residualBoundary(claw,faces,solInterface,*solution[iGroupLeft],residuInterface);
-    faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
-  }
-  //  residu[0]->print("Boundaries");
+
+    //  residu[0]->print("Boundaries");
 }
 
 void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)

Solver/dgConservationLaw.cpp

@@ -39,6 +39,34 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
   }
 };
 
+class dgBoundarySymmetry : public dgBoundaryCondition {
+  dgConservationLaw &_claw;
+  class term : public dataCacheDouble {
+    dataCacheDouble *riemannSolver;
+    dgConservationLaw &_claw;
+    public:
+    term(dgConservationLaw &claw, dataCacheMap &cacheMapLeft):
+      dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw.nbFields()), _claw(claw)
+    {
+      riemannSolver=_claw.newRiemannSolver(cacheMapLeft,cacheMapLeft);
+      riemannSolver->addMeAsDependencyOf(this);
+    }
+
+    void _eval() {
+      if(riemannSolver){
+        for(int i=0;i<_value.size1(); i++)
+          for(int j=0;j<_value.size2(); j++)
+            _value(i,j) = (*riemannSolver)(i,j);
+      }
+    }
+  };
+  public:
+  dgBoundarySymmetry(dgConservationLaw &claw): _claw(claw) {}
+  dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+    return new term(_claw,cacheMapLeft);
+  }
+};
+
 class dgBoundaryCondition0Flux : public dgBoundaryCondition {
   dgConservationLaw &_claw;
   class term : public dataCacheDouble {
@@ -55,6 +83,9 @@ class dgBoundaryCondition0Flux : public dgBoundaryCondition {
   }
 };
 
+dgBoundaryCondition *dgConservationLaw::newSymmetryBoundary() {
+  return new dgBoundarySymmetry(*this);
+}
 dgBoundaryCondition *dgConservationLaw::newOutsideValueBoundary(const std::string outsideValueFunctionName) {
   return new dgBoundaryConditionOutsideValue(*this,outsideValueFunctionName);
 }
@@ -68,6 +99,7 @@ void dgConservationLaw::registerBindings(binding *b){
   classBinding *cb = b->addClass<dgConservationLaw>("dgConservationLaw");
   cb->addMethod("addBoundaryCondition",&dgConservationLaw::addBoundaryCondition);
   cb->addMethod("new0FluxBoundary",&dgConservationLaw::new0FluxBoundary);
+  cb->addMethod("newSymmetryBoundary",&dgConservationLaw::newSymmetryBoundary);
   cb->addMethod("newOutsideValueBoundary",&dgConservationLaw::newOutsideValueBoundary);
 }
 

Solver/dgConservationLaw.h

@@ -53,6 +53,7 @@ class dgConservationLaw {
   //a generic boundary condition using the Riemann solver of the conservation Law
   dgBoundaryCondition *newOutsideValueBoundary(std::string outsideValueFunctionName);
   dgBoundaryCondition *new0FluxBoundary();
+  dgBoundaryCondition *newSymmetryBoundary();
 
   static void registerBindings(binding *b);
 };

Solver/dgGroupOfElements.cpp

@@ -229,6 +229,7 @@ void dgGroupOfFaces::init(int pOrder) {
   for (size_t i=0; i<_closuresRight.size(); i++)
     _integrationPointsRight.push_back(dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsRight,_closuresRight[i]));
   double f[256];
+  
   for (int j=0;j<_integration->size1();j++) {
     _fsFace->f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
     const double weight = (*_integration)(j,3);
@@ -237,6 +238,7 @@ void dgGroupOfFaces::init(int pOrder) {
       (*_collocation)(j,k) = f[k];
     }
   }
+  
   for (int i=0;i<_faces.size();i++){
     MElement *f = _faces[i];
     for (int j=0;j< _integration->size1() ; j++ ){
@@ -245,7 +247,7 @@ void dgGroupOfFaces::init(int pOrder) {
       (*_interfaceSurface)(i,0) += (*_integration)(j,3)*(*_detJac)(j,i);
     }
   }
-
+  
   // compute data on quadrature points : normals and dPsidX
   double g[256][3];
   _normals = new fullMatrix<double> (3,_integration->size1()*_faces.size());
@@ -255,58 +257,59 @@ void dgGroupOfFaces::init(int pOrder) {
   }
   int index = 0;
   for (size_t i=0; i<_faces.size();i++){
-    const std::vector<int> &closureLeft=getClosureLeft(i);
-    fullMatrix<double> &intLeft=_integrationPointsLeft[_closuresIdLeft[i]];
-    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++){
-        int inode=closureLeft[k];
-        nu += g[inode][0];
-        nv += g[inode][1];
-        nw += g[inode][2];
+  
+      const std::vector<int> &closureLeft=getClosureLeft(i);
+      fullMatrix<double> &intLeft=_integrationPointsLeft[_closuresIdLeft[i]];
+      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++){
+	  int inode=closureLeft[k];
+	  nu += g[inode][0];
+	  nv += g[inode][1];
+	  nw += g[inode][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++;
       }
-      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++;
-    }
-    // there is nothing on the right for boundary groups
-    if(_fsRight){
-      fullMatrix<double> &intRight=_integrationPointsRight[_closuresIdRight[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];
-        }
+      // there is nothing on the right for boundary groups
+      if(_fsRight){
+	fullMatrix<double> &intRight=_integrationPointsRight[_closuresIdRight[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];
+	  }
+	}
       }
-    }
   }
 }
 
@@ -399,7 +402,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
           if(vertexMap.find(vertex) == vertexMap.end()){
             vertexMap[vertex] = i;
           }else{
-            addVertex(vertex,vertexMap[vertex],i);
+	    addVertex(vertex,vertexMap[vertex],i);
           }
         }
       }
@@ -454,7 +457,6 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroup
       std::map<MVertex*,int> vertexMap1;
       _closuresLeft = _fsLeft->vertexClosure;
       _closuresRight = _fsRight->vertexClosure;
-
       for(int i=0; i<elGroup1.getNbElements(); i++){
         MElement &el = *elGroup1.getElement(i);
         for (int j=0; j<el.getNumVertices(); j++){
@@ -597,7 +599,7 @@ void dgGroupOfFaces::mapFromInterface ( int nFields,
         for(size_t j =0; j < closureLeft.size(); j++)
           vLeft(closureLeft[j], _left[i]*nFields + iField) += v(j, i*2*nFields + iField);
         for(size_t j =0; j < closureRight.size(); j++)
-          vRight(closureRight[j], _right[i]*nFields + iField) += v(j, (i*2+1)*nFields + iField);
+	  vRight(closureRight[j], _right[i]*nFields + iField) += v(j, (i*2+1)*nFields + iField);
       }
     }
   }

Solver/dgGroupOfElements.h

@@ -148,7 +148,7 @@ public:
   inline double getElementVolumeLeft(int iFace) const {return _groupLeft.getElementVolume(_left[iFace]);}
   inline double getElementVolumeRight(int iFace) const {return _groupRight.getElementVolume(_right[iFace]);}
   inline MElement* getFace (int iElement) const {return _faces[iElement];}  
-  inline const std::vector<int> &getClosureLeft(int iFace) const{ return _closuresLeft[_closuresIdLeft[iFace]];}
+  inline const std::vector<int> &getClosureLeft(int iFace) const{ return _closuresLeft[_closuresIdLeft[iFace]]; }
   inline const std::vector<int> &getClosureRight(int iFace) const{ return _closuresRight[_closuresIdRight[iFace]];}
   inline fullMatrix<double> &getIntegrationOnElementLeft(int iFace) { return _integrationPointsLeft[_closuresIdLeft[iFace]];}
   inline fullMatrix<double> &getIntegrationOnElementRight(int iFace) { return _integrationPointsRight[_closuresIdRight[iFace]];}

Solver/dgSlopeLimiter.cpp

@@ -9,7 +9,8 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
 			     std::vector<dgGroupOfFaces*> &fGroups) 
 {    
 
-  //WARNING FOR ONLY 1 GROUP OF FACES    
+  //WARNING: ONLY FOR 1 GROUP OF FACES 
+  //TODO: make this more general   
   dgGroupOfFaces* group = fGroups[0];   
   fullMatrix<double> &SolLeft = *(solution._dataProxys[0]);
   fullMatrix<double> &SolRight = *(solution._dataProxys[0]);    
@@ -104,7 +105,6 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
         locMin = std::min (locMin, Temp (i,k)); 
       }	
       AVG /= (double) fSize;  
-      // printf("AVG %e   LocMax %e   Locmin %e\n",AVG,locMax,locMin);
 
       //SLOPE LIMITING DG
       //-------------------  
@@ -116,24 +116,10 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
       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 (%e) elem =%d \n", slopeLimiterValue, iElement);
-//      slopeLimiterValue=0;
-
       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;