Skip to content
Snippets Groups Projects
Commit 2c4561b5 authored by Stefan Girke's avatar Stefan Girke
Browse files

adjust flux discrete model, still not compiling (advdiff-test)

parent 475ecc9a
No related branches found
No related tags found
No related merge requests found
...@@ -167,7 +167,7 @@ namespace Fem ...@@ -167,7 +167,7 @@ namespace Fem
void analyticalFlux(const LocalEvaluation& local, void analyticalFlux(const LocalEvaluation& local,
JacobianRangeType& f) JacobianRangeType& f)
{ {
model_.jacobian( local.entity(), local.time(), local.position(), local.values()[ uVar ], f); model_.jacobian( local, local.values()[ uVar ], f);
} }
/** /**
...@@ -178,41 +178,27 @@ namespace Fem ...@@ -178,41 +178,27 @@ namespace Fem
RangeType& gLeft, RangeType& gLeft,
JacobianRangeType& gDiffLeft ) const JacobianRangeType& gDiffLeft ) const
{ {
return boundaryFluxImpl( left.intersection(), left.time(), return boundaryFluxImpl( left, left.values()[ uVar ], gLeft, gDiffLeft );
left.quadrature(), left.index(),
left.values()[ uVar ], gLeft, gDiffLeft );
} }
protected: protected:
template <class QuadratureImp, template <class LocalEvaluation, class UType >
class UType> double boundaryFluxImpl( const LocalEvaluation& left,
double boundaryFluxImpl(const Intersection& it, const UType& uLeft,
double time, RangeType& gLeft,
const QuadratureImp& faceQuadInner, JacobianRangeType& gDiffLeft ) const
const int quadPoint,
const UType& uLeft,
RangeType& gLeft,
JacobianRangeType& gDiffLeft ) const
{ {
const FaceDomainType& x = faceQuadInner.localPosition( quadPoint ); const FaceDomainType& x = left.localPosition();
const DomainType normal = left.intersection().integrationOuterNormal( x );
UType uRight; UType uRight;
if( model_.hasBoundaryValue(it,time,x) ) if( model_.hasBoundaryValue( left ) )
{ model_.boundaryValue( left, uLeft, uRight);
model_.boundaryValue(it, time, x, uLeft, uRight);
}
else else
{
uRight = uLeft; uRight = uLeft;
}
return gradientFlux_.gradientBoundaryFlux(it, inside(), return gradientFlux_.gradientBoundaryFlux( left, uLeft, uRight, gLeft, gDiffLeft );
time, faceQuadInner, quadPoint,
uLeft,
uRight,
gLeft,
gDiffLeft );
} }
private: private:
...@@ -344,7 +330,7 @@ namespace Fem ...@@ -344,7 +330,7 @@ namespace Fem
if (diffusion) if (diffusion)
{ {
const double dtStiff = const double dtStiff =
model_.stiffSource( local.entity(), local.time(), local.position(), local.values()[uVar], uJac, s ); model_.stiffSource( local, local.values()[uVar], uJac, s );
dtEst = ( dtStiff > 0 ) ? dtStiff : dtEst; dtEst = ( dtStiff > 0 ) ? dtStiff : dtEst;
} }
...@@ -352,7 +338,7 @@ namespace Fem ...@@ -352,7 +338,7 @@ namespace Fem
{ {
RangeType sNonStiff(0); RangeType sNonStiff(0);
const double dtNon = const double dtNon =
model_.nonStiffSource( local.entity(), local.time(), local.position(), local.values()[uVar], uJac, sNonStiff ); model_.nonStiffSource( local, local.values()[uVar], uJac, sNonStiff );
s += sNonStiff; s += sNonStiff;
...@@ -410,8 +396,7 @@ namespace Fem ...@@ -410,8 +396,7 @@ namespace Fem
{ {
RangeType dLeft, dRight; RangeType dLeft, dRight;
diffTimeStep = diffTimeStep =
diffFlux_.numericalFlux(left.intersection(), *this, diffFlux_.numericalFlux(left,
left.time(), left.quadrature(), right.quadrature(), left.index(),
left.values()[ uVar ], right.values()[ uVar ], left.values()[ uVar ], right.values()[ uVar ],
left.values() [ sigmaVar ], right.values()[ sigmaVar ], left.values() [ sigmaVar ], right.values()[ sigmaVar ],
dLeft, dRight, dLeft, dRight,
...@@ -445,20 +430,18 @@ namespace Fem ...@@ -445,20 +430,18 @@ namespace Fem
double diffTimeStep = 0.0; double diffTimeStep = 0.0;
bool hasBoundaryValue = bool hasBoundaryValue = model_.hasBoundaryValue( left );
model_.hasBoundaryValue( left );
if( diffusion && hasBoundaryValue ) if( diffusion && hasBoundaryValue )
{ {
// diffusion boundary flux for Dirichlet boundaries // diffusion boundary flux for Dirichlet boundaries
RangeType dLeft ( 0 ); RangeType dLeft ( 0 );
diffTimeStep = diffFlux_.boundaryFlux( left.intersection(), diffTimeStep = diffFlux_.boundaryFlux( left,
*this, left.values()[ uVar ],
left.time(), left.quadrature(), left.index(), uBnd_, // is set during call of BaseType::boundaryFlux
left.values()[ uVar ], uBnd_, // is set during call of BaseType::boundaryFlux left.values()[ sigmaVar ],
left.values()[ sigmaVar ], dLeft,
dLeft, gDiffLeft);
gDiffLeft);
gLeft += dLeft; gLeft += dLeft;
} }
else if ( diffusion ) else if ( diffusion )
...@@ -468,8 +451,7 @@ namespace Fem ...@@ -468,8 +451,7 @@ namespace Fem
typedef typename DiffusionFluxType::GradientRangeType GradientRangeType; typedef typename DiffusionFluxType::GradientRangeType GradientRangeType;
Dune::Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType > uJac( left.values()[ sigmaVar ] ); Dune::Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType > uJac( left.values()[ sigmaVar ] );
model_.diffusionBoundaryFlux( left.intersection(), left.time(), left.localPosition(), model_.diffusionBoundaryFlux( left, left.values()[uVar], uJac, diffBndFlux );
left.values()[uVar], uJac, diffBndFlux );
gLeft += diffBndFlux; gLeft += diffBndFlux;
} }
else else
...@@ -498,7 +480,7 @@ namespace Fem ...@@ -498,7 +480,7 @@ namespace Fem
typedef typename DiffusionFluxType::GradientRangeType GradientRangeType; typedef typename DiffusionFluxType::GradientRangeType GradientRangeType;
Dune::Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType > uJac( local.values()[ sigmaVar ] ); Dune::Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType > uJac( local.values()[ sigmaVar ] );
model_.diffusion( local.entity(), local.time(), local.position(), local.values()[ uVar ], uJac, diffmatrix); model_.diffusion( local, local.values()[ uVar ], uJac, diffmatrix);
// ldg case // ldg case
f += diffmatrix; f += diffmatrix;
} }
......
...@@ -144,6 +144,7 @@ namespace Fem ...@@ -144,6 +144,7 @@ namespace Fem
if( thetaLeft > 0 ) if( thetaLeft > 0 )
{ {
//TODO use diffusionTimeStep
diffTimeStep = diffTimeStep =
/* central differences (might be suboptimal) */ /* central differences (might be suboptimal) */
model_.diffusion(left, uLeft, diffmatrix ); model_.diffusion(left, uLeft, diffmatrix );
...@@ -204,6 +205,7 @@ namespace Fem ...@@ -204,6 +205,7 @@ namespace Fem
if( thetaRight > 0 ) if( thetaRight > 0 )
uVal.axpy( thetaRight, uBnd ); uVal.axpy( thetaRight, uBnd );
//TODO use diffusionTimeStep
const double diffTimeStep = const double diffTimeStep =
model_.diffusion(left, uVal, diffmatrix ); model_.diffusion(left, uVal, diffmatrix );
...@@ -234,8 +236,6 @@ namespace Fem ...@@ -234,8 +236,6 @@ namespace Fem
const LocalEvaluation& right, const LocalEvaluation& right,
const RangeType& uLeft, const RangeType& uLeft,
const RangeType& uRight, const RangeType& uRight,
const JacobianRangeType& jacLeft,
const JacobianRangeType& jacRight,
const GradientRangeType& sigmaLeft, const GradientRangeType& sigmaLeft,
const GradientRangeType& sigmaRight, const GradientRangeType& sigmaRight,
RangeType& gLeft, RangeType& gLeft,
...@@ -309,8 +309,10 @@ namespace Fem ...@@ -309,8 +309,10 @@ namespace Fem
/* Diffusion (Pass 2) */ /* Diffusion (Pass 2) */
/****************************/ /****************************/
JacobianRangeType diffmatrix; JacobianRangeType diffmatrix;
//TODO use diffusionTimeStep
double diffTimeStep = double diffTimeStep =
model_.diffusion(left, uLeft, sigmaLeft, diffmatrix); model_.diffusion(left, uLeft, sigmaLeft, diffmatrix);
diffmatrix.mv(normal, gLeft); diffmatrix.mv(normal, gLeft);
// add penalty term // add penalty term
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment