copy poisson implementation from dglagrange.

Not working yet.

dune/fem-dg/test/poisson/Makefile.am

+# install these headers 
+poissondir=$(includedir)/test/advdiff
+poisson_HEADERS = models.hh  problemcreator.hh steppertraits.hh
+
+LDADD = $(ALL_PKG_LDFLAGS) $(ALL_PKG_LIBS) $(LOCAL_LIBS) $(DUNEMPILDFLAGS) $(DUNEMPILIBS)
+
+BASEDIR=../../main
+
+#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
+EXTRA_PROGRAMS = advdiff pulse # quasiadvdiff quasiadvdiff12 plaplace plaplace12
+check_PROGRAMS = advdiffonep pulseonep # quasiadvdiff quasiadvdiff12 plaplace plaplace12
+
+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
+
+pulse_SOURCES = $(BASEDIR)/main.cc $(BASEDIR)/main_1.cc $(BASEDIR)/main_2.cc \
+                $(BASEDIR)/main_0.cc $(BASEDIR)/main_3.cc
+pulseonep_SOURCES = $(BASEDIR)/main.cc $(BASEDIR)/main_pol.cc
+
+AM_CPPFLAGS = $(USE_OMP) -I../../problems/advdiff  $(ALL_PKG_CPPFLAGS) -DGRIDDIM=$(GRIDDIM) \
+              -D$(GRIDTYPE) $(DUNEMPICPPFLAGS) \
+              -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)
+
+pulse_CPPFLAGS = $(AM_CPPFLAGS)
+pulseonep_CPPFLAGS = $(advdiff_CPPFLAGS) -DONLY_ONE_P -DPOLORDER=$(POLORDER)
+
+DISTCHECK_CONFIGURE_FLAGS = DOXYGEN="true"
+
+EXTRA_DIST = parameter
+
+CLEANFILES=manager.*.log eoc_*.tex
+
+PROG=advdiffonep
+# 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

dune/fem-dg/test/poisson/benchmarkproblem.hh

