From 49f431281b7f885e80e6abfaa52916c72366db8e Mon Sep 17 00:00:00 2001
From: Robert K <robertk@posteo.org>
Date: Sat, 20 Feb 2016 16:42:11 +0100
Subject: [PATCH] [cleanup] Unified LDG and BR1 fluxes.
---
.../operator/fluxes/diffusion/averageflux.hh | 355 --------------
.../operator/fluxes/diffusion/fluxes.hh | 38 +-
.../operator/fluxes/diffusion/ldgflux.hh | 458 +++++++-----------
.../operator/fluxes/diffusion/parameters.hh | 8 +-
4 files changed, 218 insertions(+), 641 deletions(-)
delete mode 100644 dune/fem-dg/operator/fluxes/diffusion/averageflux.hh
diff --git a/dune/fem-dg/operator/fluxes/diffusion/averageflux.hh b/dune/fem-dg/operator/fluxes/diffusion/averageflux.hh
deleted file mode 100644
index c6425147..00000000
--- a/dune/fem-dg/operator/fluxes/diffusion/averageflux.hh
+++ /dev/null
@@ -1,355 +0,0 @@
-#ifndef DUNE_FEM_DG_LDGAVERAGEFLUX_HH
-#define DUNE_FEM_DG_LDGAVERAGEFLUX_HH
-
-// local includes
-#include "fluxbase.hh"
-
-namespace Dune
-{
-namespace Fem
-{
-
- /**
- * \brief A diffusion flux for the LDG scheme.
- *
- * \ingroup DiffusionFluxes
- */
- template <class DiscreteFunctionSpaceImp,
- class ModelImp,
- class FluxParameterImp >
- class LDGAverageDiffusionFlux :
- public LDGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp >
- {
- typedef LDGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp > BaseType;
- public:
- typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType;
-
- enum { dimDomain = DiscreteFunctionSpaceType :: dimDomain };
- enum { dimRange = DiscreteFunctionSpaceType :: dimRange };
-
- typedef typename DiscreteFunctionSpaceType :: DomainType DomainType;
- typedef typename DiscreteFunctionSpaceType :: RangeFieldType RangeFieldType;
- typedef typename DiscreteFunctionSpaceType :: DomainFieldType DomainFieldType;
- typedef typename DiscreteFunctionSpaceType :: RangeType RangeType;
- typedef typename DiscreteFunctionSpaceType :: JacobianRangeType JacobianRangeType;
-
- typedef FieldVector< DomainFieldType, dimDomain-1 > FaceDomainType;
-
- typedef typename DiscreteFunctionSpaceType :: GridPartType GridPartType;
- typedef typename GridPartType :: IntersectionIteratorType IntersectionIterator;
- typedef typename IntersectionIterator :: Intersection Intersection;
- typedef typename GridPartType :: GridType GridType;
- typedef typename DiscreteFunctionSpaceType :: EntityType EntityType;
- enum { dimGradRange = dimDomain * dimRange };
- enum { polOrd = DiscreteFunctionSpaceType :: polynomialOrder };
-
- typedef typename BaseType :: ParameterType ParameterType;
-
- // type of gradient space
- typedef typename DiscreteFunctionSpaceType ::
- template ToNewDimRange< dimGradRange > :: Type DiscreteGradientSpaceType;
-
- typedef typename DiscreteGradientSpaceType :: RangeType GradientRangeType;
- typedef typename DiscreteGradientSpaceType :: JacobianRangeType GradientJacobianType;
-
- // jacobians of the functions do not have to be evaluated for this flux
- enum { evaluateJacobian = false };
-
- protected:
- using BaseType::determineDirection;
- using BaseType::model_;
- using BaseType::cflDiffinv_;
- using BaseType::numericalFlux ;
- using BaseType::parameter ;
- using BaseType::nonconformingFactor_;
- using BaseType::dimensionFactor_;
-
- public:
- /**
- * \brief constructor
- */
- LDGAverageDiffusionFlux(GridPartType& gridPart,
- const ModelImp& mod,
- const ParameterType& param ) :
- BaseType( gridPart, mod, param ),
- penalty_( parameter().penalty() ),
- // Set CFL number for penalty term (compare diffusion in first pass)
- penaltyTerm_( std::abs( penalty_ ) > 0 )
- {
- if( Fem::Parameter::verbose () )
- {
- std::cout << "LDGAverageDiffusionFlux: penalty = " << penalty_ << std::endl;
- }
- }
-
- LDGAverageDiffusionFlux(const LDGAverageDiffusionFlux& other)
- : BaseType( other ),
- penalty_( other.penalty_ ),
- penaltyTerm_( other.penaltyTerm_ )
- {}
-
- //! returns true if lifting has to be calculated
- const bool hasLifting () const { return false; }
-
- protected:
- enum { realLDG = true };
- double theta( const DomainType& normal ) const
- {
- if( realLDG )
- {
- // LDG thete is 1 or 0
- if( determineDirection( normal ) )
- return 1.0;
- else
- return 0.0;
- }
- else
- // Average fluxes
- return 0.5;
- }
-
- public:
- /**
- * \brief flux function on interfaces between cells
- *
- * \param left local evaluation
- * \param right local evaluation
- * \param uLeft DOF evaluation on this side of \c intersection
- * \param uRight DOF evaluation on the other side of \c intersection
- * \param gLeft result for this side of \c intersection
- * \param gRight result for the other side of \c intersection
- * \return wave speed estimate (multiplied with the integration element of the intersection).
- * To estimate the time step |T|/wave is used
- */
- template <class LocalEvaluation>
- double gradientNumericalFlux(const LocalEvaluation& left,
- const LocalEvaluation& right,
- const RangeType& uLeft,
- const RangeType& uRight,
- GradientRangeType& gLeft,
- GradientRangeType& gRight,
- GradientJacobianType& gDiffLeft,
- GradientJacobianType& gDiffRight) const
- {
- const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
-
- // get factor for each side
- const double thetaLeft = theta( normal );
- const double thetaRight = 1.0 - thetaLeft;
-
- GradientJacobianType diffmatrix;
-
- double diffTimeStep = 0.0;
-
- // select left or right or average
- if( thetaLeft > 0 )
- {
- //TODO use diffusionTimeStep
- diffTimeStep = 0;
-
- /* central differences (might be suboptimal) */
- model_.jacobian(left, uLeft, diffmatrix );
-
- diffmatrix.mv(normal, gLeft );
- gLeft *= thetaLeft ;
- }
- else
- gLeft = 0;
-
- if( thetaRight > 0 )
- {
- //const double diffStepRight = 0;
-
- // right jacobian
- model_.jacobian( right, uRight, diffmatrix );
-
- diffmatrix.mv(normal, gRight);
-
- // add to flux
- gLeft.axpy( thetaRight, gRight );
-
- // diffTimeStep = std::max( diffTimeStep, diffStepRight );
- }
-
- // copy flux
- gRight = gLeft;
-
-#ifndef NDEBUG
- gDiffLeft = 0;
- gDiffRight = 0;
-#endif
-
- // upper bound for the next time step length
- return diffTimeStep; // * cflDiffinv_;
- }
-
-
- template <class LocalEvaluation>
- double gradientBoundaryFlux(const LocalEvaluation& left,
- const RangeType& uLeft,
- const RangeType& uBnd,
- GradientRangeType& gLeft,
- GradientJacobianType& gDiffLeft) const
- {
- const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
-
- // get factor for each side
- const double thetaLeft = theta( normal );
- const double thetaRight = 1.0 - thetaLeft;
-
- GradientJacobianType diffmatrix;
-
- // calculate uVal
- RangeType uVal( 0 );
-
- // select u from uLeft and uRight
- if( thetaLeft > 0 )
- uVal.axpy( thetaLeft , uLeft );
- if( thetaRight > 0 )
- uVal.axpy( thetaRight, uBnd );
-
- // compute jacobian of u
- model_.jacobian(left, uVal, diffmatrix );
-
- diffmatrix.mv(normal, gLeft);
-
-#ifndef NDEBUG
- gDiffLeft = 0;
-#endif
-
- //TODO use diffusionTimeStep
- const double diffTimeStep = 0; // * cflDiffinv_;
- return diffTimeStep;
- }
-
-
- /**
- * \brief flux function on interfaces between cells
- *
- * \param left local evaluation
- * \param right local evaluation
- * \param uLeft DOF evaluation on this side of \c intersection
- * \param uRight DOF evaluation on the other side of \c intersection
- * \param gLeft result for this side of \c intersection
- * \param gRight result for the other side of \c intersection
- * \return wave speed estimate (multiplied with the integration element of the intersection).
- * To estimate the time step |T|/wave is used
- */
- template <class LocalEvaluation>
- double numericalFlux(const LocalEvaluation& left,
- const LocalEvaluation& right,
- const RangeType& uLeft,
- const RangeType& uRight,
- const JacobianRangeType& sigmaLeft,
- const JacobianRangeType& sigmaRight,
- RangeType& gLeft,
- RangeType& gRight,
- JacobianRangeType& gDiffLeft, // not used here (only for primal passes)
- JacobianRangeType& gDiffRight )
- {
- const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
- const double faceLengthSqr = normal.two_norm2();
-
- /**********************************
- * Diffusion sigma Flux (Pass 2) *
- **********************************/
- JacobianRangeType diffmatrix;
-
- // diffusion for uLeft
- model_.diffusion( left, uLeft, sigmaLeft, diffmatrix);
-
- RangeType diffflux;
- diffmatrix.mv(normal, diffflux);
-
- // diffusion for uRight
- model_.diffusion( right, uRight, sigmaRight, diffmatrix);
- // add to diffflux
- diffmatrix.umv(normal, diffflux);
- diffflux *= 0.5;
-
- const double faceVolumeEstimate = dimensionFactor_ *
- (left.intersection().conforming() ? faceLengthSqr
- : (nonconformingFactor_ * faceLengthSqr));
-
- const double diffTimeLeft =
- model_.diffusionTimeStep( left, faceVolumeEstimate, uLeft );
-
- const double diffTimeRight =
- model_.diffusionTimeStep( right, faceVolumeEstimate, uRight );
-
- double diffTimeStep = std::max( diffTimeLeft, diffTimeRight );
-
- // add penalty factor
- const double factor = penalty_ * diffTimeStep ;
-
- RangeType jump( uLeft );
- jump -= uRight;
- diffflux.axpy(factor, jump);
-
- gLeft = diffflux;
- gRight = diffflux;
-
-#ifndef NDEBUG
- gDiffLeft = 0;
- gDiffRight = 0;
-#endif
-
- // timestep restict to diffusion timestep
- // WARNING: reconsider this
- diffTimeStep *= cflDiffinv_;
- return diffTimeStep;
- }
-
-
- /**
- * \brief same as numericalFlux() but for fluxes over boundary interfaces
- */
- template <class LocalEvaluation>
- double boundaryFlux(const LocalEvaluation& left,
- const RangeType& uLeft,
- const RangeType& uRight,
- const JacobianRangeType& sigmaLeft,
- RangeType& gLeft,
- JacobianRangeType& gDiffLeft )
- {
- // get local point
- const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
- const double faceLengthSqr = normal.two_norm2();
-
- /****************************/
- /* Diffusion (Pass 2) */
- /****************************/
- JacobianRangeType diffmatrix;
-
- model_.diffusion(left, uLeft, sigmaLeft, diffmatrix);
-
- diffmatrix.mv(normal, gLeft);
-
- const double faceVolumeEstimate = dimensionFactor_ * faceLengthSqr;
-
- double diffTimeStep =
- model_.diffusionTimeStep( left, faceVolumeEstimate, uLeft );
-
- // add penalty term
- const double factor = penalty_ * diffTimeStep;
-
- RangeType jump( uLeft );
- jump -= uRight;
- gLeft.axpy(factor,jump);
-
-#ifndef NDEBUG
- gDiffLeft = 0;
-#endif
-
- diffTimeStep *= cflDiffinv_;
- return diffTimeStep;
- }
- protected:
- const double penalty_;
- const bool penaltyTerm_;
-
- }; // end LDGAverageDiffusionFlux
-
-} // end namespace
-} // end namespace
-#endif
diff --git a/dune/fem-dg/operator/fluxes/diffusion/fluxes.hh b/dune/fem-dg/operator/fluxes/diffusion/fluxes.hh
index 12ad3073..3a4494e5 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/fluxes.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/fluxes.hh
@@ -4,7 +4,6 @@
#include "fluxbase.hh"
#include "dgprimalfluxes.hh"
#include "ldgflux.hh"
-#include "averageflux.hh"
namespace Dune
{
@@ -227,9 +226,9 @@ namespace Fem
template <class DiscreteFunctionSpaceImp,
class Model>
class DGDualDiffusionFlux< DiscreteFunctionSpaceImp, Model, DualDiffusionFlux::Enum::average >
- : public LDGAverageDiffusionFlux< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
+ : public LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
{
- typedef LDGAverageDiffusionFlux< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
+ typedef LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
BaseType;
public:
@@ -243,7 +242,7 @@ namespace Fem
DGDualDiffusionFlux( GridPartType& gridPart,
const Model& model,
const ParameterType& parameters = ParameterType() )
- : BaseType( gridPart, model, parameters )
+ : BaseType( gridPart, model, parameters, DualDiffusionFlux::Enum::average )
{
}
};
@@ -251,9 +250,9 @@ namespace Fem
template <class DiscreteFunctionSpaceImp,
class Model>
class DGDualDiffusionFlux< DiscreteFunctionSpaceImp, Model, DualDiffusionFlux::Enum::ldg >
- : public LDGDiffusionFlux< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
+ : public LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
{
- typedef LDGDiffusionFlux< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
+ typedef LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
BaseType;
public:
@@ -267,18 +266,35 @@ namespace Fem
DGDualDiffusionFlux( GridPartType& gridPart,
const Model& model,
const ParameterType& parameters = ParameterType() )
- : BaseType( gridPart, model, parameters )
+ : BaseType( gridPart, model, parameters, DualDiffusionFlux::Enum::ldg )
{
}
};
- //TODO implement general case...
template <class DiscreteFunctionSpaceImp,
class Model>
class DGDualDiffusionFlux< DiscreteFunctionSpaceImp, Model, DualDiffusionFlux::Enum::general >
- {};
+ : public LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
+ {
+ typedef LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
+ BaseType;
+
+ public:
+ typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType;
+ typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
+ typedef typename BaseType::ParameterType ParameterType;
+ /**
+ * \brief constructor reading parameters
+ */
+ DGDualDiffusionFlux( GridPartType& gridPart,
+ const Model& model,
+ const ParameterType& parameters = ParameterType() )
+ : BaseType( gridPart, model, parameters, DualDiffusionFlux::Enum::general )
+ {
+ }
+ };
-} // end namespace
-} // end namespace
+} // end namespace Fem
+} // end namespace Dune
#endif
diff --git a/dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh b/dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh
index ba794bf3..a6b0a2ca 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh
@@ -1,5 +1,5 @@
-#ifndef DUNE_FEM_DG_LDGFLUX_HH
-#define DUNE_FEM_DG_LDGFLUX_HH
+#ifndef DUNE_FEM_DG_LDGAVERAGEFLUX_HH
+#define DUNE_FEM_DG_LDGAVERAGEFLUX_HH
// local includes
#include "fluxbase.hh"
@@ -10,17 +10,17 @@ namespace Fem
{
/**
- * \brief diffusion flux
+ * \brief Implementation of Bassi-Rebay 1 and Local DG fluxes.
*
* \ingroup DiffusionFluxes
*/
template <class DiscreteFunctionSpaceImp,
class ModelImp,
class FluxParameterImp>
- class LDGDiffusionFlux :
- public DGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp >
+ class LDGDiffusionFluxImpl :
+ public LDGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp >
{
- typedef DGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp > BaseType;
+ typedef LDGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp > BaseType;
public:
typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType;
@@ -46,7 +46,8 @@ namespace Fem
typedef typename BaseType :: ParameterType ParameterType;
// type of gradient space
- typedef typename BaseType :: DiscreteGradientSpaceType DiscreteGradientSpaceType;
+ typedef typename DiscreteFunctionSpaceType ::
+ template ToNewDimRange< dimGradRange > :: Type DiscreteGradientSpaceType;
typedef typename DiscreteGradientSpaceType :: RangeType GradientRangeType;
typedef typename DiscreteGradientSpaceType :: JacobianRangeType GradientJacobianType;
@@ -54,394 +55,309 @@ namespace Fem
// jacobians of the functions do not have to be evaluated for this flux
enum { evaluateJacobian = false };
- private:
- // no copying
- LDGDiffusionFlux(const LDGDiffusionFlux& other);
protected:
- using BaseType :: determineDirection;
- using BaseType :: model_;
- using BaseType :: cflDiffinv_;
- using BaseType :: dimensionFactor_;
- using BaseType :: nonconformingFactor_;
- using BaseType :: parameter;
+ typedef typename BaseType::IdEnum EnumType;
+
+ using BaseType::determineDirection;
+ using BaseType::model_;
+ using BaseType::cflDiffinv_;
+ using BaseType::numericalFlux ;
+ using BaseType::parameter ;
+ using BaseType::nonconformingFactor_;
+ using BaseType::dimensionFactor_;
public:
/**
* \brief constructor
*/
- LDGDiffusionFlux(GridPartType& gridPart,
- const ModelImp& model,
- const ParameterType& parameters = ParameterType() )
- : BaseType( model, true, parameters ),
+ LDGDiffusionFluxImpl(GridPartType& gridPart,
+ const ModelImp& mod,
+ const ParameterType& param,
+ const EnumType staticMethod ) :
+ BaseType( gridPart, mod, param ),
+ method_( staticMethod == EnumType::general ? param.getMethod() : staticMethod ),
penalty_( parameter().penalty() ),
// Set CFL number for penalty term (compare diffusion in first pass)
penaltyTerm_( std::abs( penalty_ ) > 0 )
{
if( Fem::Parameter::verbose () )
{
- std::cout << "LDGDiffusionFlux: penalty = " << penalty_ << std::endl;
+ std::cout << "LDGDiffusionFluxImpl: penalty = " << penalty_ << std::endl;
}
}
+ LDGDiffusionFluxImpl(const LDGDiffusionFluxImpl& other)
+ : BaseType( other ),
+ method_( other.method_ ),
+ penalty_( other.penalty_ ),
+ penaltyTerm_( other.penaltyTerm_ )
+ {}
+
//! returns true if lifting has to be calculated
- bool hasLifting () const { return false; }
+ const bool hasLifting () const { return false; }
- void diffusionFluxPenalty( std::ostream& out ) const
+ protected:
+ enum { realLDG = true };
+ double theta( const DomainType& normal ) const
{
- out <<penalty_;
+ if( method_ == EnumType::ldg )
+ {
+ // LDG thete is 1 or 0
+ if( determineDirection( normal ) )
+ return 1.0;
+ else
+ return 0.0;
+ }
+ else
+ {
+ assert( method_ == EnumType::br1 );
+ // Average fluxes (Bassi-Rebay 1)
+ return 0.5;
+ }
}
public:
/**
* \brief flux function on interfaces between cells
*
- * @param intersection intersection
- * @param time current time given by TimeProvider
- * @param x coordinate of required evaluation local to \c intersection
- * @param uLeft DOF evaluation on this side of \c intersection
- * @param uRight DOF evaluation on the other side of \c intersection
- * @param gLeft result for this side of \c intersection
- * @param gRight result for the other side of \c intersection
+ * \param left local evaluation
+ * \param right local evaluation
+ * \param uLeft DOF evaluation on this side of \c intersection
+ * \param uRight DOF evaluation on the other side of \c intersection
+ * \param gLeft result for this side of \c intersection
+ * \param gRight result for the other side of \c intersection
* \return wave speed estimate (multiplied with the integration element of the intersection).
* To estimate the time step |T|/wave is used
*/
- template <class QuadratureImp>
- double gradientNumericalFlux(
- const Intersection& intersection,
- const EntityType& inside,
- const EntityType& outside,
- const double time,
- const QuadratureImp& faceQuadInner,
- const QuadratureImp& faceQuadOuter,
- const int quadPoint,
- const RangeType& uLeft,
- const RangeType& uRight,
- GradientRangeType& gLeft,
- GradientRangeType& gRight,
- GradientJacobianType& gDiffLeft,
- GradientJacobianType& gDiffRight) const
+ template <class LocalEvaluation>
+ double gradientNumericalFlux(const LocalEvaluation& left,
+ const LocalEvaluation& right,
+ const RangeType& uLeft,
+ const RangeType& uRight,
+ GradientRangeType& gLeft,
+ GradientRangeType& gRight,
+ GradientJacobianType& gDiffLeft,
+ GradientJacobianType& gDiffRight) const
{
- const FaceDomainType& x = faceQuadInner.localPoint( quadPoint );
- const DomainType normal = intersection.integrationOuterNormal( x );
+ const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
- // determine side
- const bool useInterior = determineDirection( normal );
+ // get factor for each side
+ const double thetaLeft = theta( normal );
+ const double thetaRight = 1.0 - thetaLeft;
GradientJacobianType diffmatrix;
- if( useInterior )
+ double diffTimeStep = 0.0;
+
+ // select left or right or average
+ if( thetaLeft > 0 )
{
- model_.jacobian(inside, /* inside entity */
- time, /* for time dependent diffusion */
- faceQuadInner.point( quadPoint ), /* inside point on intersection */
- uLeft, /* { u_(x^-) } */
- diffmatrix /* return diffusion tensor */
- );
+ //TODO use diffusionTimeStep
+ diffTimeStep = 0;
+
+ /* central differences (might be suboptimal) */
+ model_.jacobian(left, uLeft, diffmatrix );
+
+ diffmatrix.mv(normal, gLeft );
+ gLeft *= thetaLeft ;
}
else
+ gLeft = 0;
+
+ if( thetaRight > 0 )
{
- model_.jacobian(outside, /* outside entity */
- time, /* for time dependent diffusion */
- faceQuadOuter.point( quadPoint ), /* outside point on intersection */
- uRight, /* { u_(x^+) } */
- diffmatrix /* return diffusion tensor */
- );
- }
+ //const double diffStepRight = 0;
- // mutliply with normal
- diffmatrix.mv(normal, gLeft);
+ // right jacobian
+ model_.jacobian( right, uRight, diffmatrix );
+
+ diffmatrix.mv(normal, gRight);
+
+ // add to flux
+ gLeft.axpy( thetaRight, gRight );
+
+ // diffTimeStep = std::max( diffTimeStep, diffStepRight );
+ }
// copy flux
gRight = gLeft;
+#ifndef NDEBUG
gDiffLeft = 0;
gDiffRight = 0;
+#endif
- // time step is set in 2nd pass
- return 0.0;
+ // upper bound for the next time step length
+ return diffTimeStep; // * cflDiffinv_;
}
- /*
- * \brief numerical flux for u given as u_h
- */
- template <class QuadratureImp>
- double gradientBoundaryFlux(const Intersection& intersection,
- const EntityType& inside,
- const double time,
- const QuadratureImp& faceQuadInner,
- const int quadPoint,
+
+ template <class LocalEvaluation>
+ double gradientBoundaryFlux(const LocalEvaluation& left,
const RangeType& uLeft,
const RangeType& uBnd,
GradientRangeType& gLeft,
GradientJacobianType& gDiffLeft) const
{
- const FaceDomainType& x = faceQuadInner.localPoint( quadPoint );
- const DomainType normal = intersection.integrationOuterNormal(x);
- const DomainType& xglInside = faceQuadInner.point( quadPoint );
+ const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
+
+ // get factor for each side
+ const double thetaLeft = theta( normal );
+ const double thetaRight = 1.0 - thetaLeft;
+
GradientJacobianType diffmatrix;
- // determine side
- //const bool useInterior = determineDirection( normal );
+ // calculate uVal
+ RangeType uVal( 0 );
- // get apropriate value
- //const RangeType& uVal = ( useInterior ) ? uLeft : uBnd ;
+ // select u from uLeft and uRight
+ if( thetaLeft > 0 )
+ uVal.axpy( thetaLeft , uLeft );
+ if( thetaRight > 0 )
+ uVal.axpy( thetaRight, uBnd );
- // for the numerical diffusion flux \tilde u
- // one uses \tilde u = g_D where g_D is Dirichlet boundary data
- model_.jacobian( inside,
- time,
- xglInside,
- uBnd,
- diffmatrix
- );
+ // compute jacobian of u
+ model_.jacobian(left, uVal, diffmatrix );
- // apply normal
diffmatrix.mv(normal, gLeft);
- // gDiffLeft is not needed for dual formalation (LDG)
+#ifndef NDEBUG
gDiffLeft = 0;
+#endif
- // time step is set in 2nd pass
- return 0.0;
+ //TODO use diffusionTimeStep
+ const double diffTimeStep = 0; // * cflDiffinv_;
+ return diffTimeStep;
}
/**
* \brief flux function on interfaces between cells
*
- * @param intersection intersection
- * @param time current time given by TimeProvider
- * @param x coordinate of required evaluation local to \c intersection
- * @param uLeft DOF evaluation on this side of \c intersection
- * @param uRight DOF evaluation on the other side of \c intersection
- * @param gLeft result for this side of \c intersection
- * @param gRight result for the other side of \c intersection
+ * \param left local evaluation
+ * \param right local evaluation
+ * \param uLeft DOF evaluation on this side of \c intersection
+ * \param uRight DOF evaluation on the other side of \c intersection
+ * \param gLeft result for this side of \c intersection
+ * \param gRight result for the other side of \c intersection
* \return wave speed estimate (multiplied with the integration element of the intersection).
* To estimate the time step |T|/wave is used
*/
- template <class DiscreteModelImp,
- class QuadratureImp>
- inline double numericalFlux(const Intersection& intersection,
- const DiscreteModelImp& discreteModel,
- const double time,
- const QuadratureImp& faceQuadInner,
- const QuadratureImp& faceQuadOuter,
- const int quadPoint,
- const RangeType& uLeft,
- const RangeType& uRight,
- const GradientRangeType& sigmaLeft,
- const GradientRangeType& sigmaRight,
- RangeType& gLeft,
- RangeType& gRight,
- JacobianRangeType& gDiffLeft, // not used here (only for primal passes)
- JacobianRangeType& gDiffRight ) const
- {
- Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType> jacLeft( sigmaLeft );
- Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType> jacRight( sigmaRight );
- return numericalFlux( intersection, discreteModel, time,
- faceQuadInner, faceQuadOuter, quadPoint,
- uLeft, uRight, jacLeft, jacRight,
- gLeft, gRight, gDiffLeft, gDiffRight );
- }
-
-
- template< class DiscreteModelImp,
- class QuadratureImp, class JacobianRangeTypeImp >
- double numericalFlux(const Intersection& intersection,
- const DiscreteModelImp& discreteModel,
- const double time,
- const QuadratureImp& faceQuadInner,
- const QuadratureImp& faceQuadOuter,
- const int quadPoint,
+ template <class LocalEvaluation>
+ double numericalFlux(const LocalEvaluation& left,
+ const LocalEvaluation& right,
const RangeType& uLeft,
const RangeType& uRight,
- const JacobianRangeTypeImp& sigmaLeft,
- const JacobianRangeTypeImp& sigmaRight,
+ const JacobianRangeType& sigmaLeft,
+ const JacobianRangeType& sigmaRight,
RangeType& gLeft,
RangeType& gRight,
JacobianRangeType& gDiffLeft, // not used here (only for primal passes)
- JacobianRangeType& gDiffRight ) const
+ JacobianRangeType& gDiffRight )
{
- const FaceDomainType& x = faceQuadInner.localPoint( quadPoint );
- const DomainType normal = intersection.integrationOuterNormal(x);
-
- const EntityType& inside = discreteModel.inside();
- const EntityType& outside = discreteModel.outside();
+ const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
+ const double faceLengthSqr = normal.two_norm2();
- /**********************************
- * Diffusion sigma Flux (Pass 2) *
- **********************************/
+ /**********************************
+ * Diffusion sigma Flux (Pass 2) *
+ **********************************/
JacobianRangeType diffmatrix;
- // determine side (needs to be opposite of above)
- const bool useExterior = ! determineDirection( normal );
+ // diffusion for uLeft
+ model_.diffusion( left, uLeft, sigmaLeft, diffmatrix);
+ RangeType diffflux;
+ diffmatrix.mv(normal, diffflux);
- if( useExterior )
- {
- model_.diffusion( inside, time,
- faceQuadInner.point( quadPoint ),
- uLeft, sigmaLeft, diffmatrix);
- }
- else
- {
- model_.diffusion( outside, time,
- faceQuadOuter.point( quadPoint ),
- uRight, sigmaRight, diffmatrix);
- }
-
- // apply normal
- diffmatrix.mv(normal, gLeft);
-
-
- //////////////////////////////////////////////////////////
- //
- // --Time step calculation
- //
- //////////////////////////////////////////////////////////
- const double faceLengthSqr = normal.two_norm2();
+ // diffusion for uRight
+ model_.diffusion( right, uRight, sigmaRight, diffmatrix);
+ // add to diffflux
+ diffmatrix.umv(normal, diffflux);
+ diffflux *= 0.5;
const double faceVolumeEstimate = dimensionFactor_ *
- ( intersection.conforming() ) ? faceLengthSqr : ( nonconformingFactor_ * faceLengthSqr );
+ (left.intersection().conforming() ? faceLengthSqr
+ : (nonconformingFactor_ * faceLengthSqr));
const double diffTimeLeft =
- model_.diffusionTimeStep( intersection,
- discreteModel.enVolume(),
- faceVolumeEstimate,
- time, x, uLeft );
+ model_.diffusionTimeStep( left, faceVolumeEstimate, uLeft );
const double diffTimeRight =
- model_.diffusionTimeStep( intersection,
- discreteModel.nbVolume(),
- faceVolumeEstimate,
- time, x, uRight );
-
- //////////////////////////////////////////////////////////
- //
- // --Penalty Term
- //
- //////////////////////////////////////////////////////////
-
- // take minimum to proceed
- const double diffTimeStep = std::max( diffTimeLeft, diffTimeRight );
- if( penaltyTerm_ )
- {
- // add penalty term ( enVolume() is available since we derive from
- // DiscreteModelDefaultWithInsideOutside)
- const double factor = penalty_ * diffTimeStep ;
+ model_.diffusionTimeStep( right, faceVolumeEstimate, uRight );
- RangeType jump( uLeft );
- jump -= uRight;
- gLeft.axpy(factor, jump);
- }
+ double diffTimeStep = std::max( diffTimeLeft, diffTimeRight );
+
+ // add penalty factor
+ const double factor = penalty_ * diffTimeStep ;
- gRight = gLeft ;
+ RangeType jump( uLeft );
+ jump -= uRight;
+ diffflux.axpy(factor, jump);
- // gDiffLeft should be 0 in case of LDG
+ gLeft = diffflux;
+ gRight = diffflux;
+
+#ifndef NDEBUG
gDiffLeft = 0;
gDiffRight = 0;
+#endif
// timestep restict to diffusion timestep
// WARNING: reconsider this
- return diffTimeStep * cflDiffinv_;
+ diffTimeStep *= cflDiffinv_;
+ return diffTimeStep;
}
/**
* \brief same as numericalFlux() but for fluxes over boundary interfaces
*/
- template <class DiscreteModelImp,
- class QuadratureImp>
- inline double boundaryFlux(const Intersection& intersection,
- const DiscreteModelImp& discreteModel,
- const double time,
- const QuadratureImp& faceQuadInner,
- const int quadPoint,
+ template <class LocalEvaluation>
+ double boundaryFlux(const LocalEvaluation& left,
const RangeType& uLeft,
const RangeType& uRight,
- const GradientRangeType& sigmaLeft,
+ const JacobianRangeType& sigmaLeft,
RangeType& gLeft,
- JacobianRangeType& gDiffLeft) const /*@LST0E@*/
- {
- Fem::FieldMatrixConverter< GradientRangeType, JacobianRangeType> jacLeft( sigmaLeft );
- return boundaryFlux( intersection, discreteModel, time,
- faceQuadInner, quadPoint, uLeft, uRight,
- jacLeft, gLeft, gDiffLeft );
- }
-
- /**
- * \brief same as numericalFlux() but for fluxes over boundary interfaces
- */
- template< class DiscreteModelImp,
- class QuadratureImp, class JacobianRangeTypeImp >
- double boundaryFlux(const Intersection& intersection,
- const DiscreteModelImp& discreteModel,
- const double time,
- const QuadratureImp& faceQuadInner,
- const int quadPoint,
- const RangeType& uLeft,
- const RangeType& uRight,
- const JacobianRangeTypeImp& sigmaLeft,
- RangeType& gLeft,
- JacobianRangeType& gDiffLeft) const /*@LST0E@*/
+ JacobianRangeType& gDiffLeft )
{
// get local point
- const FaceDomainType& x = faceQuadInner.localPoint( quadPoint );
- const DomainType normal = intersection.integrationOuterNormal(x);
-
- const EntityType& inside = discreteModel.inside();
+ const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
+ const double faceLengthSqr = normal.two_norm2();
/****************************/
/* Diffusion (Pass 2) */
/****************************/
JacobianRangeType diffmatrix;
- model_.diffusion(inside, time,
- faceQuadInner.point( quadPoint ),
- uLeft, sigmaLeft, diffmatrix);
+ model_.diffusion(left, uLeft, sigmaLeft, diffmatrix);
diffmatrix.mv(normal, gLeft);
+ const double faceVolumeEstimate = dimensionFactor_ * faceLengthSqr;
- //////////////////////////////////////////////////////////
- //
- // --Time step calculation
- //
- //////////////////////////////////////////////////////////
- const double faceVolumeEstimate = normal.two_norm2();
-
- const double diffTimeStep =
- model_.diffusionTimeStep( intersection,
- discreteModel.enVolume(),
- faceVolumeEstimate,
- time, x, uLeft );
-
- //////////////////////////////////////////////////////////
- //
- // --Penalty Term
- //
- //////////////////////////////////////////////////////////
- if( penaltyTerm_ )
- {
- // add penalty term
- const double factor = penalty_ * diffTimeStep;
+ double diffTimeStep =
+ model_.diffusionTimeStep( left, faceVolumeEstimate, uLeft );
- RangeType jump( uLeft );
- jump -= uRight;
- gLeft.axpy(factor, jump);
- }
+ // add penalty term
+ const double factor = penalty_ * diffTimeStep;
+
+ RangeType jump( uLeft );
+ jump -= uRight;
+ gLeft.axpy(factor,jump);
- // gDiffLeft should be 0 in case of LDG
+#ifndef NDEBUG
gDiffLeft = 0;
+#endif
- return diffTimeStep * cflDiffinv_;
+ diffTimeStep *= cflDiffinv_;
+ return diffTimeStep;
}
-
protected:
+ const EnumType method_;
const double penalty_;
const bool penaltyTerm_;
- }; // end LDGDiffusionFlux
+ }; // end LDGDiffusionFluxImpl
} // end namespace
} // end namespace
diff --git a/dune/fem-dg/operator/fluxes/diffusion/parameters.hh b/dune/fem-dg/operator/fluxes/diffusion/parameters.hh
index 25231908..6605282a 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/parameters.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/parameters.hh
@@ -40,7 +40,7 @@ namespace Fem
//! Contains all known enums for primal diffusion fluxes which can be chosen via parameter file.
const Enum _enums[] = { Enum::cdg2, Enum::cdg, Enum::br2, Enum::ip, Enum::nipg, Enum::bo };
- //! Contains all known names of primal diffusion fluxes which can be chosen via parameter file.
+ //! Contains all known names of primal DG diffusion fluxes which can be chosen via parameter file.
const std::string _strings[] = { "CDG2", "CDG" , "BR2", "IP" , "NIPG", "BO" };
//! Number of known primal diffusion fluxes which can be chosen via parameter file.
static const int _size = 6;
@@ -50,15 +50,15 @@ namespace Fem
namespace DualDiffusionFlux
{
/**
- * \brief Enum of all known primal diffusion flux implementations.
+ * \brief Enum of all known local DG diffusion flux implementations.
*
* \ingroup FemDGParameter
*/
enum class Enum
{
- //! average theta flux.
+ //! Bassi-Rebay flux.
average,
- //! ldg flux.
+ //! Local DG flux.
ldg,
//! general flux: Parameter selection is done via parameter file!
general,
--
GitLab