fix bug in clipping

-This line, and those below, will be ignored--

M    Common/LuaBindings.cpp
M    Geo/MVertex.cpp
M    Geo/GVertex.h
M    Geo/MVertex.h
M    Solver/dgAlgorithm.cpp
M    Solver/TESTCASES/ForwardFacingStep.lua
M    Solver/dgSlopeLimiter.cpp
M    Solver/dgSystemOfEquations.cpp
M    Solver/dgConservationLawPerfectGas.cpp

Common/LuaBindings.cpp

@@ -3,6 +3,8 @@
 #include <iostream>
 #include <string>
 #include "Gmsh.h"
+#include "MVertex.h"
+#include "MElement.h"
 #include "Bindings.h"
 #include "dgSystemOfEquations.h"
 #include "luaFunction.h"
@@ -169,6 +171,7 @@ binding::binding(){
   function::registerBindings(this);
   functionLua::registerBindings(this);
   function::registerDefaultFunctions();
+  MVertex::registerBindings(this);
 }
 binding *binding::_instance=NULL;
 #endif

Geo/GVertex.h

@@ -78,6 +78,7 @@ class GVertex : public GEntity
   bool isOnSeam(const GFace *gf) const;
 
   std::vector<MPoint*> points;
+  
 };
 
 #endif

Geo/MVertex.cpp

@@ -396,3 +396,16 @@ bool reparamMeshVertexOnEdge(const MVertex *v, const GEdge *ge, double &param)
   if(param < 1.e6) return true;
   return false;
 }
+
+#include "Bindings.h"
+
+void MVertex::registerBindings(binding *b)
+{
+  classBinding *cb = b->addClass<MVertex>("MVertex");
+  methodBinding *cm;
+  cm = cb->addMethod("getNum",&MVertex::getNum);
+  //  cm = cb->addMethod("x",&MVertex::x);
+  //  cm = cb->addMethod("y",&MVertex::y);
+  //  cm = cb->addMethod("z",&MVertex::z);
+  cm = cb->setConstructor<MVertex,double,double,double>();
+}

Geo/MVertex.h

@@ -15,6 +15,7 @@ class GEntity;
 class GEdge;
 class GFace;
 class MVertex;