+#ifndef  DUNE_BENCHMARK_PROBLEM_HH
+#define  DUNE_BENCHMARK_PROBLEM_HH
+ 
+#include <dune/fem/io/parameter.hh>
+#include <dune/fem/space/common/functionspace.hh>
+
+// local includes
+#include "../common/probleminterfaces.hh"
+#include "benchmark.hh"
+
+
+namespace Dune {
+
+template <class GridType>                                          /*@LST0S@*/
+struct BenchmarkProblems : public ProblemInterface<
+                    Dune :: Fem :: FunctionSpace< double, double, GridType::dimension, DIMRANGE> >
+{                                                                  /*@LST0E@*/
+public:
+  typedef ProblemInterface<
+                 Dune :: Fem::FunctionSpace< double, double,
+                                        GridType::dimension, DIMRANGE >
+                          >                    BaseType;
+
+  enum{ dimDomain = BaseType :: dimDomain };
+  enum{ dim = dimDomain };
+  typedef typename BaseType :: DomainType                DomainType;
+  typedef typename BaseType :: RangeType                 RangeType;
+  typedef typename BaseType :: JacobianRangeType         JacobianRangeType;
+  typedef typename BaseType :: DiffusionMatrixType   DiffusionMatrixType;
+  typedef typename GridType :: ctype   FieldType ;
+
+  typedef DataFunctionIF< dimDomain, FieldType, FieldType > DataFunctionType;
+
+  /**
+   * @brief define problem parameters
+   */
+  BenchmarkProblems (const int problemNumber) :   
+    BaseType (),
+    data_(0)
+  {            
+    FieldType shift = Dune :: Fem::Parameter :: getValue< double > ("femhowto.globalshift", 0);
+    FieldType factor = Dune :: Fem::Parameter :: getValue< double > ("femhowto.factor", 1);
+    if( problemNumber == 0 )
+    {
+      data_ = new BenchMark_1<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 1 )
+    {
+      data_ = new BenchMark_1_2<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 2 )
+    {
+      data_ = new BenchMark_2<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 3 )
+    {
+      data_ = new BenchMark_3<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 4 )
+    {
+      data_ = new BenchMark_4<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 5 )
+    {
+      data_ = new BenchMark_5<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 6 )
+    {
+      data_ = new BenchMark_6<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 7 )
+    {
+      data_ = new BenchMark_7<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 8 )
+    {
+      //DUNE_THROW(InvalidStateException,"Problem 8 not available");
+      data_ = new CurvedRidges< dim, FieldType,FieldType> (shift,factor); 
+    }
+    if( problemNumber == 9 )
+    {
+      data_ = new BenchMark_9<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 10 )
+    {
+      data_ = new SinSin<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 11 )
+    {
+      data_ = new Excercise2_3< dim, FieldType,FieldType> (shift,factor); 
+    }
+    if( problemNumber == 12 )
+    {
+      data_ = new CastilloProblem<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 13 )
+    {
+      data_ = new InSpringingCorner<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 14 )
+    {
+      data_ = new RiviereProblem<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 15 )
+    {
+      data_ = new HeatProblem<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 16 )
+    {
+      data_ = new AlbertaProblem<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 17 )
+    {
+      data_ = new SinSin<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 18 )
+    {
+      data_ = new CosCos<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 19 )
+    {
+      data_ = new SinSinSin<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 20 )
+    {
+      data_ = new Hump<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 21 )
+    {
+      data_ = new BenchMark3d_1<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 23 )
+    {
+      data_ = new BenchMark3d_3<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 24 )
+    {
+      data_ = new BenchMark3d_4<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 25 )
+    {
+      data_ = new BenchMark3d_5<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 31 )
+    {
+      data_ = new BoundaryLayerProblem<dim,FieldType,FieldType> (shift,factor);
+    }
+    if( problemNumber == 32 )
+    {
+      data_ = new FicheraCorner<dim,FieldType,FieldType> (shift,factor);
+    }
+
+    if( data_ == 0 ) 
+    {
+      std::cerr << "ERROR: wrong problem number " << std::endl; 
+      abort();
+    }
+
+    myName = "Poisson Eqn.";
+  }
+
+  const DataFunctionType& data() const { assert( data_ ); return *data_; }
+
+  //! this problem has no source term
+  bool hasStiffSource() const { return true; }
+  bool hasNonStiffSource() const { return false; }
+
+  void f(const DomainType& arg,
+         RangeType& res) const
+  {
+    res = data().rhs( &arg[ 0 ] );
+  }
+
+  double stiffSource(const DomainType& arg,
+                     const double t,
+                     const RangeType& u,
+                     RangeType& res) const
+  {
+    res = data().rhs( &arg[ 0 ] );
+    return 0;
+  }
+
+  double nonStiffSource(const DomainType& arg,
+                const double t,
+                const RangeType& u,
+                RangeType& res) const
+  {
+    abort();
+    return 0.0;
+  }
+
+  //! return start time
+  double startTime() const { return 0; }
+
+  bool constantK() const { return data().constantLocalK();  }
+
+  void K(const DomainType& x, DiffusionMatrixType& m) const 
+  {
+    double k[ dimDomain ][ dimDomain ]; 
+    data().K( &x[0], k );
+
+    for(int i=0; i<dimDomain; ++i)
+      for(int j=0; j<dimDomain; ++j)
+         m[i][j] = k[i][j];
+  }
+
+  /**
+   * @brief getter for the velocity
+   */
+  void velocity(const DomainType& x, DomainType& v) const
+  {
+    double vv[ dimDomain ];
+    data().velocity( &x[0], vv );
+    for (int i=0;i<dimDomain; ++i)
+      v[i] = vv[i];
+  }
+
+  void u(const DomainType& arg, RangeType& res) const         /*@LST0S@@LST0E@*/
+  {
+    evaluate(arg, res );
+  }
+
+  /**
+   * @brief evaluates \f$ u_0(x) \f$
+   */
+  void evaluate(const DomainType& arg, RangeType& res) const         /*@LST0S@@LST0E@*/
+  {
+    evaluate(arg, 0, res);
+  }
+
+  /**
+   * @brief evaluate exact solution
+   */
+  void evaluate(const DomainType& arg, const double t, RangeType& res) const /*@LST0S@@LST0E@*/
+  {
+    res = data().exact( &arg[ 0 ] );
+  }
+
+  void gradient(const DomainType& x,
+                JacobianRangeType& grad) const 
+  { 
+    for (int i=0;i<RangeType::dimension;++i)
+      data().gradExact( &x[ 0 ], &grad[ i ][ 0 ] );
+  }
+
+  /**
+   * @brief latex output for EocOutput
+   */
+  std::string description() const
+  {
+    std::ostringstream ofs;
+
+    ofs << "Problem: " << myName ;
+    ofs << ", End time: " << Dune:: Fem ::  Parameter::getValue<double>("femhowto.endTime");
+
+    return ofs.str();
+  }
+
+protected:
+  DataFunctionType* data_;
+  std::string myName;
+};
+
+}
+#endif  /*DUNE_PROBLEM_HH__*/
+

