diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index c75451330afcf186c39bd641130a4928abef311d..62d1d285ecd1637a74fbdf26eb976f974c147392 100644
--- a/Common/LuaBindings.cpp
+++ b/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
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index c1b6cc76d378c7b5157cfca6b4770ebf59da2ab5..29656e1845d32d7d59c9094379158035a2513675 100644
--- a/Geo/GVertex.h
+++ b/Geo/GVertex.h
@@ -78,6 +78,7 @@ class GVertex : public GEntity
bool isOnSeam(const GFace *gf) const;
std::vector<MPoint*> points;
+
};
#endif
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 908221ba922b69aee7ffd2c8b2592569eb078817..bc12cecde6a78593e7fe7bdb9837b1a6db6ff6e4 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -396,3 +396,16 @@ bool reparamMeshVertexOnEdge(const MVertex *v, const GEdge *ge, double ¶m)
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>();
+}
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 7e4287106dd96c001c0f067692b3ed3b11346a46..6ac10df95c688caed533b236066147507c82b371 100644
--- a/Geo/MVertex.h
+++ b/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{
diff --git a/Solver/TESTCASES/ForwardFacingStep.lua b/Solver/TESTCASES/ForwardFacingStep.lua
index 9f2a69f7f90894e4ab224878a0cc1db9888a494e..4c9165ce23ec5ec32919fd52f94be32144104abe 100644
--- a/Solver/TESTCASES/ForwardFacingStep.lua
+++ b/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))
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 796cc258a9178ab90181a2d46f391ffe0d65375e..7477ee22ddfcc58f239d01b160a2c51a2a181c62 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/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;
}
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index 49a2aaafa52624b340309097b0a9f71cb51c4c5f..22e974fde1585e71426d24c757df928122fe275c 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/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);
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index 026b5e990e02274e8fa266bfbcc94ade6f51ae89..f4ce4cb54ba0c732525ebb2fc39c83f21c18a9ff 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/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;
-
+
}
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 8616b5bf532e10f64b007623d35caa1405a47ea4..09d0d3fd2bae95e6855b8b600843aa0f90a05341 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/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