diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index f292d494b5e14569f672e28bfd5f51a59c4707f7..c41a265c848fa0fd7d8ca5f2825af209385a48b9 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -171,16 +171,16 @@ void CreateOutputFile(std::string fileName, int format)
break;
case FORMAT_MSH:
+ GModel::current()->writeMSH
+ (fileName, CTX::instance()->mesh.mshFileVersion,
+ CTX::instance()->mesh.binary, CTX::instance()->mesh.saveAll,
+ CTX::instance()->mesh.saveParametric, CTX::instance()->mesh.scalingFactor);
if (CTX::instance()->mesh.saveDistance){
GModel::current()->writeDistanceMSH
(fileName, CTX::instance()->mesh.mshFileVersion,
CTX::instance()->mesh.binary, CTX::instance()->mesh.saveAll,
CTX::instance()->mesh.saveParametric, CTX::instance()->mesh.scalingFactor);
}
- GModel::current()->writeMSH
- (fileName, CTX::instance()->mesh.mshFileVersion,
- CTX::instance()->mesh.binary, CTX::instance()->mesh.saveAll,
- CTX::instance()->mesh.saveParametric, CTX::instance()->mesh.scalingFactor);
break;
case FORMAT_STL:
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index c7ec58263bfb12fc99bb4491269f4c078ffd12ff..15c512c124a2bfd69400c91872ffdbd4a6360eea 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -677,15 +677,6 @@ static void getFaceClosure(int iFace, int iSign, int iRotate, std::vector<int> &
default:
int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
int order1node[4][3] = {{0, 2, 1}, {0, 1, 3}, {0, 3, 2}, {3, 1, 2}};
- // int facedofstmp[200];
- // face 0 | 0 1 2
- // face 1 | 0 4 -3
- // face 2 | 1 5 -4
- // face 3 | -2 5 -3
- // edge 0 | 4 -> 4+order-2
- // edge 1 | 4+order-1 -> 4 + 2*(order-1)-1
- // edge 2 | 4+2*(order-1)-> 4+ 3*(order-1)-1
- // edge k | 4+k*(order-1) -> 4+(k+1)*(order-1)-1
for (int i = 0; i < 3; ++i){
int k = (3 + (iSign * i) + iRotate) % 3; //- iSign * iRotate
closure[i] = order1node[iFace][k];
diff --git a/Solver/TESTCASES/ForwardFacingStep.lua b/Solver/TESTCASES/ForwardFacingStep.lua
index d0ec53d98bd660d29ed094800943b47d89c5fb12..9f2a69f7f90894e4ab224878a0cc1db9888a494e 100644
--- a/Solver/TESTCASES/ForwardFacingStep.lua
+++ b/Solver/TESTCASES/ForwardFacingStep.lua
@@ -1,4 +1,4 @@
-MACH = 3;
+MACH = 3.0;
GAMMA = 1.4;
U = 3.0
V = 0.0
@@ -26,6 +26,7 @@ end
--[[
Example of a lua program driving the DG code
--]]
+
order = 1
print'*** Loading the mesh and the model ***'
myModel = GModel ()
@@ -56,7 +57,6 @@ print'*** export ***'
DG:exportSolution('output/solution_0')
print'*** solve ***'
-
CFL = 2
for i=1,5000 do
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 8d3dc1bac928221d85011ed5356f4c301e989e8c..796cc258a9178ab90181a2d46f391ffe0d65375e 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -268,7 +268,7 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
std::vector<dgGroupOfElements*> &eGroups, // group of elements
std::vector<dgGroupOfFaces*> &fGroups, // group of interfacs
std::vector<dgGroupOfFaces*> &bGroups, // group of boundaries
- double h, // time-step
+ double h, // time-step
dgDofContainer &sol,
dgDofContainer &resd,
dgLimiter *limiter,
@@ -305,9 +305,6 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
if(j){
K._data->scale(b[j]);
K._data->axpy(*(sol._data));
- if (limiter){
- // limiter->apply(K, eGroups, fGroups);
- }
}
if (limiter){
@@ -325,7 +322,6 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
}
}
Unp._data->axpy(*(K._data),a[j]);
- //if (limiter) limiter->apply(Unp, eGroups, fGroups);
}
if (limiter) limiter->apply(Unp, eGroups, fGroups);
for (int i=0;i<sol._dataSize;i++){
@@ -546,7 +542,6 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
}
}
-
}
// ----- 3 ---- do the redistribution at face nodes using BLAS3
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
diff --git a/Solver/dgConservationLawAdvection.h b/Solver/dgConservationLawAdvection.h
index f56894bb21ca4e7348df329a6bfd7b0c7d2fbcb7..1251f992b0f9c922d96992344047290116b6549f 100644
--- a/Solver/dgConservationLawAdvection.h
+++ b/Solver/dgConservationLawAdvection.h
@@ -1,6 +1,8 @@
#ifndef _DG_CONSERVATION_LAW_ADVECTION_H
#define _DG_CONSERVATION_LAW_ADVECTION_H
#include "dgConservationLaw.h"
+//Advection diffusion equation
+//dc/dt + v div(c) - nu lapl(c) = 0
class dgConservationLawAdvection : public dgConservationLaw {
std::string _vFunctionName,_nuFunctionName;
class advection;
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index d3cecab9427e3cf63998176841d4fa1ff8132e28..cf95a2987f08ddb736cea089ac220232ec0e866a 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -213,11 +213,14 @@ static inline void _ROE2D (const double &_GAMMA,
const double *solR,
double *FLUX){
+ //sol cons var: rho, rhou, rhov, rhoE
+ //u prim var : rho, u, v, p
const double gamma = _GAMMA;
const double GM1 = gamma - 1.;
double uL[4], uR[4];
- const double overRhoL = 1./solL[0];
+ const double overRhoL = 1./solL[0];
+
uL[0] = solL[0];
uL[1] = solL[1]*overRhoL;
uL[2] = solL[2]*overRhoL;
@@ -228,38 +231,33 @@ static inline void _ROE2D (const double &_GAMMA,
uR[1] = solR[1]*overRhoR;
uR[2] = solR[2]*overRhoR;
uR[3] = GM1*(solR[3] - 0.5* (solR[1]*solR[1]+solR[2]*solR[2])*overRhoR);
- // central contributions
-
- double halfrhoun;
-
- /* --- left contributions ---*/
+
+ //* --- central contributions ---*/
+
+ double halfrhoun;
+ // --- left contributions
halfrhoun = 0.5*(uL[0]*(uL[1]*nx + uL[2]*ny));
- double HL = gamma/GM1* uL[3]/uL[0]+0.5*(uL[1]*uL[1]+
- uL[2]*uL[2]);
-
+ double HL = gamma/GM1* uL[3]/uL[0]+0.5*(uL[1]*uL[1]+uL[2]*uL[2]);
FLUX[0] = halfrhoun;
FLUX[1] = halfrhoun*uL[1] + .5*uL[3]*nx;
FLUX[2] = halfrhoun*uL[2] + .5*uL[3]*ny;
FLUX[3] = halfrhoun*HL;
- /* --- right contributions ---*/
-
+ // --- right contributions
halfrhoun = 0.5*(uR[0]*(uR[1]*nx+uR[2]*ny));
- double HR = gamma/GM1* uR[3]/uR[0]+0.5*(uR[1]*uR[1]+uR[2]*uR[2]);
-
+ double HR = gamma/GM1* uR[3]/uR[0]+0.5*(uR[1]*uR[1]+uR[2]*uR[2]);
FLUX[0] += halfrhoun;
FLUX[1] += halfrhoun*uR[1] + .5*uR[3]*nx;
FLUX[2] += halfrhoun*uR[2] + .5*uR[3]*ny;
FLUX[3] += halfrhoun*HR;
- /* --- add dissipation ---*/
+ /* --- add rhoe dissipation ---*/
double sqr_rhoL = sqrt(uL[0]);
- double sqr_rhoR = sqrt(uR[0]);
-
+ double sqr_rhoR = sqrt(uR[0]);
double invz1 = 1./ (sqr_rhoL + sqr_rhoR);
- // double rho = sqr_rhoL * sqr_rhoR;
+ //rhoe average state
double u = ( sqr_rhoL* uL[1] + sqr_rhoR * uR[1] ) * invz1;
double v = ( sqr_rhoL* uL[2] + sqr_rhoR * uR[2] ) * invz1;
double H = ( sqr_rhoL* HL + sqr_rhoR * HR ) * invz1;
@@ -281,17 +279,14 @@ static inline void _ROE2D (const double &_GAMMA,
double g1OnC2 = g1*oC2;
double TtOnT = (1.0 - u2*g1OnC2);
- // matrix of left eigenvectors
-
+ // matrix of left eigenvectors
double L[16] = {
nx*TtOnT ,nx*u*g1OnC2 ,nx*v*g1OnC2 , -nx*g1OnC2, // L1
- tet , ny ,-nx , 0, // L3
g1*u2 - c*un , c*nx - g1*u , c*ny - g1*v, g1, // L3
g1*u2 + c*un ,-c*nx - g1*u ,-c*ny - g1*v, g1}; // L4
-
// characteristic decomposition of differences
-
double dW[4] = {0,0,0,0};
int idx = 0;
for (int i=0;i<4;i++)
@@ -299,7 +294,6 @@ static inline void _ROE2D (const double &_GAMMA,
dW[i] += L[idx++]*dU[j];
// matrix of right eigenvectors
-
double R[16] = {
//R1 //R2 //R3 //R4
nx ,0 ,0.5*oC2 ,0.5*oC2,
@@ -307,11 +301,9 @@ static inline void _ROE2D (const double &_GAMMA,
v*nx ,- nx ,0.5*(ny*oC + v*oC2),0.5*(-ny*oC + v*oC2),
u2*nx,tet ,0.5*(un*oC + H*oC2),0.5*(-un*oC + H*oC2)};
-
// eigenvalues
- // KH : shouldn't we take into account an entropy correction ?
- // absorb half the surface : scaling wrt central term
-
+ // TODO : shouldn't we take into account an entropy correction ?
+ // absorb half the surface : scaling wrt central term
const double A = 0.5;
double eps = 1.e-6;
@@ -330,8 +322,7 @@ static inline void _ROE2D (const double &_GAMMA,
lflux -= lA[j]*dW[j]*R[index++];
FLUX[k] = -lflux;
}
-}
-// perfect gas law, GAMMA is the only parameter
+};
class dgPerfectGasLaw2d::advection : public dataCacheDouble {
dataCacheDouble /
@@ -345,7 +336,6 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
_value=fullMatrix<double>(nQP,8);
const double GM1 = GAMMA - 1.0;
for (size_t k = 0 ; k < nQP; k++ ){
- // printf("%d %g %g %g %g\n",k,sol(k,0),sol(k,1),sol(k,2),sol(k,3));
const double invrho = 1./sol(k,0);
const double q12 = sol(k,1)*sol(k,2)*invrho;
@@ -365,10 +355,6 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
_value(k,2+4) = q22+p;
_value(k,3+4) = sol(k,2)*qq;
- /* _value(k,8) = 0;
- _value(k,9) = 0;
- _value(k,10) = 0;
- _value(k,11) = 0;*/
}
}
};
@@ -455,6 +441,32 @@ public:
}
};
+class dgPerfectGasLaw2d::clipToPhysics : public dataCacheDouble {
+ dataCacheDouble /
+public:
+ clipToPhysics(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this))
+ {};
+ 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;
+ 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 ;
+ }
+ }
+};
+
class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
@@ -621,6 +633,9 @@ dataCacheDouble *dgPerfectGasLaw2d::newSourceTerm (dataCacheMap &cacheMap) const
else
return new source(cacheMap,_sourceFunctionName);
}
+dataCacheDouble *dgPerfectGasLaw2d::newClipToPhysics( dataCacheMap &cacheMap) const {
+ return new clipToPhysics(cacheMap);
+}
dgPerfectGasLaw2d::dgPerfectGasLaw2d()
{
diff --git a/Solver/dgConservationLawPerfectGas.h b/Solver/dgConservationLawPerfectGas.h
index 8eb66fd4b01e880027c209afbed128cf6696ad20..3cbe06c1c44f0dfd55e7f596e2c006fb78f255fe 100644
--- a/Solver/dgConservationLawPerfectGas.h
+++ b/Solver/dgConservationLawPerfectGas.h
@@ -1,6 +1,11 @@
#ifndef _DG_CONSERVATION_LAW_PERFECT_GAS_H_
#define _DG_CONSERVATION_LAW_PERFECT_GAS_H_
#include "dgConservationLaw.h"
+// Compressible Navier-Stokes equations
+// d(rho)/dt + d(rhou)/dx + d(rhov)/dy = 0
+// d(rhou)/dt + d(rhou^2+p)/dx + d(rhouv)/dy + d(t_xx)/dx + d(t_xy)/dy= 0
+// d(rhov)/dt + d(rhouv)/dx + d(rhov^2+p)/dy + d(t_xy)/dx + d(t_yy)/dy= 0
+// d(rhoE)/dt + d(rhouH)/dx + d(rhovH)/dy + d(ut_xx+vt_xy-q_x)/dx + d(ut_xy+vt_yy-q_y)/dy= 0
class dgPerfectGasLaw2d : public dgConservationLaw {
class advection;
class diffusion;
@@ -9,6 +14,7 @@ class dgPerfectGasLaw2d : public dgConservationLaw {
class source;
class maxConvectiveSpeed;
class maxDiffusivity;
+ class clipToPhysics;
// the name of the functions for
// viscosity (_muFunctionName)
// thermal conductivity (_kappaFunctionName)
@@ -22,6 +28,7 @@ class dgPerfectGasLaw2d : public dgConservationLaw {
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
+ dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const;
dgPerfectGasLaw2d();
dgBoundaryCondition *newWallBoundary() ;
dgBoundaryCondition *newSlipWallBoundary() ;
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index e7c27a94b790c8e6032bb37aa822477a6fd7335b..30d34e3a2ee835b27194c438534ff6042e4730b2 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -1,9 +1,7 @@
#include "dgConservationLaw.h"
#include "dgConservationLawWaveEquation.h"
#include "function.h"
-// dp/dt - rho*c^2 div(u,v) = 0
-// du/dt + 1/rho dp/dx = 0
-// dv/dt + 1/rho dp/dy = 0
+
static double c=1;
static double rho=1;
diff --git a/Solver/dgConservationLawWaveEquation.h b/Solver/dgConservationLawWaveEquation.h
index 0202905d09a34ffa14a68b8d23b316be2cc33e04..8086b61ae7f73fc21899ea1d66c5cbfedd0cb61e 100644
--- a/Solver/dgConservationLawWaveEquation.h
+++ b/Solver/dgConservationLawWaveEquation.h
@@ -1,6 +1,10 @@
#ifndef _DG_CONSERVATION_LAW_WAVE_EQUATION_H_
#define _DG_CONSERVATION_LAW_WAVE_EQUATION_H_
#include "dgConservationLaw.h"
+// Linear acoustics
+// dp/dt - rho*c^2 div(u,v) = 0
+// du/dt + 1/rho dp/dx = 0
+// dv/dt + 1/rho dp/dy = 0
class methodBinding;
class constructorBinding;
class dgConservationLawWaveEquation : public dgConservationLaw {
diff --git a/Solver/dgLimiter.h b/Solver/dgLimiter.h
index 7d8cab78fa279020485a3bd572bc475f0b8fb0e3..860dd4282252b4cd9a8a82cd9a40aaf7a685e0c1 100644
--- a/Solver/dgLimiter.h
+++ b/Solver/dgLimiter.h
@@ -6,17 +6,20 @@
struct dgDofContainer;
class dgGroupOfElements;
class dgGroupOfFaces;
+class dgConservationLaw;
class dgLimiter{
+protected:
+ dgConservationLaw *_claw;
public:
- dgLimiter () {}
+ dgLimiter (dgConservationLaw *claw) : _claw(claw) {}
virtual bool apply ( dgDofContainer &sol, std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroup ) = 0;
};
class dgSlopeLimiter : public dgLimiter{
public :
- dgSlopeLimiter () : dgLimiter () {}
+ dgSlopeLimiter (dgConservationLaw *claw) : dgLimiter (claw) {}
virtual bool apply ( dgDofContainer &solution,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroup);
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index ba4f410d282cafc80c2a8e35ef1298612a0433d8..d8a205f18a2651b3b0df5a0c3f8a2cc82cd17692 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/Solver/dgSlopeLimiter.cpp
@@ -1,7 +1,7 @@
#include "dgLimiter.h"
#include "dgGroupOfElements.h"
#include "dgSystemOfEquations.h"
-
+#include "function.h"
//----------------------------------------------------------------------------------
bool dgSlopeLimiter::apply ( dgDofContainer &solution,
@@ -13,13 +13,11 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
//TODO: make this more general
dgGroupOfFaces* group = fGroups[0];
- fullMatrix<double> &SolLeft = *(solution._dataProxys[0]);
- fullMatrix<double> &SolRight = *(solution._dataProxys[0]);
- int nbFields = solution.getNbFields();
+ fullMatrix<double> &solleft = solution.getGroupProxy(0); //*(solution._dataProxys[0]);
+ fullMatrix<double> &solright = solution.getGroupProxy(0); //*(solution._dataProxys[0]);
+ int nbFields =_claw->nbFields();
int totNbElems = solution.getNbElements();
-
- // --- CLIPPING: check unphysical values
// first compute max and min of all fields for all stencils
//----------------------------------------------------------
@@ -29,14 +27,14 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
MAX.setAll (-1.e22);
int iElementL, iElementR, fSize;
- for(int iFace=0 ; iFace<group->getNbElements();iFace++)
- {
+ for(int iFace=0 ; iFace<group->getNbElements();iFace++) {
+
iElementL = group->getElementLeftId(iFace);
iElementR = group->getElementRightId(iFace);
fullMatrix<double> TempL, TempR;
- TempL.setAsProxy(SolLeft, nbFields*iElementL, nbFields );
- TempR.setAsProxy(SolRight, nbFields*iElementR, nbFields );
+ TempL.setAsProxy(solleft, nbFields*iElementL, nbFields );
+ TempR.setAsProxy(solright, nbFields*iElementR, nbFields );
fSize = TempL.size1();
for (int k=0; k< nbFields; ++k){
@@ -55,16 +53,13 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
}
}
- // printf("fSize = %d\n",fSize);
-
- //----------------------------------------------------------
+ //----------------------------------------------------------
// then limit the solution
//----------------------------------------------------------
- for (int iElement=0 ; iElement<totNbElems; ++iElement)
- {
+ for (int iElement=0 ; iElement<totNbElems; ++iElement) {
fullMatrix<double> Temp;
- Temp.setAsProxy(SolRight, nbFields*iElement, nbFields );
+ Temp.setAsProxy(solleft, nbFields*iElement, nbFields );
for (int k=0; k<nbFields; ++k)
{
double AVG = 0.;
@@ -90,9 +85,7 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
if (AVG < neighMin) slopeLimiterValue = 0;
if (AVG > neighMax) slopeLimiterValue = 0;
- // if (slopeLimiterValue != 1.0){
- // printf("LIMTING %g\n",slopeLimiterValue);
- // }
+ // if (slopeLimiterValue != 1.0) printf("LIMTING %g\n",slopeLimiterValue);
// slopeLimiterValue = 0.0;
for (int i=0; i<fSize; ++i) Temp(i,k) = AVG + Temp(i,k)*slopeLimiterValue;
@@ -100,31 +93,50 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
}
}
- #if 0
- double rhomin = 1.e-3;
- double presmin= 1.e-3;
- for (int iElement=0;iElement<totNbElems;iElement++){
- fullMatrix<double> solElem;
- solElem.setAsProxy(SolRight, nbFields*iElement, nbFields );
- for (int k=0;k< solElem.size1() ;k++){
- if (solElem(k,0) < rhomin) {
- solElem(k,0) = rhomin;
- }
- double rhoV2 = 0;
- for (int j=0;j<2;j++) {
- double rhov = solElem(k,j+1);
- rhoV2 += rhov*rhov;
- }
- rhoV2 /= solElem(k,0);
- const double p = (1.4-1)*solElem(k,3) - 0.5*(1.4-1)* rhoV2;
- if (p < presmin) {
- solElem(k,3) = 0.5 *rhoV2 + presmin / (1.4-1);
- // printf("negative pressure %g cliiped !!\n",p);
- }
- }
- }
- #endif
+ // --- 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);
+
+ //#if 0
+// double rhomin = 1.e-3;
+// double presmin= 1.e-3;
+// for (int iElement=0;iElement<totNbElems;iElement++){
+// fullMatrix<double> solElem;
+// solElem.setAsProxy(solleft, nbFields*iElement, nbFields );
+// for (int k=0;k< solElem.size1() ;k++){
+// if (solElem(k,0) < rhomin) {
+// solElem(k,0) = rhomin;
+// }
+// double rhoV2 = 0;
+// for (int j=0;j<2;j++) {
+// double rhov = solElem(k,j+1);
+// rhoV2 += rhov*rhov;
+// }
+// rhoV2 /= solElem(k,0);
+// const double p = (1.4-1)*solElem(k,3) - 0.5*(1.4-1)* rhoV2;
+// if (p < presmin) {
+// solElem(k,3) = 0.5 *rhoV2 + presmin / (1.4-1);
+// }
+// }
+// }
+ //#endif
return true;
}
+
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 2d30b2923b6ae9dccf078263d6e787e47a048007..a4d074086caf3e53c7e922f415c6e500e3b018ee 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -103,8 +103,7 @@ double dgSystemOfEquations::RK44(double dt){
return _solution->_data->norm();
}
-double dgSystemOfEquations::computeInvSpectralRadius(){
-
+double dgSystemOfEquations::computeInvSpectralRadius(){
double sr = 1.e22;
for (int i=0;i<_elementGroups.size();i++){
std::vector<double> DTS;
@@ -115,7 +114,7 @@ double dgSystemOfEquations::computeInvSpectralRadius(){
}
double dgSystemOfEquations::RK44_limiter(double dt){
- dgLimiter *sl = new dgSlopeLimiter();
+ dgLimiter *sl = new dgSlopeLimiter(_claw);
_algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, sl);
delete sl;
return _solution->_data->norm();
@@ -135,7 +134,7 @@ void dgSystemOfEquations::exportSolution(std::string outputFile){
}
void dgSystemOfEquations::limitSolution(){
- dgLimiter *sl = new dgSlopeLimiter();
+ dgLimiter *sl = new dgSlopeLimiter(_claw);
sl->apply(*_solution,_elementGroups,_faceGroups);
delete sl;
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 22b11b795c8b68790a440308f8b9a70f9dc6bf78..f92ae8b44e75d4224788e468387373808bab85b3 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -23,7 +23,6 @@ public:
dgDofContainer (std::vector<dgGroupOfElements*> &groups, const dgConservationLaw &claw);
~dgDofContainer ();
int getNbElements() {return totalNbElements;}
- int getNbFields() {return nbFields;}
};
class binding;