dune/fem-dg/test/poisson/models.hh

+#ifndef POIS_MODELS_HH
+#define POIS_MODELS_HH
+
+#include <dune/fem/version.hh>
+#include <dune/fem/misc/fmatrixconverter.hh>
+#include <dune/fem/io/parameter.hh>
+
+#if HAVE_UGGRID
+#include <dune/grid/uggrid.hh>
+#endif
+
+namespace Dune
+{
+template <class ModelType>
+class UpwindFlux;
+}
+
+/**********************************************
+ * Analytical model                           *
+ *********************************************/
+/**
+ * @brief Traits class for PoissonModel
+ */
+template <class GridPart,int dimRange2,
+          int dimRange1=dimRange2* GridPart::GridType::dimensionworld>
+class PoissonModelTraits
+{
+public:
+  typedef GridPart                                                   GridPartType;
+  typedef typename GridPartType :: GridType                          GridType;
+  static const int dimDomain = GridType::dimensionworld;
+  static const int dimRange = dimRange2;
+  static const int dimGradRange = dimRange1;
+  // 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. PoissonModel)
+  typedef Dune::FieldMatrix< double, dimRange, dimDomain >           FluxRangeType;
+  typedef Dune::FieldMatrix< double, dimRange, dimDomain >           JacobianRangeType;
+  typedef Dune::FieldMatrix< double, dimGradRange, dimDomain >       DiffusionRangeType;
+  typedef typename GridType :: template Codim< 0 > :: Entity         EntityType;
+  typedef typename GridPartType :: IntersectionIteratorType          IntersectionIterator;
+  typedef typename IntersectionIterator :: Intersection              IntersectionType;
+
+};
+
+/**
+ * @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 PoissonModel
+{
+public:
+  typedef ProblemImp  ProblemType ;
+
+  typedef typename GridPartType :: GridType                          GridType;
+  static const int dimDomain = GridType::dimensionworld;
+  static const int dimRange = DIMRANGE;
+  typedef PoissonModelTraits< GridPartType, dimRange,
+                              dimRange*dimDomain >                   Traits;
+  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 :: FaceDomainType                          FaceDomainType;
+  typedef typename Traits :: JacobianRangeType                       JacobianRangeType;
+
+  typedef typename ProblemType :: DiffusionMatrixType   DiffusionMatrixType ;
+
+  typedef typename Traits :: EntityType                       EntityType;
+  typedef typename Traits :: IntersectionType                 IntersectionType;
+
+public:
+  static const int ConstantVelocity = false;
+  /**
+   * @brief Constructor
+   *
+   * initializes model parameter
+   *
+   * @param problem Class describing the initial(t=0) and exact solution
+   */
+  PoissonModel(const ProblemType& problem) : problem_(problem)
+  {
+  }
+
+  inline bool hasFlux() const { return true ; }
+
+  inline bool hasStiffSource() const { return true ; }
+  inline bool hasNonStiffSource() const { return false ; }
+
+  inline double nonStiffSource( const EntityType& en,
+                        const double time,
+                        const DomainType& x,
+                        const RangeType& u,
+                        const GradientType& du,
+                        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,
+                        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
+  {
+    abort();
+    DomainType xgl = en.geometry().global( x );
+    problem_.f( xgl, s );
+    return 0;
+  }
+
+  inline double stiffSource( const EntityType& en,
+                        const double time,
+                        const DomainType& x,
+                        const RangeType& u,
+                        const GradientType& du,
+                        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,
+                        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 );
+    problem_.f( xgl, s );
+    return 0;
+  }
+
+  /**
+   * @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, v );
+
+    //f = uV;
+    for( int r=0; r<dimRange; ++r )
+      for( int d=0; d<dimDomain; ++d )
+        f[r][d] = v[ d ] * u[ r ];
+  }
+
+  bool hasDirichletBoundary () const 
+  {
+    return true ;
+  }
+
+  bool isDirichletPoint( const DomainType& global ) const 
+  { 
+    return true ;
+  }
+
+protected:
+
+  /**
+   * @brief velocity calculation, is called by advection()
+   */
+  inline  void velocity(const EntityType& en,
+                        const double time,
+                        const DomainType& x,
+                        DomainType& v) const
+  {
+    problem_.velocity(en.geometry().global(x), v);
+  }
+
+public:
+
+  /**
+   * @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];
+  }
+
+  template <class T> 
+  T SQR( const T& a ) const 
+  {
+    return (a * a);
+  }
+
+  inline void eigenValues(const EntityType& en,
+                          const double time,
+                          const DomainType& x,
+                          const RangeType& u,
+                          RangeType& maxValue) const
+  {
+    DiffusionMatrixType K ; 
+    DomainType xgl = en.geometry().global( x );
+    problem_.K( xgl, K );
+
+    DomainType values ;
+    // calculate eigenvalues 
+    Dune::FMatrixHelp :: eigenValues( K, values );
+
+    maxValue = SQR(values[ dimDomain -1 ]) /values[0];
+    return ;
+
+    // take max eigenvalue 
+    maxValue = values.infinity_norm();
+  }
+
+  inline double lambdaK( const DiffusionMatrixType& K ) const 
+  {
+    DomainType values ;
+    // calculate eigenvalues 
+    Dune::FMatrixHelp :: eigenValues( K, values );
+
+    // value[ 0 ] is smallest ev 
+    return SQR(values[ dimDomain -1 ]) / values[ 0 ];
+  }
+
+  inline double penaltyFactor( const EntityType& inside,
+                               const EntityType& outside, 
+                               const double time, 
+                               const DomainType& xInside, 
+                               const RangeType& uLeft, 
+                               const RangeType& uRight ) const 
+  {
+    DiffusionMatrixType K ;
+    double betaK, betaInside, betaOutside ;
+    if( problem_.constantK() ) 
+    {
+      DiffusionMatrixType Kinside ; 
+      DiffusionMatrixType Koutside; 
+
+      const DomainType xglIn = inside.geometry().center();
+      problem_.K( xglIn , Kinside );
+      const DomainType xglOut = outside.geometry().center();
+      problem_.K( xglOut , Koutside );
+
+      K = Kinside ;
+      K += Koutside ; 
+      K *= 0.5 ; 
+    
+      betaInside  = lambdaK( Kinside );
+      betaOutside = lambdaK( Koutside );
+
+      betaK = lambdaK( K );
+    }
+    else 
+    {
+      const DomainType xgl = inside.geometry().global( xInside );
+      problem_.K( xgl , K );
+      
+      betaK = lambdaK( K );
+      betaInside = betaOutside = betaK;
+    }
+
+    const double jump = std::tanh( std::abs( betaInside - betaOutside ) );
+
+    // only for small values of betS apply betS in smooth cases 
+    const double betaN = std :: min( betaK , 1.0 );
+
+    // betS becomes 1 if the eigen values of both matrices are the same 
+    betaK = betaK * jump + (1.0 - jump) * betaN;
+
+    return betaK ;
+  }
+
+  inline double penaltyBoundary( const EntityType& inside,
+                                 const double time, 
+                                 const DomainType& xInside, 
+                                 const RangeType& uLeft ) const 
+  {
+    return penaltyFactor( inside, inside, time, xInside, uLeft, uLeft );
+  }
+
+  /**
+   * @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
+  {
+    // for constant K evalute at center (see Problem 4)
+    const DomainType xgl = ( problem_.constantK() ) ? 
+      en.geometry().center () : en.geometry().global( x )  ;
+
+    DiffusionMatrixType K ; 
+
+    // fill diffusion matrix 
+    problem_.K( xgl, K );
+
+    // apply diffusion 
+    for( int r =0; r<dimRange; ++r ) 
+      K.mv( jac[ r ] , A[ r ] );
+  }
+
+  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 0;
+  }
+
+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.;
+  }
+
+  /**
+   * @brief diffusion boundary flux
+   */
+  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 in this model" <<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
+  {
+    if (it.boundaryId() == 99) // Dirichlet zero boundary conditions
+    {
+      uRight = 0;
+    }
+    else
+    {
+      DomainType xgl = it.geometry().global( x );
+      problem_.g(xgl, uRight);
+    }
+  }
+  inline  void boundaryValue(const DomainType& xgl,
+                             RangeType& uRight) const
+  {
+    problem_.g(xgl, uRight);
+  }
+
+  const ProblemType& problem () const { return problem_; }
+
+ protected:
+  const ProblemType& problem_;
+  friend class Dune::UpwindFlux<PoissonModel>;
+};
+
+#endif

