Commit c2b6a404 authored by bertrand's avatar bertrand
Browse files

dg : function -> functor

git-svn-id: https://geuz.org/svn/gmsh/trunk@19207 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent 7e5f3846
......@@ -6,6 +6,7 @@
#include "dgFunctionIntegratorInterface.h"
#include "dgMeshJacobian.h"
#include "slimMovingBathWettingDrying.h"
#include "functorMember.h"
static double g = 9.80616;
......@@ -385,104 +386,93 @@ class dgConservationLawShallowWater2d::diffusivity : public function {
}
};
class dgConservationLawShallowWater2d::riemann:public function {
fullMatrix<double> solL, solR, normals, ip, _xyz;
fullMatrix<double> bath, movingBathFactorL, movingBathFactorR;
bool _linear, _isDiffusive;
double _R, _upwindFactor;
public:
riemann (const functor *bathymetry, const functor *ipTerm, bool linear, double R, const functor *movingBathFactorF, const double upwindFactorRiemann, const functor *xyz): function(3){
setArgument(solL,function::getSolution(), 0);
setArgument(solR,function::getSolution(), 1);
setArgument(bath,bathymetry, 0);
setArgument(normals,function::getNormals(), 0);
setArgument(movingBathFactorL, movingBathFactorF,0);
setArgument(movingBathFactorR, movingBathFactorF,1);
_isDiffusive = ipTerm;
if (_isDiffusive) {
setArgument(ip, ipTerm);
}
_linear = linear;
_R = R;
if(_R>0) {
setArgument(_xyz, xyz ? xyz : function::getCoordinates());
void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<double> &val) const {
const size_t nQP = cache.nEvaluationPoint();
val.resize(nQP,3);
const fullMatrix<double> &solL = cache.get(*function::getSolution(),0);
const fullMatrix<double> &solR = cache.get(*function::getSolution(),1);
const fullMatrix<double> &bath = cache.get(*_bathymetry,0);
const fullMatrix<double> &normals = cache.get(*function::getNormals(),0);
const fullMatrix<double> &movingBathFactorL = cache.get(*_movingBathFactor,0);
const fullMatrix<double> &movingBathFactorR = cache.get(*_movingBathFactor,1);
bool _isDiffusive = _ipTerm;
double _upwindFactor = _upwindFactorRiemann;
const fullMatrix<double> *ip = _isDiffusive ? &cache.get(*_ipTerm) : NULL;
const fullMatrix<double> *xyz = _R>0 ? &cache.get(*_xyz) : NULL;
int rev = 1;
double J = 1.0;
const dgGroupOfFaces &iGroup = *cache.interfaceGroup();
int connectionR = iGroup.nConnection() == 2 ? 1 : 0;
const dgGroupOfElements *groupR = &iGroup.elementGroup(connectionR);
const dgGroupOfElements *groupL = &iGroup.elementGroup(0);
for(size_t i=0; i< nQP; i++) {
rev = 1;
if(_R > 0)
{
double x = (*xyz)(i,0);
double y = (*xyz)(i,1);
J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
J = 1.0;
if(groupR->getPhysicalTag() != groupL->getPhysicalTag())
rev = -1;
}
_upwindFactor = upwindFactorRiemann;
};
void call(dataCacheMap *m, fullMatrix<double> &val) {
size_t nQP = solL.size1();
int rev = 1;
double J = 1.0;
const dgGroupOfFaces &iGroup = *m->interfaceGroup();
int connectionR = iGroup.nConnection() == 2 ? 1 : 0;
const dgGroupOfElements *groupR = &iGroup.elementGroup(connectionR);
const dgGroupOfElements *groupL = &iGroup.elementGroup(0);
for(size_t i=0; i< nQP; i++) {
rev = 1;
if(_R > 0)
{
double x = _xyz(i,0);
double y = _xyz(i,1);
J = 4.0*_R*_R/(4.0*_R*_R+x*x+y*y);
J = 1.0;
if(groupR->getPhysicalTag() != groupL->getPhysicalTag())
rev = -1;
}
double Fun, Fut, Feta;
double nx = normals(i,0);
double ny = normals(i,1);
double Fun, Fut, Feta;
double nx = normals(i,0);
double ny = normals(i,1);
/* if(rev<0)
{
double RR = sqrt(_xyz.get(i,0)*_xyz.get(i,0) + _xyz.get(i,1)*_xyz.get(i,1));
printf("nx = %.16g nx2 = %.16g, ny = %.16g ny2 = %.16g, 2R = %.16g , RR = %.16g\n ", nx,_xyz.get(i,0)/RR,ny,_xyz.get(i,1)/RR, 2*_R, RR);
{
double RR = sqrt(_xyz.get(i,0)*_xyz.get(i,0) + _xyz.get(i,1)*_xyz.get(i,1));
printf("nx = %.16g nx2 = %.16g, ny = %.16g ny2 = %.16g, 2R = %.16g , RR = %.16g\n ", nx,_xyz.get(i,0)/RR,ny,_xyz.get(i,1)/RR, 2*_R, RR);
}*/
double uL = nx*solL(i,1) + ny*solL(i,2), uR = nx*solR(i,1) + ny*solR(i,2);
double vL = -ny*solL(i,1) + nx*solL(i,2), vR = - ny*solR(i,1) + nx*solR(i,2);
double h = bath(i,0);
uR *= rev;
if (_linear) {
//linear equations
double etaR = solR(i,0);
double etaL = solL(i,0);
double sq_g_h = sqrt(g/h);
double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2;
double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2;
Fun = -g*etaRiemann*J;
Fut = 0;
Feta = -(h)*unRiemann*J;
}else{
double HR = solR(i,0) + h, HL = solL(i,0) + h;
if(HR<0 || HL<0) printf(" HR = %e HL =%e\n", HR,HL);
double u,v,H,Amb;
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor);
double eta = H-h;
Fun = -g*eta*J - u*u*J;
Fut = -u*v*J;
Feta = -H*u*J/Amb;
}
double uL = nx*solL(i,1) + ny*solL(i,2), uR = nx*solR(i,1) + ny*solR(i,2);
double vL = -ny*solL(i,1) + nx*solL(i,2), vR = - ny*solR(i,1) + nx*solR(i,2);
double h = bath(i,0);
uR *= rev;
if (_linear) {
//linear equations
double etaR = solR(i,0);
double etaL = solL(i,0);
double sq_g_h = sqrt(g/h);
double etaRiemann = (etaL+etaR + (uL-uR)/sq_g_h)/2;
double unRiemann = (uL+uR + (etaL-etaR)*sq_g_h)/2;
Fun = -g*etaRiemann*J;
Fut = 0;
Feta = -(h)*unRiemann*J;
}else{
double HR = solR(i,0) + h, HL = solL(i,0) + h;
if(HR<0 || HL<0) printf(" HR = %e HL =%e\n", HR,HL);
double u,v,H,Amb;
roeSolver(uL, uR, vL, vR, HL, HR, u, v, H, Amb, movingBathFactorL(i,0), movingBathFactorR(i,0), _upwindFactor);
double eta = H-h;
Fun = -g*eta*J - u*u*J;
Fut = -u*v*J;
Feta = -H*u*J/Amb;
}
val(i,0) = Feta;
val(i,1) = Fun * nx - Fut * ny;
val(i,2) = Fun * ny + Fut * nx;
/*val(i,3) = -val(i,0);
val(i,4) = - Fun * nx * rev + Fut * ny;
val(i,5) = - Fun * ny * rev - Fut * nx;*/
//printf("%f\n", val(i,0));
val(i,0) = Feta;
val(i,1) = Fun * nx - Fut * ny;
val(i,2) = Fun * ny + Fut * nx;
/*val(i,3) = -val(i,0);
val(i,4) = - Fun * nx * rev + Fut * ny;
val(i,5) = - Fun * ny * rev - Fut * nx;*/
//printf("%f\n", val(i,0));
if (_isDiffusive) {
val(i,1)+=ip(i,1);
val(i,2)+=ip(i,2);
/*val(i,4)+=ip(i,4)*rev;
val(i,5)+=ip(i,5)*rev;*/
}
if (_isDiffusive) {
val(i,1)+=(*ip)(i,1);
val(i,2)+=(*ip)(i,2);
/*val(i,4)+=ip(i,4)*rev;
val(i,5)+=ip(i,5)*rev;*/
}
}
};
}
void roeSolver(double uL, double uR, double vL, double vR, double HL, double HR,
double &uStar, double &vStar, double &HStar, double &AStar,
......@@ -584,11 +574,11 @@ void dgConservationLawShallowWater2d::setup() {
_sourceTerm = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _bathymetry, _bathymetryGradient, _nu, false || _linear, _R, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_gradPsiTerm = new gradPsiTerm(_bathymetry, _diffusiveFlux, false || _linear,_R, _movingBathFactor, _xyz);
_riemannTerm = new riemann(_bathymetry, _ipTerm, false || _linear, _R, _movingBathFactor, _upwindFactorRiemann, _xyz);
_riemannTerm = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeed = new maxConvectiveSpeed(_bathymetry, _linear, _R, false);
_sourceTermLin = new source(_linearDissipation, _quadraticDissipation, _coriolisFactor, _source, _bathymetry, _bathymetryGradient, _nu, true, _R, _movingBathFactor, _movingBathFactorGradient, _xyz, NULL, NULL);
_gradPsiTermLin = new gradPsiTerm(_bathymetry, _diffusiveFlux, true,_R, _movingBathFactor, _xyz);
_riemannTermLin = new riemann(_bathymetry, _ipTerm, true, _R, _movingBathFactor, _upwindFactorRiemann, _xyz);
_riemannTermLin = newFunctorMember(*this,&dgConservationLawShallowWater2d::riemann);
_maxSpeedFilt = new maxConvectiveSpeed(_bathymetry, _linear, _R, true);
setImexMode(_imexMode);
setLinearFilterMode(_linearFilter);
......
......@@ -21,7 +21,7 @@ class dgConservationLawShallowWater2d : public dgConservationLawFunction {
class source;
class diffusiveFlux;
class diffusivity;
class riemann;
void riemann(functorCache &cache,fullMatrix<double> &val) const;
class boundaryWall;
class boundaryForcedDischarge;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment