diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 7b04811ea66a5a59bf656b9e6c7d0ee915f63c1c..4426632c904453b067cfd7da944252ae9f8cdefa 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -97,7 +97,6 @@
#define MSH_QUA_16I 40
#define MSH_QUA_20 41
-
#define MSH_NUM_TYPE 35
// Geometric entities
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index a4f906ffa2d6625a70db2832551844d2ca96053e..3f5d79bcc3d6c5c06eb0b0b68911286003a4cf29 100644
--- a/Geo/MLine.cpp
+++ b/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"
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 04be90f28b0a152f7bdf67e49e63c8d1dfc4d7d4..ab7171141eb1960a537cfd53b2d6a06dbdda0c8e 100644
--- a/Geo/MTriangle.cpp
+++ b/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"
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 22216602059f4cfe8515828c6952e18a5de5ba23..c7ec58263bfb12fc99bb4491269f4c078ffd12ff 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/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)
diff --git a/Solver/TESTCASES/Advection1D.lua b/Solver/TESTCASES/Advection1D.lua
index 8df7e25056d5180acb96a72e5fdca6112291f35b..0caaa663be1564ab5767a42199d8968ef374bb53 100644
--- a/Solver/TESTCASES/Advection1D.lua
+++ b/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))
diff --git a/Solver/TESTCASES/Advection3D.lua b/Solver/TESTCASES/Advection3D.lua
index e7513854a4fb695ccb379cd5b9abccfa4c4dfd5b..c00bb6d9d609db046abb0469be95c50786e08df1 100644
--- a/Solver/TESTCASES/Advection3D.lua
+++ b/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))
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index ff7dcf9c5e046aa284c8e50bd1bca38fe22550aa..a86a0d4305f7da8090feaf1670215968bcee8253 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/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)
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index d66c9dfe824e4f1d6f52d249eaf2b4f689444b96..937866f44fb3fbad503bf9444137cc3ea279da0f 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/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);
}
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 6a4460a637fd105f470da4842085a5f1823e7add..2bc50dc249558680184dc570257d51bdd5b5cd22 100644
--- a/Solver/dgConservationLaw.h
+++ b/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);
};
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index c57c47a8ccb326325ef41b300b0aded3c9e8a02b..f93e120e72eb265980e920562eb2acbb7e88e7f0 100644
--- a/Solver/dgGroupOfElements.cpp
+++ 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);
}
}
}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 7d97db7081fbb150fe7103d645a6aee34a48382e..c1f7de0c26fa2e8cc963dc7c24feadbe8ccfb2fe 100644
--- a/Solver/dgGroupOfElements.h
+++ b/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]];}
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index b0b11effaf54420a6b140ffedc875c0de31f7b7c..33cd53fb712e46d994f3da868a3ccd9c3059acc9 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/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;