dune/fem-dg/test/poisson/passtraits.hh

+#ifndef DUNE_FEM_DG_PASSTRAITS_HH
+#define DUNE_FEM_DG_PASSTRAITS_HH
+
+// Dune includes
+#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+
+// Dune-Fem includes
+#include <dune/fem/space/finitevolume.hh>
+#include <dune/fem/space/discontinuousgalerkin.hh>
+#include <dune/fem/space/lagrange.hh>
+#include <dune/fem/space/padaptivespace.hh>
+#include <dune/fem/pass/localdg/discretemodel.hh>
+#include <dune/fem/function/adaptivefunction.hh>
+#include <dune/fem/quadrature/cachingquadrature.hh>
+#include <dune/fem-dg/operator/adaptation/adaptation.hh>
+#include <dune/fem/function/blockvectorfunction.hh>
+
+#if HAVE_PETSC
+#include <dune/fem/function/petscdiscretefunction.hh>
+#endif
+
+
+namespace Dune {
+
+  //PassTraits
+  //----------
+
+  template <class Model,int dimRange,int polOrd>
+  class PassTraits
+  {
+  public:
+    typedef typename Model :: Traits                                 ModelTraits;
+    typedef typename ModelTraits :: GridPartType                     GridPartType;
+    typedef typename GridPartType :: GridType                        GridType;
+    typedef typename GridType :: ctype                               ctype;
+    static const int dimDomain = Model :: Traits :: dimDomain;
+
+    template <class Grid, int dim >
+    struct FaceQuadChooser
+    {
+      typedef Dune::Fem::CachingQuadrature< GridPartType, 1 > Type;
+    };
+
+#if HAVE_UG
+    template <int dim>
+    struct FaceQuadChooser< const Dune::UGGrid< dim >, dim > 
+    {
+      typedef Dune::Fem::ElementQuadrature< GridPartType, 1 > Type;
+    };
+#endif
+
+    // CACHING
+    typedef typename FaceQuadChooser< const GridType, GridType::dimension > :: Type  FaceQuadratureType;
+    typedef Dune::Fem::CachingQuadrature< GridPartType, 0 >     VolumeQuadratureType; 
+
+    // Allow generalization to systems
+    typedef Dune::Fem::FunctionSpace< ctype, double, dimDomain, dimRange >      FunctionSpaceType;
+#if DGSCHEME
+#if ONB
+    #warning using DG space with ONB 
+    typedef Dune::Fem::DiscontinuousGalerkinSpace
+#elif LAG
+    #warning using DG space with Lagrange base functions
+    typedef Dune::Fem::LagrangeDiscontinuousGalerkinSpace
+#elif LEG
+    #warning using DG space with hierarich Legendre base functions
+    typedef Dune::Fem::HierarchicLegendreDiscontinuousGalerkinSpace
+#else
+#define USE_STRONG_BC
+    #warning using p-adaptive DG space
+    typedef Dune::Fem::PAdaptiveDGSpace
+#define PADAPTSPACE
+#endif
+#else
+#define USE_STRONG_BC
+#if LAGRANGESPACE
+    #warning using Lagrange space 
+    typedef Dune::Fem::LagrangeDiscreteFunctionSpace
+#else
+#define USE_STRONG_BC
+    #warning using p-adaptive Lagrange space
+    typedef Dune::Fem::PAdaptiveLagrangeSpace
+#define PADAPTSPACE
+#endif
+#endif
+       < FunctionSpaceType, GridPartType, polOrd, Dune::Fem::CachingStorage > DiscreteFunctionSpaceType;
+    
+#if WANT_ISTL
+  typedef Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpaceType > DestinationType;
+#elif WANT_PETSC
+  typedef Dune::Fem::PetscDiscreteFunction< DiscreteFunctionSpaceType >           DestinationType;
+#else
+  typedef Dune::Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceType >        DestinationType;
+#endif
+
+    // Indicator for Limiter
+    typedef Dune::Fem::FunctionSpace< ctype, double, dimDomain, 3> FVFunctionSpaceType;
+    typedef Dune::Fem::FiniteVolumeSpace<FVFunctionSpaceType,GridPartType, 0, Dune::Fem::SimpleStorage> IndicatorSpaceType;
+    typedef Dune::Fem::AdaptiveDiscreteFunction<IndicatorSpaceType> IndicatorType;
+
+    typedef AdaptationHandler< GridType, FunctionSpaceType >  AdaptationHandlerType ;
+  };
+
+} // end namespace Dune 
+#endif