+class binding;
 
 class MVertexLessThanLexicographic{
  public:
@@ -108,6 +109,7 @@ class MVertex{
   void writeMESH(FILE *fp, double scalingFactor=1.0);
   void writeBDF(FILE *fp, int format=0, double scalingFactor=1.0);
   void writeDIFF(FILE *fp, bool binary, double scalingFactor=1.0);
+  static void registerBindings(binding *b);
 };
 
 class MEdgeVertex : public MVertex{

Solver/TESTCASES/ForwardFacingStep.lua

@@ -59,12 +59,15 @@ DG:exportSolution('output/solution_0')
 print'*** solve ***'
 CFL = 2
 
+local x = os.clock()
+
 for i=1,5000 do
     dt = CFL * DG:computeInvSpectralRadius();    
-    norm = DG:RK44_limiter(dt)
+--    norm = DG:RK44_limiter(dt)
+    norm = DG:ForwardEuler(dt)
     DG:limitSolution()
     if (i % 10 == 0) then 
-       print('*** ITER ***',i,norm,dt)
+       print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
     end
     if (i % 100 == 0) then 
        DG:exportSolution(string.format("output/solution-%06d", i)) 

Solver/dgAlgorithm.cpp

@@ -281,7 +281,7 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw,			// conservation l
    double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
    double b[4] = {0.,h/2.0,h/2.0,h};
    
-   if (orderRK = 1.0){
+   if (orderRK == 1){
      a[0] = h;
    }
 

Solver/dgConservationLawPerfectGas.cpp

@@ -446,7 +446,7 @@ class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
   double _presMin, _rhoMin;
 public:
   clipToPhysics(dataCacheMap &cacheMap, double presMin, double rhoMin):
-    sol(cacheMap.get("SolToClip",this))
+    sol(cacheMap.get("Solution",this))
   {
     _presMin=presMin;
     _rhoMin=rhoMin;
@@ -457,7 +457,7 @@ public:
     for (size_t k = 0 ; k < nDofs; k++ ){
       _value(k,0) = sol(k,0);
       _value(k,1) = sol(k,1);
-      _value(k,2) = sol(k,1);
+      _value(k,2) = sol(k,2);
       _value(k,3) = sol(k,3);
       if (sol(k,0) < _rhoMin){
 	//printf("CL: clip rho min =%g \n", _rhoMin);

Solver/dgSlopeLimiter.cpp

@@ -92,42 +92,29 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
 
     }
   }  
-
   //  --- CLIPPING: check unphysical values
   for (int iG = 0; iG < eGroups.size(); iG++){
     dgGroupOfElements* egroup = eGroups[iG];  
     fullMatrix<double> &solGroup = solution.getGroupProxy(iG);
 
     dataCacheMap cacheMap(egroup->getNbNodes());//nbdofs for each element
-
-    dataCacheDouble &solutionE = cacheMap.provideData("SolToClip");
-    solutionE.set(fullMatrix<double>(egroup->getNbNodes(),nbFields));//e.g 3x4 for perfect gas
-
-    dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
+    dataCacheDouble &solutionE = cacheMap.provideData("Solution");
     dataCacheElement &cacheElement = cacheMap.getElement();
-    for (int iElement=0 ; iElement<egroup->getNbElements() ;++iElement) {
-
-      cacheElement.set(egroup->getElement(iElement));
-      solutionE.setAsProxy(solGroup,iElement*nbFields,nbFields);
-      solutionE.set((*solutionEClipped)());
-
-//       fullMatrix<double> solElem;
-//       double min = 1.e-3;
-//       solElem.setAsProxy(solGroup, nbFields*iElement, nbFields );   
-//       for (int k=0;k< solElem.size1() ;k++){
-// 	if (solElem(k,0) < min) solElem(k,0) = min;
-// 	double rhoV2 = solElem(k,1)*solElem(k,1)+solElem(k,2)*solElem(k,2);
-// 	rhoV2 /= solElem(k,0);
-//        const double p = (1.4-1)*(solElem(k,3) - 0.5* rhoV2);
-//        if (p < 1.e-3)  solElem(k,3) = 0.5 *rhoV2 + min / (1.4-1);
-//       }
-      
+    dataCacheDouble *solutionEClipped = _claw->newClipToPhysics(cacheMap);
+    if (solutionEClipped){
+      for (int iElement=0 ; iElement<egroup->getNbElements() ;++iElement) {
+	solutionE.setAsProxy(solGroup, iElement*nbFields, nbFields );
+	fullMatrix<double> Temp;  
+	Temp.setAsProxy(solGroup, nbFields*iElement, nbFields );    	
+	cacheElement.set(egroup->getElement(iElement));
+	for (int K=0;K<Temp.size1();K++)
+	  for (int L=0;L<Temp.size2();L++)
+	    Temp(K,L) = (*solutionEClipped)(K,L);      
+      }
+      delete solutionEClipped;
     }
-    delete solutionEClipped;
-    
-  }
-
+  }  
   return true; 
-
+  
 }
 

Solver/dgSystemOfEquations.cpp

@@ -66,12 +66,20 @@ void dgSystemOfEquations::registerBindings(binding *b) {
   cm = cb->addMethod("setOrder",&dgSystemOfEquations::setOrder);
   cm->setArgNames("order",NULL);
   cm->setDescription("set the polynpolynomialomial order of the lagrange shape functions");
-  cb->addMethod("limitSolution",&dgSystemOfEquations::limitSolution);
-  cb->addMethod("computeInvSpectralRadius",&dgSystemOfEquations::computeInvSpectralRadius);
-  cb->addMethod("RK44_limiter",&dgSystemOfEquations::RK44_limiter);
-  cb->addMethod("multirateRK43",&dgSystemOfEquations::multirateRK43);
-  cb->addMethod("saveSolution",&dgSystemOfEquations::saveSolution);
-  cb->addMethod("loadSolution",&dgSystemOfEquations::loadSolution);
+  cm = cb->addMethod("limitSolution",&dgSystemOfEquations::limitSolution);
+  cm->setDescription("apply a slope limiter to the solution (only if polynomial order p = 1 for now).");
+  cm = cb->addMethod("computeInvSpectralRadius",&dgSystemOfEquations::computeInvSpectralRadius);
+  cm->setDescription("returns the inverse of the spectral radius (largest eigenvalue) of L(u). Useful for computing stable explicit time step");
+  cm = cb->addMethod("RK44_limiter",&dgSystemOfEquations::RK44_limiter);
+  cm->setArgNames("dt",NULL);
+  cm->setDescription("do one RK44 time step with the slope limiter (only for p=1)");
+  cm = cb->addMethod("multirateRK43",&dgSystemOfEquations::multirateRK43);
+  cm = cb->addMethod("saveSolution",&dgSystemOfEquations::saveSolution);
+  cm->setArgNames("filename",NULL);
+  cm->setDescription("dump the solution in binary format");
+  cm = cb->addMethod("loadSolution",&dgSystemOfEquations::loadSolution);
+  cm->setArgNames("filename",NULL);
+  cm->setDescription("reload a solution in binary format");
 }
 
 // do a L2 projection