From 4469c33414f7e81df242142905bb6ad9ab4b7700 Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Fri, 11 Dec 2009 10:56:47 +0000
Subject: [PATCH] dg: IP
---
Numeric/fullMatrix.cpp | 1 +
Numeric/fullMatrix.h | 1 +
.../{Advection.lua => AdvectionDiffusion.lua} | 6 +-
Solver/TESTCASES/Diffusion.lua | 43 ++++++++++
Solver/dgAlgorithm.cpp | 82 ++++++++++++-------
Solver/dgAlgorithm.h | 2 +-
Solver/dgConservationLaw.h | 4 +-
Solver/dgConservationLawAdvection.cpp | 47 +++++++++--
Solver/dgGroupOfElements.cpp | 8 +-
Solver/dgGroupOfElements.h | 12 +++
Solver/dgSystemOfEquations.cpp | 3 +-
11 files changed, 165 insertions(+), 44 deletions(-)
rename Solver/TESTCASES/{Advection.lua => AdvectionDiffusion.lua} (91%)
create mode 100644 Solver/TESTCASES/Diffusion.lua
diff --git a/Numeric/fullMatrix.cpp b/Numeric/fullMatrix.cpp
index 6b601c9d8a..b7d4c32011 100644
--- a/Numeric/fullMatrix.cpp
+++ b/Numeric/fullMatrix.cpp
@@ -289,6 +289,7 @@ bool fullMatrix<double>::svd(fullMatrix<double> &V, fullVector<double> &S)
return false;
}
+
#if defined(HAVE_LUA)
template<>
int fullMatrix<double>::gemm(lua_State *L)
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index 94fd3c8f07..d6b65e37db 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -134,6 +134,7 @@ class fullMatrix
lua_pushnumber(L, (*this)(r,c));
return 1;
}
+
int set(lua_State *L){
int r = luaL_checkint(L, 1);
int c = luaL_checkint(L, 2);
diff --git a/Solver/TESTCASES/Advection.lua b/Solver/TESTCASES/AdvectionDiffusion.lua
similarity index 91%
rename from Solver/TESTCASES/Advection.lua
rename to Solver/TESTCASES/AdvectionDiffusion.lua
index 5d4f9f0204..68cb64479b 100644
--- a/Solver/TESTCASES/Advection.lua
+++ b/Solver/TESTCASES/AdvectionDiffusion.lua
@@ -11,7 +11,11 @@ v=fullMatrix(3,1);
v:set(0,0,0.15)
v:set(1,0,0.05)
v:set(2,0,0)
-dg:setConservationLaw('AdvectionDiffusion',createFunction.constant(v))
+-- diffusivity
+nu=fullMatrix(1,1);
+nu:set(0,0,0.001)
+
+dg:setConservationLaw('AdvectionDiffusion',createFunction.constant(v),createFunction.constant(nu))
-- boundary condition
outside=fullMatrix(1,1)
diff --git a/Solver/TESTCASES/Diffusion.lua b/Solver/TESTCASES/Diffusion.lua
new file mode 100644
index 0000000000..99f20e84de
--- /dev/null
+++ b/Solver/TESTCASES/Diffusion.lua
@@ -0,0 +1,43 @@
+model = GModel ()
+model:load ('square.geo')
+model:load ('square.msh')
+dg = dgSystemOfEquations (model)
+dg:setOrder(5)
+
+
+-- conservation law
+-- advection speed
+nu=fullMatrix(1,1);
+nu:set(0,0,0.01)
+dg:setConservationLaw('AdvectionDiffusion','',createFunction.constant(nu))
+
+-- boundary condition
+outside=fullMatrix(1,1)
+outside:set(0,0,0.)
+dg:addBoundaryCondition('Border','0Flux')
+
+dg:setup()
+
+-- initial condition
+function initial_condition( _x , _f )
+ xyz = fullMatrix(_x)
+ f = fullMatrix(_f)
+ for i=0,xyz:size1()-1 do
+ x = xyz:get(i,0)
+ y = xyz:get(i,1)
+ z = xyz:get(i,2)
+ f:set (i, 0, math.exp(-100*((x-0.2)^2 +(y-0.3)^2)))
+ end
+end
+dg:L2Projection(createFunction.lua(1,'initial_condition','XYZ'))
+
+dg:exportSolution('output/Diffusion_00000')
+
+-- main loop
+for i=1,10000 do
+ norm = dg:RK44(0.001)
+ if (i % 50 == 0) then
+ print('iter',i,norm)
+ dg:exportSolution(string.format("output/Diffusion-%05d", i))
+ end
+end
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 718d2cd8f4..3c32bb7d70 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -16,7 +16,7 @@
void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group,
- const fullMatrix<double> &solution,
+ fullMatrix<double> &solution,
fullMatrix<double> &residual // the residual
)
{
@@ -45,27 +45,26 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
// provided dataCache
cacheMap.provideData("UVW").set(group.getIntegrationPointsMatrix());
dataCacheDouble &solutionQPe = cacheMap.provideData("Solution");
- dataCacheDouble &gradSolutionQPe = cacheMap.provideData("gradSolution");
+ dataCacheDouble &gradientSolutionQPe = cacheMap.provideData("SolutionGradient");
+ gradientSolutionQPe.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
// needed dataCache
dataCacheDouble *sourceTerm = claw.newSourceTerm(cacheMap);
// fluxes in XYZ coordinates;
dataCacheDouble *convectiveFlux = claw.newConvectiveFlux(cacheMap);
dataCacheDouble *diffusiveFlux = claw.newDiffusiveFlux(cacheMap);
- // compute the gradient of solution if necessary
- fullMatrix<double> *gradientSolutionQP= 0;
- if (gradSolutionQPe.somethingDependOnMe()) {
- gradientSolutionQP = new fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields * 3);
- //group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP);
- }
fullMatrix<double> fuvwe;
fullMatrix<double> source;
+ fullMatrix<double> dPsiDx,dofs;
// ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
solutionQPe.setAsProxy(solutionQP, iElement*nbFields, nbFields );
- if (gradSolutionQPe.somethingDependOnMe())
- gradSolutionQPe.setAsProxy(*gradientSolutionQP, 3*iElement*nbFields,3*nbFields );
+ if(gradientSolutionQPe.somethingDependOnMe()){
+ dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
+ dofs.setAsProxy(solution, nbFields*iElement, nbFields);
+ dPsiDx.mult(dofs, gradientSolutionQPe.set());
+ }
cacheElement.set(group.getElement(iElement));
if(convectiveFlux || diffusiveFlux) {
// ----- 2.3.3 --- compute fluxes in UVW coordinates
@@ -84,8 +83,9 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
for (int k=0;k<nbFields;k++){
if(convectiveFlux)
fuvwe(iPt,k) += ( (*convectiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
- if(diffusiveFlux)
+ if(diffusiveFlux){
fuvwe(iPt,k) += ( (*diffusiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
+ }
}
}
}
@@ -100,8 +100,6 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
}
}
}
- if(gradientSolutionQP)
- delete gradientSolutionQP;
}
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
if(convectiveFlux || diffusiveFlux){
@@ -158,8 +156,8 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
dataCacheDouble &uvwRight = cacheMapRight.provideData("UVW");
dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution");
- dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("GradientSolution");
- dataCacheDouble &gradientSolutionRight = cacheMapRight.provideData("GradientSolution");
+ dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient");
+ dataCacheDouble &gradientSolutionRight = cacheMapRight.provideData("SolutionGradient");
//resize gradientSolutionLeft and Right
gradientSolutionLeft.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
gradientSolutionRight.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
@@ -167,9 +165,14 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
// A2 ) ask the law to provide the fluxes (possibly NULL)
dataCacheDouble *riemannSolver = claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
dataCacheDouble *diffusiveFluxLeft = claw.newDiffusiveFlux(cacheMapLeft);
+ dataCacheDouble *maximumDiffusivityLeft = claw.newMaximumDiffusivity(cacheMapLeft);
+ dataCacheDouble *maximumDiffusivityRight = claw.newMaximumDiffusivity(cacheMapRight);
dataCacheDouble *diffusiveFluxRight = claw.newDiffusiveFlux(cacheMapRight);
fullMatrix<double> dPsiLeftDx,dPsiRightDx,dofsLeft,dofsRight,normalFluxQP;
+ int p = group.getGroupLeft().getOrder();
+ int dim = group.getGroupLeft().getElement(0)->getDim();
+
for (int iFace=0 ; iFace < nbFaces ; ++iFace) {
// B1 ) adjust the proxies for this element
// NB : the B1 section will almost completely disapear
@@ -183,7 +186,7 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
normals.setAsProxy(group.getNormals(), iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
// proxies needed to compute the gradient of the solution and the IP term
dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*nbNodesLeft,nbNodesLeft);
- dPsiRightDx.setAsProxy(DPsiLeftDx,iFace*nbNodesRight,nbNodesRight);
+ dPsiRightDx.setAsProxy(DPsiRightDx,iFace*nbNodesRight,nbNodesRight);
dofsLeft.setAsProxy(solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
dofsRight.setAsProxy(solutionRight, nbFields*group.getElementRightId(iFace), nbFields);
uvwLeft.setAsProxy(group.getIntegrationOnElementLeft(iFace));
@@ -198,27 +201,44 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
}
// B3 ) compute fluxes
-
- //JF : here you can test (diffusiveFluxLeft!=NULL) to know if the law is diffusive
- //you can access the values by (*diffusiveFluxLeft)(iQP*3+iXYZ, iField)
- //use gradientSolutionLef(IQP*3+IXYZ, iField) and dPsiLeftDx(iQP*3+iXYZ, iPsi)
-
- for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
- const double detJ = group.getDetJ (iFace, iPt);
- for (int k=0;k<nbFields*2;k++)
- normalFluxQP(iPt,k) = (*riemannSolver)(iPt,k)*detJ;
+ if (diffusiveFluxLeft) {
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iFace, iPt);
+ //just for the lisibility :
+ const fullMatrix<double> &dfl = (*diffusiveFluxLeft)();
+ const fullMatrix<double> &dfr = (*diffusiveFluxRight)();
+ for (int k=0;k<nbFields;k++) {
+ double meanNormalFlux =
+ ((dfl(iPt,k*3 )+dfr(iPt,k*3 )) * normals(0,iPt)
+ +(dfl(iPt,k*3+1)+dfr(iPt,k*3+1)) * normals(1,iPt)
+ +(dfl(iPt,k*3+2)+dfr(iPt,k*3+2)) * normals(2,iPt))/2;
+ double minl = std::min(group.getElementVolumeLeft(iFace),
+ group.getElementVolumeRight(iFace)
+ )/group.getInterfaceSurface(iFace);
+ double nu = std::max((*maximumDiffusivityRight)(iPt,0),(*maximumDiffusivityLeft)(iPt,0));
+ double mu = (p+1)*(p+dim)/dim*nu/minl;
+ double solutionJumpPenalty = (solutionQPLeft(iPt,k)-solutionQPRight(iPt,k))/2*mu;
+ normalFluxQP(iPt,k) -= (meanNormalFlux+solutionJumpPenalty)*detJ;
+ normalFluxQP(iPt,k+nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ;
+ }
+ }
+ }
+ if (riemannSolver) {
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iFace, iPt);
+ for (int k=0;k<nbFields*2;k++)
+ normalFluxQP(iPt,k) += (*riemannSolver)(iPt,k)*detJ;
+ }
}
}
// C ) redistribute the flux to the residual (at Faces nodes)
- residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
+ if(riemannSolver || diffusiveFluxLeft)
+ residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
// D ) delete the dataCacheDouble provided by the law
- delete riemannSolver;
- if(diffusiveFluxLeft) {
- delete diffusiveFluxLeft;
- delete diffusiveFluxRight;
- }
+ if(riemannSolver)
+ delete riemannSolver;
}
void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index fd9c321490..a8ad3eb857 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -15,7 +15,7 @@ class dgAlgorithm {
static void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group,
- const fullMatrix<double> &solution,
+ fullMatrix<double> &solution,
fullMatrix<double> &residual);
static void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 3e3c75a642..dbd7cdf68b 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -32,6 +32,8 @@ class dgConservationLaw {
virtual dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newDiffusiveFlux (dataCacheMap &cacheMap) const {return NULL;}
+ virtual dataCacheDouble *newDiffusiveInterfaceTerm (dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {return NULL;}
+ virtual dataCacheDouble *newMaximumDiffusivity (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newConvectiveFlux (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newMaxConvectiveSpeed (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newRiemannSolver (dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {return NULL;}
@@ -51,7 +53,7 @@ class dgConservationLaw {
};
-dgConservationLaw *dgNewConservationLawAdvection(const std::string vname);
+dgConservationLaw *dgNewConservationLawAdvection(const std::string vname,const std::string nuname);
dgConservationLaw *dgNewConservationLawShallowWater2d();
dgConservationLaw *dgNewConservationLawWaveEquation();
dgConservationLaw *dgNewPerfectGasLaw2d();
diff --git a/Solver/dgConservationLawAdvection.cpp b/Solver/dgConservationLawAdvection.cpp
index b29bdb89ce..ac44af08c3 100644
--- a/Solver/dgConservationLawAdvection.cpp
+++ b/Solver/dgConservationLawAdvection.cpp
@@ -6,7 +6,7 @@
#include "function.h"
class dgConservationLawAdvection : public dgConservationLaw {
- std::string _vFunctionName;
+ std::string _vFunctionName,_nuFunctionName;
class advection : public dataCacheDouble {
dataCacheDouble &sol, &v;
public:
@@ -46,23 +46,56 @@ class dgConservationLawAdvection : public dgConservationLaw {
}
}
};
+ class diffusion : public dataCacheDouble {
+ dataCacheDouble &solgrad, ν
+ public:
+ diffusion(std::string nuFunctionName, dataCacheMap &cacheMap):
+ solgrad(cacheMap.get("SolutionGradient",this)),
+ nu(cacheMap.get(nuFunctionName,this))
+ {};
+ void _eval () {
+ if(_value.size1() != solgrad().size1()/3)
+ _value=fullMatrix<double>(solgrad().size1()/3,3);
+ for(int i=0; i< solgrad().size1()/3; i++) {
+ _value(i,0) = -solgrad(i*3,0)*nu(i,0);
+ _value(i,1) = -solgrad(i*3+1,0)*nu(i,0);
+ _value(i,2) = -solgrad(i*3+1,0)*nu(i,0);
+ }
+ }
+ };
public:
dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const {
- return new advection(_vFunctionName,cacheMap);
+ if( !_vFunctionName.empty())
+ return new advection(_vFunctionName,cacheMap);
+ else
+ return NULL;
+ }
+ dataCacheDouble *newMaximumDiffusivity( dataCacheMap &cacheMap) const {
+ if( !_nuFunctionName.empty())
+ return &cacheMap.get(_nuFunctionName);
+ else
+ return NULL;
}
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
- return new riemann(_vFunctionName,cacheMapLeft, cacheMapRight);
+ if( !_vFunctionName.empty())
+ return new riemann(_vFunctionName,cacheMapLeft, cacheMapRight);
+ else
+ return NULL;
}
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const {
- return 0;
+ if( !_nuFunctionName.empty())
+ return new diffusion(_nuFunctionName,cacheMap);
+ else
+ return NULL;
}
- dgConservationLawAdvection(std::string vFunctionName)
+ dgConservationLawAdvection(std::string vFunctionName, std::string nuFunctionName)
{
_vFunctionName = vFunctionName;
+ _nuFunctionName = nuFunctionName;
_nbf = 1;
}
};
-dgConservationLaw *dgNewConservationLawAdvection(std::string vFunctionName) {
- return new dgConservationLawAdvection(vFunctionName);
+dgConservationLaw *dgNewConservationLawAdvection(std::string vFunctionName, std::string nuFunctionName) {
+ return new dgConservationLawAdvection(vFunctionName,nuFunctionName);
}
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index e3bd1a6fdc..62778f75d2 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -47,7 +47,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
_fs(*_elements[0]->getFunctionSpace(polyOrder)),
_integration(dgGetIntegrationRule (_elements[0], polyOrder))
{
-
+ _order=polyOrder;
_dimUVW = _dimXYZ = e[0]->getDim();
// this is the biggest piece of data ... the mappings
int nbPsi = _fs.coefficients.size1();
@@ -59,6 +59,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
_mapping = new fullMatrix<double> (e.size(), 10 * _integration->size1());
_imass = new fullMatrix<double> (nbPsi,nbPsi*e.size());
_dPsiDx = new fullMatrix<double> ( _integration->size1()*3,nbPsi*e.size());
+ _elementVolume = new fullMatrix<double> (e.size(),1);
double g[256][3],f[256];
for (int i=0;i<_elements.size();i++){
@@ -70,6 +71,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
double jac[3][3],ijac[3][3],detjac;
(*_mapping)(i,10*j + 9) =
e->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
+ (*_elementVolume)(i,0) += (*_mapping)(i,10*j+9)*(*_integration)(j,3);
const double weight = (*_integration)(j,3);
detjac=inv3x3(jac,ijac);
(*_mapping)(i,10*j + 0) = ijac[0][0];
@@ -123,6 +125,7 @@ dgGroupOfElements::~dgGroupOfElements(){
delete _mapping;
delete _collocation;
delete _imass;
+ delete _elementVolume;
}
@@ -204,6 +207,7 @@ void dgGroupOfFaces::init(int pOrder) {
_redistribution = new fullMatrix<double> (_fsFace->coefficients.size1(),_integration->size1());
_collocation = new fullMatrix<double> (_integration->size1(), _fsFace->coefficients.size1());
_detJac = new fullMatrix<double> (_integration->size1(), _faces.size());
+ _interfaceSurface = new fullMatrix<double>(_faces.size(),1);
for (size_t i=0; i<_closuresLeft.size(); i++)
_integrationPointsLeft.push_back(dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,_closuresLeft[i]));
for (size_t i=0; i<_closuresRight.size(); i++)
@@ -222,6 +226,7 @@ void dgGroupOfFaces::init(int pOrder) {
for (int j=0;j< _integration->size1() ; j++ ){
double jac[3][3];
(*_detJac)(j,i) = f->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
+ (*_interfaceSurface)(i,0) += (*_integration)(j,3)*(*_detJac)(j,i);
}
}
@@ -296,6 +301,7 @@ dgGroupOfFaces::~dgGroupOfFaces()
delete _normals;
delete _dPsiLeftDxOnQP;
delete _dPsiRightDxOnQP;
+ delete _interfaceSurface;
}
dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices):
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 793eb98bf0..ecf50829ee 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -63,6 +63,11 @@ class dgGroupOfElements {
// dimension of the parametric space and of the real space
// may be different if the domain is a surface in 3D (manifold)
int _dimUVW, _dimXYZ;
+ // polynomial order of the interpolation
+ int _order;
+ //
+ // volume/surface/length of the element (sum_qp w_i detJ_i)
+ fullMatrix<double> *_elementVolume;
// forbid the copy
// dgGroupOfElements (const dgGroupOfElements &e, int order) {}
// dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
@@ -79,11 +84,13 @@ public:
inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
inline const fullMatrix<double> & getFluxRedistributionMatrix (int i) const {return *_redistributionFluxes[i];}
inline const fullMatrix<double> & getSourceRedistributionMatrix () const {return *_redistributionSource;}
+ inline double getElementVolume (int iElement)const {return (*_elementVolume)(iElement,0);}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_mapping)(iElement, 10*iGaussPoint + 9);}
inline double getInvJ (int iElement, int iGaussPoint, int i, int j) const {return (*_mapping)(iElement, 10*iGaussPoint + i + 3*j);}
inline fullMatrix<double> & getDPsiDx() const { return *_dPsiDx;}
inline fullMatrix<double> &getInverseMassMatrix () const {return *_imass;}
inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
+ inline int getOrder() const {return _order;}
};
class dgGroupOfFaces {
@@ -127,6 +134,8 @@ class dgGroupOfFaces {
//fullMatrix<double> *_collocationLeft, *_collocationRight;
// redistribution matrices \psi_i (GP_j) * weight_j
fullMatrix<double> *_redistribution;
+ // surface/length/1 of the interface element (sum_qp w_i detJ_i)
+ fullMatrix<double> *_interfaceSurface;
//common part of the 3 constructors
void init(int pOrder);
public:
@@ -136,6 +145,8 @@ public:
inline int getElementRightId (int i) const {return _right[i];};
inline MElement* getElementLeft (int i) const {return _groupLeft.getElement(_left[i]);}
inline MElement* getElementRight (int i) const {return _groupRight.getElement(_right[i]);}
+ 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> &getClosureRight(int iFace) const{ return _closuresRight[_closuresIdRight[iFace]];}
@@ -160,6 +171,7 @@ public:
inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
inline const fullMatrix<double> & getRedistributionMatrix () const {return *_redistribution;}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_detJac)(iGaussPoint,iElement);}
+ inline double getInterfaceSurface (int iFace)const {return (*_interfaceSurface)(iFace,0);}
//keep this outside the Algorithm because this is the only place where data overlap
void mapToInterface(int nFields, const fullMatrix<double> &vLeft, const fullMatrix<double> &vRight, fullMatrix<double> &v);
void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight);
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 5609f35610..12add17a2b 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -69,8 +69,7 @@ int dgSystemOfEquations::setConservationLaw(lua_State *L){
else if (_cLawName == "PerfectGas2d")
_claw = dgNewPerfectGasLaw2d();
else if (_cLawName == "AdvectionDiffusion"){
- std::string advFunction = luaL_checkstring(L,2);
- _claw = dgNewConservationLawAdvection(advFunction);
+ _claw = dgNewConservationLawAdvection(luaL_checkstring(L,2),luaL_checkstring(L,3));
}
if (!_claw)throw;
return 0;
--
GitLab