dune/fem-dg/test/poisson/problemcreator.hh

+#undef ENABLE_MPI
+
+#ifndef FEMHOWTO_POISSONSTEPPER_HH
+#define FEMHOWTO_POISSONSTEPPER_HH
+#include <config.h>
+
+#include <dune/fem/main/codegen.hh>
+#include "passtraits.hh"
+
+// dune-grid includes
+#include <dune/grid/io/file/dgfparser/dgfparser.hh>
+
+// dune-fem includes
+#include <dune/fem/io/parameter.hh>
+
+// dune-fem-dg includes
+#include <dune/fem-dg/operator/fluxes/diffusionflux.hh>
+
+
+#include <dune/fem-dg/operator/dg/dgoperatorchoice.hh>
+#include <dune/fem-dg/assemble/primalmatrix.hh>
+#include "../common/inverseop.hh"
+
+// local includes
+#include "benchmarkproblem.hh"
+#include "models.hh"
+#include "estimator1.hh"
+
+#define NS_ELLIPTIC_OPERATOR 
+#ifndef TESTOPERATOR
+#define TESTOPERATOR 
+#endif
+
+template <class GridType> 
+struct ProblemGenerator 
+{
+  typedef Dune :: ProblemInterface<
+             Dune::Fem::FunctionSpace< double, double, GridType::dimension, DIMRANGE> >  ProblemType;
+
+  template <class GridPart>
+  struct Traits
+  {
+    typedef ProblemType  InitialDataType;
+    typedef PoissonModel< GridPart, InitialDataType > ModelType;
+    typedef Dune::NoFlux< ModelType > FluxType;
+    // choice of diffusion flux (see diffusionflux.hh for methods)
+    static const Dune :: DGDiffusionFluxIdentifier PrimalDiffusionFluxId 
+      = Dune :: method_general ;
+  };
+
+
+  static inline std::string diffusionFluxName()
+  {
+#if (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 
+    //GridPtr<GridType> gridptr = initialize< GridType >( description );
+    std::string filename = Dune::Fem::Parameter::getValue< std::string >(Dune::Fem::IOInterface::defaultGridKey(GridType :: dimension, false)); 
+  
+    // initialize grid
+    Dune::GridPtr< GridType > gridptr = initialize< GridType >( description );
+
+    std::string::size_type position = std::string::npos;
+    if( GridType::dimension == 2) 
+      position = filename.find("mesh3");
+    if( GridType::dimension == 3 )
+      position = filename.find("_locraf.msh" );
+
+    std::string::size_type locraf = filename.find("locrafgrid_");
+
+    std::string::size_type checkerboard = filename.find("heckerboard" );;
+
+    GridType& grid = *gridptr ;
+
+    const int refineelement = 1 ;
+
+    bool nonConformOrigin = Dune::Fem::Parameter::getValue< bool > ( "poisson.nonConformOrigin",false );
+
+    if ( nonConformOrigin )
+    {
+      std::cout << "Create local refined grid" << std::endl;
+      if( grid.comm().rank() == 0)
+      {
+        typedef typename GridType:: template Codim<0>::LeafIterator IteratorType;
+        IteratorType endit = grid.template leafend<0>();
+        for(IteratorType it = grid.template leafbegin<0>(); it != endit ; ++it)
+        {
+          const typename IteratorType :: Entity & entity = *it ;
+          const typename IteratorType :: Entity :: Geometry& geo = entity.geometry();
+          if (geo.center().two_norm() < 0.5)
+            grid.mark(refineelement, entity );
+        }
+      }
+    }
+    else if ( position != std::string::npos || 
+              locraf   != std::string::npos ) 
+    {
+      std::cout << "Create local refined grid" << std::endl;
+      if( grid.comm().rank() == 0)
+      {
+        typedef typename GridType:: template Codim<0>::LeafIterator IteratorType;
+        Dune::FieldVector<double,GridType::dimension> point(1.0);
+        if( locraf != std::string::npos ) point[ 0 ] = 3.0;
+        const int loops = ( GridType::dimension == 2 ) ? 2 : 1;
+        for (int i=0; i < loops; ++i)
+        {
+          point /= 2.0 ;
+
+          IteratorType endit = grid.template leafend<0>();
+          for(IteratorType it = grid.template leafbegin<0>(); it != endit ; ++it)
+          {
+            const typename IteratorType :: Entity & entity = *it ;
+            const typename IteratorType :: Entity :: Geometry& geo = entity.geometry();
+            bool inside = true ; 
+            const int corners = geo.corners();
+            for( int c = 0; c<corners; ++ c) 
+            {
+              for( int i = 0; i<GridType::dimension; ++i )
+              {
+                if( geo.corner( c )[ i ] - point[ i ] > 1e-8 ) 
+                {
+                  inside = false ; 
+                  break ;
+                }
+              }
+            }
+            if( inside ) grid.mark(refineelement, entity );
+          }
+        }
+      }
+    }
+    else if ( checkerboard != std::string::npos ) 
+    {
+      std::cout << "Create Checkerboard mesh" << std::endl;
+      int moduloNum = 32 ;
+      for( int i = 2; i<32; i *= 2 ) 
+      {
+        std::stringstream key; 
+        key << i << "x" << i << "x" << i;
+        std::string::size_type counter = filename.find( key.str() );
+        if( counter != std::string::npos ) 
+        {
+          moduloNum = i;
+          break ;
+        }
+      }
+
+      std::cout << "Found checkerboard cell number " << moduloNum << std::endl;
+
+      if( grid.comm().rank() == 0)
+      {
+        typedef typename GridType:: template Codim<0>::LeafIterator IteratorType;
+        IteratorType endit = grid.template leafend<0>();
+        int count = 1;
+        int modul = 1; // start with 1, such that the fine cube is at (0,0,0)
+        for(IteratorType it = grid.template leafbegin<0>(); it != endit ; ++it, ++ count )
+        {
+          const typename IteratorType :: Entity & entity = *it ;
+          if( count % 2 == modul ) 
+            grid.mark(refineelement, entity );
+
+          if( count % moduloNum == 0 ) 
+            modul = 1 - modul ;
+
+          if( count % (moduloNum * moduloNum) == 0 )
+            modul = 1 - modul ;
+        }
+      }
+    }
+
+    // adapt the grid 
+    grid.preAdapt();
+    grid.adapt();
+    grid.postAdapt();
+
+    //Dune :: GrapeGridDisplay< GridType > grape( grid );
+    //grape.display ();
+
+    return gridptr ;
+  }
+
+  static ProblemType* problem()
+  {
+    // choice of benchmark problem 
+    int probNr = Dune::Fem::Parameter::getValue< int > ( "femhowto.problem" );
+    return new Dune :: BenchmarkProblems< GridType > ( probNr );
+  }
+};
+
+#endif // FEMHOWTO_POISSONSTEPPER_HH