clipping for emilie

Solver/dgConservationLawPerfectGas.cpp

@@ -443,26 +443,29 @@ public:
 
 class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
   dataCacheDouble /
+  double _presMin, _rhoMin;
 public:
-  clipToPhysics(dataCacheMap &cacheMap):
+  clipToPhysics(dataCacheMap &cacheMap, double presMin, double rhoMin):
     sol(cacheMap.get("Solution",this))
-  {};
+  {
+    _presMin=presMin;
+    _rhoMin=rhoMin;
+    _value=fullMatrix<double>(cacheMap.getNbEvaluationPoints(),4);
+  };
   void _eval () { 
-    double rhomin = 1.e-3;
-    double presmin= 1.e-3;
     const int nQP = sol().size1();      
     for (size_t k = 0 ; k < nQP; k++ ){
       _value(k,0) = sol(k,0);
       _value(k,1) = sol(k,1);
       _value(k,2) = sol(k,1);
       _value(k,3) = sol(k,3);
-      if (sol(k,0) < rhomin)  
-	_value(k,0) = rhomin;
+      if (sol(k,0) < _rhoMin)  
+	_value(k,0) = _rhoMin;
       double rhoV2 = sol(k,1)*sol(k,1)+sol(k,2)*sol(k,2);
       rhoV2 /= sol(k,0);
       const double p = (GAMMA-1)*(sol(k,3) - 0.5*rhoV2);
-      if (p < presmin) 
-	_value(k,3) = presmin / (GAMMA-1) +  0.5 *rhoV2 ; 
+      if (p < _presMin) 
+	_value(k,3) = _presMin / (GAMMA-1) +  0.5 *rhoV2 ; 
      }
   }
 };
@@ -634,7 +637,7 @@ dataCacheDouble *dgPerfectGasLaw2d::newSourceTerm (dataCacheMap &cacheMap) const
     return new source(cacheMap,_sourceFunctionName);    
 }
 dataCacheDouble *dgPerfectGasLaw2d::newClipToPhysics( dataCacheMap &cacheMap) const {
-  return new clipToPhysics(cacheMap);
+  return new clipToPhysics(cacheMap,1e-3,1e-3);
 }
 
 dgPerfectGasLaw2d::dgPerfectGasLaw2d() 

Solver/dgSlopeLimiter.cpp

@@ -92,26 +92,25 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
 
     }
   }  
-
+/*
   //  --- CLIPPING: check unphysical values
- //  fullMatrix<double> solutionQP (group->getNbIntegrationPoints(),group->getNbElements() * nbFields);
-//   group->getCollocationMatrix().mult(solleft, solutionQP);
-//   //dataCacheDouble &solutionQPe = cacheMap.provideData("Solution");
-//   //solutionQPe.set(fullMatrix<double>(group.getNbIntegrationPoints(),nbFields));
-//   //dataCacheElement &cacheElement = cacheMap.getElement();
-//   dataCacheMap cacheMap(group->getNbIntegrationPoints());
-//   dataCacheDouble *clipping = _claw->newClipToPhysics(cacheMap);
-//   fullMatrix<double> clippedSol;
-//    for (int iElement=0 ; iElement<group->getNbElements() ;++iElement) {
-//      clippedSol.setAsProxy(solutionQP, iElement*nbFields, nbFields );
-//      for (int iPt =0; iPt< group->getNbIntegrationPoints(); iPt++) {
-//        for (int k=0;k<nbFields;k++)
-// 	 clippedSol(iPt,k) =  (*clipping)(iPt,k);
-//      }
-//      //cacheElement.set(group.getElement(iElement));
-//    }
-//    solleft.gemm(group->getRedistributionMatrix(),solutionQP);
+  dataCacheMap cacheMap(group->getNbNodes());
+  dgGroupOfElements* group = eGroups[0];   
+  fullMatrix<double> &solGroup = solution.getGroupProxy(0);
+  dataCacheDouble &solutionE = cacheMap.provideData("Solution");
+  solutionE.set(fullMatrix<double>(group.getNbNodes(),nbFields));
+  dataCacheElement &cacheElement = cacheMap.getElement();
+  dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
+  for (int iElement=0 ; iElement<group->getNbElements() ;++iElement) {
+    cacheElement.set(group.getElement(iElement));
+    solutionE.setAsProxy(solGroup,iElement*nbFields,nbFields);
+    for (int iPt =0; iPt< group->getNbNodes(); iPt++) {
+      solutionE.set((*solutionEclipped)());
+    }
+  }
+  delete solutionEClipped;
 
+*/
   //#if 0
 //   double rhomin = 1.e-3;
 //   double presmin= 1.e-3;