From 8fe670922d87afc1ad8f7d02ee7e51740c439e06 Mon Sep 17 00:00:00 2001
From: Robert Kloefkorn <robertk@ucar.edu>
Date: Fri, 6 Dec 2013 12:27:27 +0100
Subject: [PATCH] added advection+diffusion problems.
---
dune/fem-dg/test/advdiff/Makefile.am | 69 +++
dune/fem-dg/test/advdiff/deformationalflow.hh | 183 ++++++
dune/fem-dg/test/advdiff/models.hh | 566 ++++++++++++++++++
dune/fem-dg/test/advdiff/parameter | 56 ++
dune/fem-dg/test/advdiff/problem.hh | 255 ++++++++
.../test/advdiff/problemQuasiHeatEqn.hh | 166 +++++
dune/fem-dg/test/advdiff/problemcreator.hh | 99 +++
7 files changed, 1394 insertions(+)
create mode 100644 dune/fem-dg/test/advdiff/Makefile.am
create mode 100644 dune/fem-dg/test/advdiff/deformationalflow.hh
create mode 100644 dune/fem-dg/test/advdiff/models.hh
create mode 100644 dune/fem-dg/test/advdiff/parameter
create mode 100644 dune/fem-dg/test/advdiff/problem.hh
create mode 100644 dune/fem-dg/test/advdiff/problemQuasiHeatEqn.hh
create mode 100644 dune/fem-dg/test/advdiff/problemcreator.hh
diff --git a/dune/fem-dg/test/advdiff/Makefile.am b/dune/fem-dg/test/advdiff/Makefile.am
new file mode 100644
index 00000000..7e95ba12
--- /dev/null
+++ b/dune/fem-dg/test/advdiff/Makefile.am
@@ -0,0 +1,69 @@
+# install these headers
+advdiffdir=$(includedir)/test/advdiff
+advdiff_HEADERS = models.hh problemcreator.hh \
+problem.hh problemQuasiHeatEqn.hh
+
+LDADD = $(ALL_PKG_LDFLAGS) $(ALL_PKG_LIBS) $(LOCAL_LIBS) $(DUNEMPILDFLAGS) $(DUNEMPILIBS)
+
+BASEDIR=../../main
+
+advdiff_SOURCES = $(BASEDIR)/main.cc $(BASEDIR)/main_1.cc $(BASEDIR)/main_2.cc \
+ $(BASEDIR)/main_0.cc $(BASEDIR)/main_3.cc
+advdiffonep_SOURCES = $(BASEDIR)/main.cc $(BASEDIR)/main_pol.cc
+advdiff12_SOURCES = $(BASEDIR)/main.cc $(BASEDIR)/main_1.cc $(BASEDIR)/main_2.cc
+
+#GRIDTYPE = ALUGRID_SIMPLEX
+#GRIDTYPE = YASPGRID
+#GRIDTYPE=PARALLELGRID_ALUGRID_SIMPLEX
+#GRIDTYPE=CARTESIANGRID_ALUGRID_CUBE
+GRIDTYPE = SPGRID
+POLORDER = 2
+GRIDDIM = 2
+DIMRANGE = 1
+FLUX = 1 # LLF 1, HLL 2
+
+#USE_OMP=-fopenmp
+
+# INFO SPACE OPERATOR:
+# 1. define PRIMALDG for IP, BR2, CDG, CDG2
+# 2. define DUALDG for LDG
+# INFO TRACK LIFTING:
+# 1. define LOCALDEBUG to calculate \sum_e\int_Omega(r_e*l_e) and
+# \sum_e\int_Omega(r_e*l_e). They will be output to std::cout from the Stepper
+# INFO TESTOPERATOR:
+# 1. define TESTOPERATOR for linear advdiff equation to output various
+# information on space operator
+check_PROGRAMS = advdiff advdiffonep advdiff12 # quasiadvdiff quasiadvdiff12 plaplace plaplace12
+
+AM_CPPFLAGS = $(USE_OMP) -I../../problems/advdiff $(ALL_PKG_CPPFLAGS) -DGRIDDIM=$(GRIDDIM) \
+ -D$(GRIDTYPE) $(DUNEMPICPPFLAGS) -DFEMTIMER \
+ -DDIMRANGE=$(DIMRANGE) -DFLUX=$(FLUX) -DPRIMALDG -DUSE_ASSERT_THROW
+AM_LDFLAGS = $(ALL_PKG_LDFLAGS) $(LAPACK_LDFLAGS) $(USE_OMP)
+
+advdiff_CPPFLAGS = $(AM_CPPFLAGS)
+advdiffonep_CPPFLAGS = $(advdiff_CPPFLAGS) -DONLY_ONE_P -DPOLORDER=$(POLORDER)
+advdiff12_CPPFLAGS = $(advdiff_CPPFLAGS) -DONLY_P1_P2
+
+DISTCHECK_CONFIGURE_FLAGS = DOXYGEN="true"
+
+EXTRA_DIST = parameter
+
+CLEANFILES=manager.*.log eoc_*.tex
+
+PROG=advdiff12
+# codegeneration
+generatecode:
+ $(MAKE) -i clean
+ $(MAKE) CXXFLAGS="-O0 -Wall -DNDEBUG -DBASEFUNCTIONSET_CODEGEN_GENERATE" $(PROG)
+ ./$(PROG) femhowto.maximaltimesteps:1 fem.io.outputformat:none codegenparameter
+
+# compile fast code
+compilecode:
+ $(MAKE) clean
+ $(MAKE) CXXFLAGS="$(CXXFLAGS) -DUSE_BASEFUNCTIONSET_CODEGEN" $(PROG)
+
+codegen:
+ $(MAKE) generatecode
+ $(MAKE) compilecode
+
+include $(top_srcdir)/am/global-rules
diff --git a/dune/fem-dg/test/advdiff/deformationalflow.hh b/dune/fem-dg/test/advdiff/deformationalflow.hh
new file mode 100644
index 00000000..f10b73bb
--- /dev/null
+++ b/dune/fem-dg/test/advdiff/deformationalflow.hh
@@ -0,0 +1,183 @@
+#ifndef DUNE_DEFORMATIONALFLOW_HH
+#define DUNE_DEFORMATIONALFLOW_HH
+
+// dune-fem includes
+#include <dune/fem/io/parameter.hh>
+#include <dune/fem/space/common/functionspace.hh>
+
+// local includes
+#include <dune/fem-dg/models/defaultprobleminterfaces.hh>
+
+namespace Dune {
+
+/**
+ * @brief describes the initial and exact solution of the advection-diffusion model
+ * for given constant velocity vector v=(v1,v2)
+ *
+ * \f[u(x,y,z,t):=\displaystyle{\sum_{i=0}^{1}} T_i(t) \cdot X_i(x) \cdot
+ * Y_i(y) \cdot Z_i(z)\f]
+ *
+ * with
+ *
+ * \f{eqnarray*}{
+ * T_0(t) &:=& e^{-\varepsilon t \pi^2 (2^2 + 1^2 + 1.3^2 )} \\
+ * X_0(x) &:=& 0.6\cdot \cos(2\pi (x-v_1t)) + 0.8\cdot \sin(2\pi (x-v_1t)) \\
+ * Y_0(y) &:=& 1.2\cdot \cos(1\pi (y-v_2t)) + 0.4\cdot \sin(1\pi (y-v_2t)) \\
+ * Z_0(z) &:=& 0.1\cdot \cos(1.3\pi (z-v_3t)) - 0.4\cdot \sin(1.3\pi (z-v_3t)) \\
+ * T_1(t) &:=& e^{-\varepsilon t \pi^2 (0.7^2 + 0.5^2 + 0.1^2 )} \\
+ * X_1(x) &:=& 0.9\cdot \cos(0.7\pi (x-v_1t)) + 0.2\cdot \sin(0.7\pi (x-v_1t)) \\
+ * Y_1(y) &:=& 0.3\cdot \cos(0.5\pi (y-v_2t)) + 0.1\cdot \sin(0.5\pi (y-v_2t))
+ * Z_1(z) &:=& -0.3\cdot \cos(0.1\pi (z-v_3t)) + 0.2\cdot \sin(0.1\pi (z-v_3t)) \\
+ * \f}
+ *
+ * This is a solution of the AdvectionDiffusionModel for \f$g_D = u|_{\partial
+ * \Omega}\f$.
+ *
+ */
+template <class GridType> /*@LST0S@*/
+struct DeformationalFlow : public EvolutionProblemInterface<
+ Fem::FunctionSpace< double, double, GridType::dimension, DIMRANGE>,
+ false >
+{ /*@LST0E@*/
+public:
+ typedef EvolutionProblemInterface<
+ Fem::FunctionSpace< double, double,
+ GridType::dimension, DIMRANGE >,
+ false > BaseType;
+
+ enum{ dimDomain = BaseType :: dimDomain };
+ typedef typename BaseType :: DomainType DomainType;
+ typedef typename BaseType :: RangeType RangeType;
+ typedef typename BaseType :: JacobianRangeType JacobianRangeType;
+
+ typedef Fem::Parameter ParameterType;
+
+ /**
+ * @brief define problem parameters
+ */
+ DeformationalFlow () : /*@LST0S@*/
+ BaseType () ,
+ center_( 0.5 ),
+ startTime_( ParameterType::getValue<double>("femhowto.startTime",0.0) ),
+ endTime_( ParameterType::getValue<double>("femhowto.endTime") ),
+ epsilon_( ParameterType::getValue<double>("femhowto.epsilon", 0.0 ) )
+ {
+ myName = "DeformationalFlow";
+ }
+
+ //! this problem has no source term
+ bool hasStiffSource() const { return false; }
+ bool hasNonStiffSource() const { return false; }
+
+ double stiffSource(const DomainType& arg,
+ const double t,
+ const RangeType& u,
+ RangeType& res) const
+ {
+ return 0.0;
+ }
+
+ double nonStiffSource(const DomainType& arg,
+ const double t,
+ const RangeType& u,
+ RangeType& res) const
+ {
+ return 0.0;
+ }
+
+ double diffusion( const RangeType& u, const JacobianRangeType& gradU ) const
+ {
+ return epsilon_;
+ }
+
+ //! return start time
+ double startTime() const { return startTime_; }
+
+ //! return start time
+ double epsilon() const { return epsilon_; }
+
+ /**
+ * @brief getter for the velocity
+ */
+ void velocity(const DomainType& x, const double time, DomainType& v) const
+ {
+ DomainType p( x );
+ p -= center_;
+
+ const double r = p.two_norm();
+ const double f1 = std::pow( (4.0*r), 6.0 );
+ const double factor = (1.0 - f1) / (1.0 + f1);
+
+ const double theta = std::atan( p[ 1 ]/p[ 0 ] );
+ const double utheta = 4.0 * M_PI * r / endTime_ * ( 1.0 - ( std::cos( 2.0 * M_PI * time / endTime_ ) * factor )) ;
+
+ v[ 0 ] = utheta * std::sin( theta );
+ v[ 1 ] = -utheta * std::cos( theta );
+ }
+
+ /**
+ * @brief evaluates \f$ u_0(x) \f$
+ */
+ void evaluate(const DomainType& arg, RangeType& res) const /*@LST0S@@LST0E@*/
+ {
+ evaluate(arg, startTime_, res);
+ }
+
+ /**
+ * @brief evaluate exact solution
+ */
+ void evaluate(const DomainType& x, const double t, RangeType& res) const /*@LST0S@@LST0E@*/
+ {
+ DomainType c2( 0.5 );
+ c2[ 0 ] = 0.3;
+
+ DomainType p( x );
+ p -= c2;
+
+ const double rtilde = p.two_norm() * 5.0;
+ if( rtilde <= 1.0 )
+ {
+ double rc = 0.5 * ( 1.0 + std::cos( M_PI * rtilde) );
+ res = rc * rc ;
+ }
+ else
+ {
+ res = 0.0;
+ }
+ }
+
+ /**
+ * @brief latex output for EocOutput
+ */
+ std::string description() const
+ {
+ std::ostringstream ofs;
+
+ ofs << "Problem: " << myName
+ << ", Epsilon: " << epsilon_ ;
+ return ofs.str();
+ }
+
+ /* \brief finalize the simulation using the calculated numerical
+ * solution u for this problem
+ *
+ * \param[in] variablesToOutput Numerical solution in the suitably chosen variables
+ * \param[in] eocloop Specific EOC loop
+ */
+ template< class DiscreteFunctionType >
+ void finalizeSimulation( DiscreteFunctionType& variablesToOutput,
+ const int eocloop) const
+ {}
+
+private:
+ DomainType center_;
+ const double startTime_;
+ const double endTime_;
+public:
+ const double epsilon_;
+ std::string myName;
+};
+
+}
+#endif /*DUNE_PROBLEM_HH__*/
+
diff --git a/dune/fem-dg/test/advdiff/models.hh b/dune/fem-dg/test/advdiff/models.hh
new file mode 100644
index 00000000..d429fa9f
--- /dev/null
+++ b/dune/fem-dg/test/advdiff/models.hh
@@ -0,0 +1,566 @@
+#ifndef HEAT_MODELS_HH
+#define HEAT_MODELS_HH
+
+#include <dune/fem/version.hh>
+#include <dune/fem/misc/fmatrixconverter.hh>
+#include <dune/fem/io/parameter.hh>
+
+
+#include <dune/fem-dg/operator/limiter/limitpass.hh>
+
+// local includes
+#include <dune/fem-dg/models/defaultmodel.hh>
+
+
+/**********************************************
+ * Analytical model *
+ *********************************************/
+
+/**
+ * @brief Traits class for HeatEqnModel
+ */
+template <class GridPart >
+class HeatEqnModelTraits
+{
+public:
+ typedef GridPart GridPartType;
+ typedef typename GridPartType :: GridType GridType;
+ static const int dimDomain = GridType::dimensionworld;
+ static const int dimRange = DIMRANGE;
+ static const int dimGradRange = dimRange * dimDomain ;
+ // Definition of domain and range types
+ typedef Dune::FieldVector< double, dimDomain > DomainType;
+ typedef Dune::FieldVector< double, dimDomain-1 > FaceDomainType;
+ typedef Dune::FieldVector< double, dimRange > RangeType;
+ typedef Dune::FieldVector< double, dimGradRange > GradientType;
+ // ATTENTION: These are matrices (c.f. HeatEqnModel)
+ typedef Dune::FieldMatrix< double, dimRange, dimDomain > FluxRangeType;
+ typedef Dune::FieldMatrix< double, dimRange, dimDomain > JacobianRangeType;
+ typedef Dune::FieldMatrix< double, dimGradRange, dimDomain > DiffusionRangeType;
+ typedef Dune::FieldMatrix< double, dimDomain, dimDomain > DiffusionMatrixType;
+ typedef typename GridType :: template Codim< 0 > :: Entity EntityType;
+ typedef typename GridPartType :: IntersectionIteratorType IntersectionIterator;
+ typedef typename IntersectionIterator :: Intersection IntersectionType;
+
+ //typedef Dune::Fem::MinModLimiter< FieldType > LimiterFunctionType ;
+ //typedef SuperBeeLimiter< FieldType > LimiterFunctionType ;
+ //typedef VanLeerLimiter< FieldType > LimiterFunctionType ;
+};
+
+/**
+ * @brief describes the analytical model
+ *
+ * This is an description class for the problem
+ * \f{eqnarray*}{ V + \nabla a(U) & = & 0 \\
+ * \partial_t U + \nabla (F(U)+A(U,V)) & = & 0 \\
+ * U & = & g_D \\
+ * \nabla U \cdot n & = & g_N \f}
+ *
+ * where each class methods describes an analytical function.
+ * <ul>
+ * <li> \f$F\f$: advection() </li>
+ * <li> \f$a\f$: diffusion1() </li>
+ * <li> \f$A\f$: diffusion2() </li>
+ * <li> \f$g_D\f$ boundaryValue() </li>
+ * <li> \f$g_N\f$ boundaryFlux1(), boundaryFlux2() </li>
+ * </ul>
+ *
+ * \attention \f$F(U)\f$ and \f$A(U,V)\f$ are matrix valued, and therefore the
+ * divergence is defined as
+ *
+ * \f[ \Delta M = \nabla \cdot (\nabla \cdot (M_{i\cdot})^t)_{i\in 1\dots n} \f]
+ *
+ * for a matrix \f$M\in \mathbf{M}^{n\times m}\f$.
+ *
+ * @param GridPart GridPart for extraction of dimension
+ * @param ProblemType Class describing the initial(t=0) and exact solution
+ */
+
+
+////////////////////////////////////////////////////////
+//
+// Analytical model for the Heat Equation
+// dx(u) + div(uV) - epsilon*lap(u)) = 0
+// where V is constant vector
+//
+////////////////////////////////////////////////////////
+template <class GridPartType,class ProblemImp>
+class HeatEqnModel : public DefaultModel < HeatEqnModelTraits< GridPartType > >
+{
+public:
+ typedef ProblemImp ProblemType ;
+
+ static const int ConstantVelocity = ProblemType :: ConstantVelocity;
+ typedef typename GridPartType :: GridType GridType;
+ typedef HeatEqnModelTraits< GridPartType > Traits;
+ static const int dimDomain = Traits :: dimDomain ;
+ static const int dimRange = Traits :: dimRange ;
+ typedef typename Traits :: DomainType DomainType;
+ typedef typename Traits :: RangeType RangeType;
+ typedef typename Traits :: GradientType GradientType;
+ typedef typename Traits :: FluxRangeType FluxRangeType;
+ typedef typename Traits :: DiffusionRangeType DiffusionRangeType;
+ typedef typename Traits :: DiffusionMatrixType DiffusionMatrixType;
+ typedef typename Traits :: FaceDomainType FaceDomainType;
+ typedef typename Traits :: JacobianRangeType JacobianRangeType;
+
+ typedef typename Traits :: EntityType EntityType;
+ typedef typename Traits :: IntersectionType IntersectionType;
+
+ HeatEqnModel(const HeatEqnModel& otehr);
+ const HeatEqnModel &operator=(const HeatEqnModel &other);
+public:
+ /**
+ * @brief Constructor
+ *
+ * initializes model parameter
+ *
+ * @param problem Class describing the initial(t=0) and exact solution
+ */
+ HeatEqnModel(const ProblemType& problem) :
+ problem_(problem),
+ epsilon_(problem.epsilon()),
+ tstepEps_( getTStepEps() ),
+ velocity_( getVelocity() )
+ {}
+
+ inline bool hasFlux() const { return true ; }
+ inline const ProblemType& problem() const { return problem_; }
+ inline bool hasStiffSource() const { return problem_.hasStiffSource(); }
+ inline bool hasNonStiffSource() const { return problem_.hasNonStiffSource(); }
+
+ inline double nonStiffSource( const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ const GradientType& du,
+ RangeType & s) const
+ {
+ //FieldMatrixConverter< GradientType, JacobianRangeType> jac( du );
+ return nonStiffSource( en, time, x, u, s );
+ }
+
+ inline double nonStiffSource( const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ const JacobianRangeType& jac,
+ RangeType & s) const
+ {
+ return nonStiffSource( en, time, x, u, s );
+ }
+
+ inline double nonStiffSource( const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ RangeType & s) const
+ {
+ DomainType xgl = en.geometry().global( x );
+ return problem_.nonStiffSource( xgl, time, u, s );
+ }
+
+ inline double stiffSource( const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ const GradientType& du,
+ RangeType & s) const
+ {
+ //FieldMatrixConverter< GradientType, JacobianRangeType> jac( du );
+ return stiffSource( en, time, x, u, s );
+ }
+
+
+ inline double stiffSource( const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ const JacobianRangeType& jac,
+ RangeType & s) const
+ {
+ return stiffSource( en, time, x, u, s );
+ }
+
+ inline double stiffSource( const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ RangeType & s) const
+ {
+ DomainType xgl = en.geometry().global( x );
+ return problem_.stiffSource( xgl, time, u, s );
+ }
+
+ /**
+ * @brief advection term \f$F\f$
+ *
+ * @param en entity on which to evaluate the advection term
+ * @param time current time of TimeProvider
+ * @param x coordinate local to entity
+ * @param u \f$U\f$
+ * @param f \f$f(U)\f$
+ */
+ inline void advection(const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ FluxRangeType & f) const
+ {
+ // evaluate velocity V
+ DomainType v;
+ velocity( en, time, x, u, v );
+
+ //f = uV;
+ for( int r=0; r<dimRange; ++r )
+ for( int d=0; d<dimDomain; ++d )
+ f[r][d] = v[ d ] * u[ r ];
+ }
+
+ /**
+ * @brief velocity calculation, is called by advection()
+ */
+ inline void velocity(const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ DomainType& v) const
+ {
+ velocity( en.geometry().global(x), time, v );
+ }
+
+ /*
+ * @brief velocity calculation, is called by advection()
+ */
+ inline void velocity(const DomainType& xGlobal,
+ const double time,
+ DomainType& v) const
+ {
+ problem_.velocity(xGlobal, time, v);
+ }
+
+ /**
+ * @brief diffusion term \f$a\f$
+ */
+ inline void jacobian(const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ DiffusionRangeType& a) const
+ {
+ a = 0;
+
+ assert( a.rows == dimRange * dimDomain );
+ assert( a.cols == dimDomain );
+
+ for (int r=0;r<dimRange;r++)
+ for (int d=0;d<dimDomain;d++)
+ a[dimDomain*r+d][d] = u[r];
+ }
+
+ inline void eigenValues(const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ RangeType& maxValue) const
+ {
+ FluxRangeType A(0);
+ maxValue = RangeType( problem_.diffusion( u,A ) );
+ return;
+ abort();
+ /*
+ DiffusionMatrixType K(0) ;
+ DomainType values ;
+ for( int r = 0; r<dimRange; ++r )
+ {
+ // setup matrix
+ for(int i=0; i<dimDomain; ++i)
+ K[i][i] = problem_.diffusion( u,A );
+
+ // calculate eigenvalues
+ Dune::FMatrixHelp :: eigenValues( K, values );
+ // take max eigenvalue
+ maxValue[ r ] = values.infinity_norm();
+ }
+ */
+ }
+
+ /**
+ * @brief diffusion term \f$A\f$
+ */
+ template <class JacobianType>
+ inline void diffusion(const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ const JacobianType& jac,
+ FluxRangeType& A) const
+ {
+ // copy v to A
+ A = jac;
+
+ // apply diffusion coefficient
+ //double d = 0.; //(1.-en.geometry().global(x)[0])*en.geometry().global(x)[0]+
+ // //(1.-en.geometry().global(x)[1])*en.geometry().global(x)[1];
+ A *= problem_.diffusion( u, A );//*(1.+d);
+ }
+
+ inline void diffusion(const EntityType& en,
+ const double time,
+ const DomainType& x,
+ const RangeType& u,
+ const GradientType& vecJac,
+ FluxRangeType& A) const
+ {
+ Dune::Fem::FieldMatrixConverter< GradientType, FluxRangeType> jac( vecJac );
+ diffusion( en, time, x, u, jac, A );
+ }
+
+ inline double diffusionTimeStep( const IntersectionType &it,
+ const double enVolume,
+ const double circumEstimate,
+ const double time,
+ const FaceDomainType& x,
+ const RangeType& u ) const
+ {
+ return tstepEps_;
+ }
+
+ /** \brief convert conservative to primitive variables
+ * \note This method is used only for additional output to hdd
+ */
+ inline void conservativeToPrimitive( const double time,
+ const DomainType& xgl,
+ const RangeType& cons,
+ RangeType& prim,
+ const bool ) const
+ {
+ prim = cons;
+ }
+
+public:
+ /**
+ * @brief checks for existence of dirichlet boundary values
+ */
+ inline bool hasBoundaryValue(const IntersectionType& it,
+ const double time,
+ const FaceDomainType& x) const
+ {
+ return true;
+ }
+
+ /**
+ * @brief neuman boundary values \f$g_N\f$ for pass2
+ */
+ inline double boundaryFlux(const IntersectionType& it,
+ const double time,
+ const FaceDomainType& x,
+ const RangeType& uLeft,
+ const GradientType& vLeft,
+ RangeType& gLeft) const
+ {
+ gLeft = 0.;
+ return 0.;
+ }
+
+ /**
+ * @brief neuman boundary values \f$g_N\f$ for pass1
+ */
+ inline double boundaryFlux(const IntersectionType& it,
+ const double time,
+ const FaceDomainType& x,
+ const RangeType& uLeft,
+ RangeType& gLeft) const
+ {
+ gLeft = 0.;
+ return 0.;
+ }
+
+ inline double diffusionBoundaryFlux( const IntersectionType& it,
+ const double time,
+ const FaceDomainType& x,
+ const RangeType& uLeft,
+ const GradientType& gradLeft,
+ RangeType& gLeft ) const
+ {
+ Dune::Fem::FieldMatrixConverter< GradientType, JacobianRangeType> jacLeft( gradLeft );
+ return diffusionBoundaryFlux( it, time, x, uLeft, jacLeft, gLeft );
+ }
+
+ /** \brief boundary flux for the diffusion part
+ */
+ template <class JacobianRangeImp>
+ inline double diffusionBoundaryFlux( const IntersectionType& it,
+ const double time,
+ const FaceDomainType& x,
+ const RangeType& uLeft,
+ const JacobianRangeImp& jacLeft,
+ RangeType& gLeft ) const
+ {
+ std::cerr <<"diffusionBoundaryFlux shouldn't be used for this testcase" <<std::endl;
+ abort();
+ }
+
+ /**
+ * @brief dirichlet boundary values
+ */
+ inline void boundaryValue(const IntersectionType& it,
+ const double time,
+ const FaceDomainType& x,
+ const RangeType& uLeft,
+ RangeType& uRight) const
+ {
+#ifdef TESTOPERATOR
+ uRight = 0;
+ return;
+#endif
+ DomainType xgl = it.geometry().global( x );
+ problem_.evaluate(xgl, time, uRight);
+ }
+
+
+ // here x is in global coordinates
+ inline void maxSpeed( const EntityType& entity,
+ const double time,
+ const DomainType& xGlobal,
+ const DomainType& normal,
+ const RangeType& u,
+ double& advspeed,
+ double& totalspeed ) const
+ {
+ DomainType v;
+ problem_.velocity( xGlobal, time, v );
+ advspeed = v * normal ;
+ totalspeed = advspeed;
+ }
+
+ protected:
+ DomainType getVelocity() const
+ {
+ DomainType velocity(0);
+ if(ConstantVelocity) {
+ problem_.velocity(velocity,velocity);
+ }
+ return velocity;
+ }
+ double getTStepEps() const
+ {
+ // if diffusionTimeStep is set to non-zero in the parameterfile, the
+ // deltaT in the timeprovider is updated according to the diffusion
+ // parameter epsilon.
+ bool diff_tstep;
+ Dune::Fem::Parameter::get("femhowto.diffusionTimeStep", diff_tstep);
+ return diff_tstep ? epsilon_ : 0;
+ }
+
+ protected:
+ const ProblemType& problem_;
+ const double epsilon_;
+ const double tstepEps_;
+ public:
+ const DomainType velocity_;
+
+};
+
+
+/**
+ * @brief defines the advective flux
+ */
+template <class ModelType>
+class UpwindFlux {
+public:
+ typedef ModelType Model;
+ typedef typename Model::Traits Traits;
+ enum { dimRange = Model::dimRange };
+ typedef typename Model :: DomainType DomainType;
+ typedef typename Model :: RangeType RangeType;
+ typedef typename Model :: FluxRangeType FluxRangeType;
+ typedef typename Model :: DiffusionRangeType DiffusionRangeType;
+ typedef typename Model :: FaceDomainType FaceDomainType;
+ typedef typename Model :: EntityType EntityType;
+ typedef typename Model :: IntersectionType IntersectionType;
+protected:
+ template <class Model, bool constVelo>
+ struct Velocity
+ {
+ /**
+ * @brief computes and returns the wind direction
+ */
+ static inline double upwind(const Model& model,
+ const IntersectionType& it,
+ const double time,
+ const FaceDomainType& x,
+ const RangeType& uLeft)
+ {
+ const DomainType normal = it.integrationOuterNormal(x);
+ DomainType velocity;
+ model.velocity(*it.inside(),time,
+ it.geometryInInside().global(x),
+ uLeft,velocity);
+ return normal*velocity;
+ }
+ };
+
+ template <class Model>
+ struct Velocity<Model,true>
+ {
+ /**
+ * @brief computes and returns the wind direction for models with
+ * constant velocity
+ */
+ static inline double upwind( const Model& model,
+ const IntersectionType& it,
+ const double time,
+ const FaceDomainType& x,
+ const RangeType& uLeft )
+ {
+ const DomainType normal = it.integrationOuterNormal(x);
+ return normal * model.velocity_;
+ }
+ };
+
+public:
+ /**
+ * @brief Constructor
+ */
+ UpwindFlux(const Model& mod) : model_(mod) {}
+
+ static std::string name () { return "UpwindFlux"; }
+
+ const Model& model() const {return model_;}
+
+ /**
+ * @brief evaluates the flux \f$g(u,v)\f$
+ *
+ * @return maximum wavespeed * normal
+ */
+ template <class QuadratureImp>
+ inline double numericalFlux( const IntersectionType& it,
+ const EntityType& inside,
+ const EntityType& outside,
+ const double time,
+ const QuadratureImp& faceQuadInner,
+ const QuadratureImp& faceQuadOuter,
+ const int quadPoint,
+ const RangeType& uLeft,
+ const RangeType& uRight,
+ RangeType& gLeft,
+ RangeType& gRight ) const
+ {
+ const FaceDomainType& x = faceQuadInner.localPoint( quadPoint );
+ const double upwind = Velocity<Model,Model::ConstantVelocity>::
+ upwind(model_,it,time,x,uLeft);
+
+ if (upwind>0)
+ gLeft = uLeft;
+ else
+ gLeft = uRight;
+ gLeft *= upwind;
+ gRight = gLeft;
+ return std::abs(upwind);
+ }
+protected:
+ const Model& model_;
+};
+
+#endif
diff --git a/dune/fem-dg/test/advdiff/parameter b/dune/fem-dg/test/advdiff/parameter
new file mode 100644
index 00000000..ce6d324d
--- /dev/null
+++ b/dune/fem-dg/test/advdiff/parameter
@@ -0,0 +1,56 @@
+# DATA WRITER
+#------------
+
+fem.prefix: ./data # specify directory for data output (is created if not exists)
+fem.io.datafileprefix: heat # prefix data data files
+fem.io.savestep: 0.1 # save data every time interval
+fem.io.savecount: -1 # save every i-th step
+
+
+# GRID SOLUTION
+#--------------
+
+gridsol.savestep: 0.01
+gridsol.firstwrite: 0.1
+gridsol.filename: heat-checkpoint
+
+
+# GENERAL
+#--------
+
+paramfile: ../parameter/paramBase
+
+
+# PROBLEM SETUP
+#--------------
+
+# problem: heat, quasi, plaplace
+femhowto.problem: deformational
+
+femhowto.endTime: 1.0
+femhowto.diffusionTimeStep: 1
+femhowto.epsilon: 0.0
+femhowto.plaplace: 3.0
+femhowto.xvelocity: 1. # the only advection part for the linear heat eqn
+femhowto.yvelocity: 3. # the only advection part for the linear heat eqn
+femhowto.zvelocity: 0. # the only advection part for the linear heat eqn
+
+
+# DOMAIN
+#-------
+
+fem.io.macroGridFile_1d: ../macrogrids/unitcube1.dgf
+fem.io.macroGridFile_2d: ../macrogrids/unitcube2_per.dgf
+#fem.io.macroGridFile_2d: ../macrogrids/unitcube2_unstr.dgf
+#fem.io.macroGridFile_2d: ../macrogrids/test2.dgf
+#fem.io.macroGridFile_2d: ../macrogrids/unitgrid2d_unstr_2.dgf
+#fem.io.macroGridFile_2d: ../macrogrids/stability.dgf
+#fem.io.macroGridFile_2d: ../macrogrids/cdgpaper_test.dgf
+fem.io.macroGridFile_3d: ../macrogrids/unitcube3.dgf
+
+
+# SOLVER
+#-------
+dgdiffusionflux.upwind: -1 -0.001
+paramfile: ../parameter/paramSolver
+femhowto.maxTimeStep: 0.5
diff --git a/dune/fem-dg/test/advdiff/problem.hh b/dune/fem-dg/test/advdiff/problem.hh
new file mode 100644
index 00000000..893607ee
--- /dev/null
+++ b/dune/fem-dg/test/advdiff/problem.hh
@@ -0,0 +1,255 @@
+#ifndef DUNE_PROBLEM_HH__
+#define DUNE_PROBLEM_HH__
+
+// dune-fem includes
+#include <dune/fem/io/parameter.hh>
+#include <dune/fem/space/common/functionspace.hh>
+
+// local includes
+#include <dune/fem-dg/models/defaultprobleminterfaces.hh>
+
+namespace Dune {
+
+/**
+ * @brief describes the initial and exact solution of the advection-diffusion model
+ * for given constant velocity vector v=(v1,v2)
+ *
+ * \f[u(x,y,z,t):=\displaystyle{\sum_{i=0}^{1}} T_i(t) \cdot X_i(x) \cdot
+ * Y_i(y) \cdot Z_i(z)\f]
+ *
+ * with
+ *
+ * \f{eqnarray*}{
+ * T_0(t) &:=& e^{-\varepsilon t \pi^2 (2^2 + 1^2 + 1.3^2 )} \\
+ * X_0(x) &:=& 0.6\cdot \cos(2\pi (x-v_1t)) + 0.8\cdot \sin(2\pi (x-v_1t)) \\
+ * Y_0(y) &:=& 1.2\cdot \cos(1\pi (y-v_2t)) + 0.4\cdot \sin(1\pi (y-v_2t)) \\
+ * Z_0(z) &:=& 0.1\cdot \cos(1.3\pi (z-v_3t)) - 0.4\cdot \sin(1.3\pi (z-v_3t)) \\
+ * T_1(t) &:=& e^{-\varepsilon t \pi^2 (0.7^2 + 0.5^2 + 0.1^2 )} \\
+ * X_1(x) &:=& 0.9\cdot \cos(0.7\pi (x-v_1t)) + 0.2\cdot \sin(0.7\pi (x-v_1t)) \\
+ * Y_1(y) &:=& 0.3\cdot \cos(0.5\pi (y-v_2t)) + 0.1\cdot \sin(0.5\pi (y-v_2t))
+ * Z_1(z) &:=& -0.3\cdot \cos(0.1\pi (z-v_3t)) + 0.2\cdot \sin(0.1\pi (z-v_3t)) \\
+ * \f}
+ *
+ * This is a solution of the AdvectionDiffusionModel for \f$g_D = u|_{\partial
+ * \Omega}\f$.
+ *
+ */
+template <class GridType> /*@LST0S@*/
+struct U0 : public EvolutionProblemInterface<
+ Fem::FunctionSpace< double, double, GridType::dimension, DIMRANGE>,
+ true >
+{ /*@LST0E@*/
+public:
+ typedef EvolutionProblemInterface<
+ Fem::FunctionSpace< double, double,
+ GridType::dimension, DIMRANGE >,
+ true > BaseType;
+
+ enum{ dimDomain = BaseType :: dimDomain };
+ typedef typename BaseType :: DomainType DomainType;
+ typedef typename BaseType :: RangeType RangeType;
+ typedef typename BaseType :: JacobianRangeType JacobianRangeType;
+
+ typedef Fem::Parameter ParameterType;
+
+ /**
+ * @brief define problem parameters
+ */
+ U0 () : /*@LST0S@*/
+ BaseType () ,
+ velocity_( 0 ),
+ startTime_( ParameterType::getValue<double>("femhowto.startTime",0.0) ),
+ epsilon_( ParameterType::getValue<double>("femhowto.epsilon",0.1) )
+ { /*@LST0E@*/
+ std::cout <<"Problem: HeatEqnWithAdvection, epsilon " << epsilon_ << "\n";
+ //std::cout <<"Problem: HeatEqnWithAdvection, epsilon_" << epsilon_ << "\n";
+
+ // an advection vector, default to 0
+ velocity_[0] = ParameterType::getValue<double>( "femhowto.xvelocity", 0. );
+ if (dimDomain > 1)
+ velocity_[1] = ParameterType::getValue<double>( "femhowto.yvelocity", 0. );
+ if (dimDomain > 2)
+ velocity_[2] = ParameterType::getValue<double>( "femhowto.zvelocity", 0. );
+
+ max_n_of_coefs_ = 2;
+
+ //x coordinate
+ common_coef_x_[0] = 2.0;
+ sin_coef_x_[0] = 0.8;
+ cos_coef_x_[0] = 0.6;
+
+ //x coordinate
+ common_coef_x_[1] = 0.7;
+ sin_coef_x_[1] = 0.2;
+ cos_coef_x_[1] = 0.9;
+
+ //y coordinate
+ common_coef_y_[0] = 1.0;
+ sin_coef_y_[0] = 0.4;
+ cos_coef_y_[0] = 1.2;
+
+
+ //y coordinate
+ common_coef_y_[1] = 0.5;
+ sin_coef_y_[1] = 0.1;
+ cos_coef_y_[1] = 0.3;
+
+
+ //z coordinate
+ common_coef_z_[0] = 1.3;
+ sin_coef_z_[0] = -0.4;
+ cos_coef_z_[0] = 0.1;
+
+ //z coordinate
+ common_coef_z_[1] = 0.1;
+ sin_coef_z_[1] = 0.2;
+ cos_coef_z_[1] = -0.3;
+
+ myName = "LinHeatEqn";
+ }
+
+ //! this problem has no source term
+ bool hasStiffSource() const { return false; }
+ bool hasNonStiffSource() const { return false; }
+
+ double stiffSource(const DomainType& arg,
+ const double t,
+ const RangeType& u,
+ RangeType& res) const
+ {
+ return 0.0;
+ }
+
+ double nonStiffSource(const DomainType& arg,
+ const double t,
+ const RangeType& u,
+ RangeType& res) const
+ {
+ return 0.0;
+ }
+
+ double diffusion( const RangeType& u, const JacobianRangeType& gradU ) const
+ {
+ return epsilon();
+ }
+
+ //! return start time
+ double startTime() const { return startTime_; }
+
+ //! return start time
+ double epsilon() const { return epsilon_; }
+
+ /**
+ * @brief getter for the velocity
+ */
+ void velocity(const DomainType& x, DomainType& v) const
+ {
+ v = velocity_;
+ }
+
+ /**
+ * @brief evaluates \f$ u_0(x) \f$
+ */
+ void evaluate(const DomainType& arg, RangeType& res) const /*@LST0S@@LST0E@*/
+ {
+ evaluate(arg, startTime_, res);
+ }
+
+ template <class T>
+ double SQR( const T& a ) const
+ {
+ return ( a * a );
+ }
+
+
+ /**
+ * @brief evaluate exact solution
+ */
+ void evaluate(const DomainType& arg, const double t, RangeType& res) const /*@LST0S@@LST0E@*/
+ {
+
+ res = 0;
+ DomainType x(arg);
+ x -= DomainType(0.5);
+ x.axpy(-t,velocity_);
+
+ for(int i=0;i<max_n_of_coefs_;++i)
+ {
+ if(dimDomain == 1)
+ res += exp(-epsilon_*t*(SQR(common_coef_x_[i]*M_PI)))\
+ *((cos_coef_x_[i]*cos(common_coef_x_[i]*M_PI*x[0])\
+ + sin_coef_x_[i]*sin(common_coef_x_[i]*M_PI*x[0])));
+ else if(dimDomain == 2)
+ res += exp(-epsilon_*t*(SQR(common_coef_x_[i]*M_PI)+\
+ SQR(common_coef_y_[i]*M_PI)))\
+ *((cos_coef_x_[i]*cos(common_coef_x_[i]*M_PI*x[0])\
+ + sin_coef_x_[i]*sin(common_coef_x_[i]*M_PI*x[0]))\
+ *(cos_coef_y_[i]*cos(common_coef_y_[i]*M_PI*x[1])\
+ + sin_coef_y_[i]*sin(common_coef_y_[i]*M_PI*x[1])));
+ else if(dimDomain == 3)
+ res += exp(-epsilon_*t*(SQR(common_coef_x_[i]*M_PI)+\
+ SQR(common_coef_y_[i]*M_PI)+\
+ SQR(common_coef_z_[i]*M_PI)))\
+ *((cos_coef_x_[i]*cos(common_coef_x_[i]*M_PI*x[0])\
+ + sin_coef_x_[i]*sin(common_coef_x_[i]*M_PI*x[0]))\
+ *(cos_coef_y_[i]*cos(common_coef_y_[i]*M_PI*x[1])\
+ + sin_coef_y_[i]*sin(common_coef_y_[i]*M_PI*x[1]))\
+ *(cos_coef_z_[i]*cos(common_coef_z_[i]*M_PI*x[2])\
+ + sin_coef_z_[i]*sin(common_coef_z_[i]*M_PI*x[2])));
+ }
+ }
+
+ /**
+ * @brief latex output for EocOutput
+ */
+ std::string description() const
+ {
+ std::ostringstream ofs;
+
+ ofs << "Problem: " << myName
+ << ", Epsilon: " << epsilon_
+ << ", Advection vector: (" <<velocity_[0];
+
+ if (dimDomain > 1)
+ ofs <<"," <<velocity_[1];
+ if (dimDomain > 2)
+ ofs <<"," <<velocity_[2];
+ ofs <<")";
+
+ ofs << ", End time: " << ParameterType::getValue<double>("femhowto.endTime");
+
+ return ofs.str();
+ }
+
+ /* \brief finalize the simulation using the calculated numerical
+ * solution u for this problem
+ *
+ * \param[in] variablesToOutput Numerical solution in the suitably chosen variables
+ * \param[in] eocloop Specific EOC loop
+ */
+ template< class DiscreteFunctionType >
+ void finalizeSimulation( DiscreteFunctionType& variablesToOutput,
+ const int eocloop) const
+ {}
+
+private:
+ DomainType velocity_;
+ int max_n_of_coefs_;
+ double common_coef_x_[2];
+ double sin_coef_x_[2];
+ double cos_coef_x_[2];
+ double common_coef_y_[2];
+ double sin_coef_y_[2];
+ double cos_coef_y_[2];
+ double common_coef_z_[2];
+ double sin_coef_z_[2];
+ double cos_coef_z_[2];
+ const double startTime_;
+public:
+ const double epsilon_;
+ std::string myName;
+};
+
+}
+#endif /*DUNE_PROBLEM_HH__*/
+
diff --git a/dune/fem-dg/test/advdiff/problemQuasiHeatEqn.hh b/dune/fem-dg/test/advdiff/problemQuasiHeatEqn.hh
new file mode 100644
index 00000000..e97fe39a
--- /dev/null
+++ b/dune/fem-dg/test/advdiff/problemQuasiHeatEqn.hh
@@ -0,0 +1,166 @@
+#ifndef DUNE_PROBLEM__QUASI_HH__
+#define DUNE_PROBLEM__QUASI_HH__
+
+#include <dune/common/version.hh>
+
+// dune-fem includes
+#include <dune/fem/misc/linesegmentsampler.hh>
+#include <dune/fem/io/parameter.hh>
+#include <dune/fem/space/common/functionspace.hh>
+
+// local includes
+#include <dune/fem-dg/models/defaultprobleminterfaces.hh>
+
+
+namespace Dune {
+
+/**
+ * @brief describes the initial and exact solution of the advection-diffusion model
+ */
+template <class GridType> /*@LST0S@*/
+struct QuasiHeatEqnSolution : public EvolutionProblemInterface<
+ Dune::Fem::FunctionSpace< double, double, GridType::dimension, DIMRANGE>,
+ true >
+{ /*@LST0E@*/
+public:
+ typedef EvolutionProblemInterface<
+ Dune::Fem::FunctionSpace< double, double,
+ GridType::dimension, DIMRANGE >,
+ true > BaseType;
+
+ enum{ dimDomain = BaseType :: dimDomain };
+ typedef typename BaseType :: DomainType DomainType;
+ typedef typename BaseType :: RangeType RangeType;
+ typedef typename BaseType :: JacobianRangeType JacobianRangeType;
+
+ typedef Fem :: Parameter ParameterType ;
+
+ /**
+ * @brief define problem parameters
+ */
+ QuasiHeatEqnSolution () : /*@LST0S@*/
+ BaseType () ,
+ velocity_( 0 ),
+ startTime_( ParameterType::getValue<double>("femhowto.startTime",0.0) ),
+ epsilon_( ParameterType::getValue<double>("femhowto.epsilon") )
+ { /*@LST0E@*/
+ if ( (DIMRANGE != 1) || (dimDomain != 2) )
+ {
+ std::cout <<"QuasiHeatEqn only supports dimRange=1 and dimDomain=2\n";
+ abort();
+ }
+ std::cout <<"Problem: QuasiHeatEqn" <<epsilon_ << " epsilon \n";
+
+ myName = "QuasiHeatEqn";
+ }
+
+ double diffusion( const RangeType& u, const JacobianRangeType& gradU ) const
+ {
+ return epsilon_*u;
+ }
+
+ double startTime() const { return startTime_; }
+
+ double epsilon() const { return epsilon_; }
+
+ /**
+ * @brief getter for the velocity
+ */
+ void velocity(const DomainType& x, DomainType& v) const
+ {
+ v = velocity_;
+ }
+
+ /**
+ * @brief evaluates \f$ u_0(x) \f$
+ */
+ void evaluate(const DomainType& arg, RangeType& res) const /*@LST0S@@LST0E@*/
+ {
+ evaluate(arg, startTime_, res);
+ }
+
+ /**
+ * @brief evaluate exact solution
+ */
+ void evaluate(const DomainType& arg, const double t, RangeType& res) const /*@LST0S@@LST0E@*/
+ {
+ res = std::sin(arg[0]) * std::sin(arg[1]) * std::exp( -epsilon_*t );
+ }
+
+ bool hasStiffSource() const { return true; }
+ bool hasNonStiffSource() const { return false; }
+
+ /**
+ * @brief evaluate stiff source function
+ */
+ double nonStiffSource( const DomainType& arg,
+ const double t,
+ const RangeType& u,
+ RangeType& res ) const /*@LST0S@@LST0E@*/
+ {
+ res = 0.;
+
+ // time step restriction from the stiff source term
+ return 0.;
+ }
+
+ /**
+ * @brief evaluate non stiff source function
+ */
+ double stiffSource( const DomainType& arg,
+ const double t,
+ const RangeType& u,
+ RangeType& res ) const /*@LST0S@@LST0E@*/
+ {
+ const double x = arg[0];
+ const double y = arg[1];
+ const double cosx2 = cos(x)*cos(x);
+ const double sinx2 = sin(x)*sin(x);
+ const double cosy2 = cos(y)*cos(y);
+ const double siny2 = sin(y)*sin(y);
+
+ res = (-std::exp(-2.*epsilon_*t )*(cosx2*siny2 + cosy2*sinx2) + 2*u*u - u);
+ res *= epsilon_;
+
+ // time step restriction from the non stiff source term
+ return 0.0;
+ }
+
+ /**
+ * @brief latex output for EocOutput
+ */
+ std::string description() const
+ {
+ std::ostringstream ofs;
+
+ ofs << "Problem: " << myName
+ << ", Epsilon: " << epsilon_
+ << ", End time: " << ParameterType::getValue<double>("femhowto.endTime");
+
+ return ofs.str();
+ }
+
+
+ /* \brief finalize the simulation using the calculated numerical
+ * solution u for this problem
+ *
+ * \param[in] variablesToOutput Numerical solution in the suitably chosen variables
+ * \param[in] eocloop Specific EOC loop
+ */
+ template< class DiscreteFunctionType >
+ void finalizeSimulation( DiscreteFunctionType& variablesToOutput,
+ const int eocloop) const
+ {}
+
+private:
+ DomainType velocity_;
+ double startTime_;
+
+public:
+ double epsilon_;
+ std::string myName;
+};
+
+}
+#endif /*DUNE_PROBLEM_HH__*/
+
diff --git a/dune/fem-dg/test/advdiff/problemcreator.hh b/dune/fem-dg/test/advdiff/problemcreator.hh
new file mode 100644
index 00000000..d6095fa3
--- /dev/null
+++ b/dune/fem-dg/test/advdiff/problemcreator.hh
@@ -0,0 +1,99 @@
+#ifndef FEMHOWTO_HEATSTEPPER_HH
+#define FEMHOWTO_HEATSTEPPER_HH
+#include <config.h>
+
+// dune-fem includes
+#include <dune/fem/io/parameter.hh>
+
+#include <dune/grid/io/file/dgfparser/dgfparser.hh>
+#include <dune/fem-dg/operator/fluxes/eulerfluxes.hh>
+
+// local includes
+#include <dune/fem-dg/operator/fluxes/diffusionflux.hh>
+#include "problem.hh"
+#include "problemQuasiHeatEqn.hh"
+#include "deformationalflow.hh"
+#include "models.hh"
+
+
+template< class GridType >
+struct ProblemGenerator
+{
+ typedef Dune :: EvolutionProblemInterface<
+ Dune::Fem::FunctionSpace< double, double, GridType::dimension,
+ DIMRANGE>,
+ false > ProblemType;
+ //typedef HeatProblemType ProblemType;
+
+ template< class GridPart >
+ struct Traits
+ {
+ typedef ProblemType InitialDataType;
+ typedef HeatEqnModel< GridPart, InitialDataType > ModelType;
+ typedef LLFFlux< ModelType > FluxType;
+ //typedef UpwindFlux< ModelType > FluxType;
+ // choice of diffusion flux (see diffusionflux.hh for methods)
+ static const Dune :: DGDiffusionFluxIdentifier PrimalDiffusionFluxId
+ = Dune :: method_general ;
+ };
+
+
+ static inline std::string advectionFluxName()
+ {
+#if (FLUX==1)
+ return "LLF";
+#elif (FLUX==2)
+ return "HLL(Dennis)";
+#elif (FLUX==3)
+ return "HLLC(Dennis)";
+#elif (FLUX==4)
+ return "HLL2C";
+#elif (FLUX==5)
+ return "HLL2";
+#elif (FLUX==6)
+ return "HLLEM(Mhd)";
+#endif
+ }
+
+
+ static inline std::string diffusionFluxName()
+ {
+#ifdef EULER
+ return "";
+#elif (defined PRIMALDG)
+ return Dune::Fem::Parameter::getValue< std::string >("dgdiffusionflux.method");
+#else
+ return "LDG";
+#endif
+ }
+
+
+ static inline Dune::GridPtr<GridType>
+ initializeGrid( const std::string description )
+ {
+ // use default implementation
+ return initialize< GridType > ( description );
+ }
+
+ static ProblemType* problem()
+ {
+ // choice of explicit or implicit ode solver
+ //static const std::string probString[] = { "heat" ,"quasi", "deformational" };
+ //int probNr = Dune::Fem::Parameter::getEnum( "femhowto.problem", probString, 0 );
+ //if( probNr == 0 )
+ // return new Dune :: U0< GridType > ();
+ //else if ( probNr == 1 )
+ // return new Dune :: QuasiHeatEqnSolution< GridType > ();
+ //else if ( probNr == 2 )
+ return new Dune :: DeformationalFlow< GridType > ();
+ //else
+ //{
+ // abort();
+ // return 0;
+ //}
+ }
+};
+
+//#include <dune/fem-dg/stepper/advectionstepper.hh>
+
+#endif // FEMHOWTO_HEATSTEPPER_HH
--
GitLab