From d5ab0ae425f11ca83dbc88c1a0056420336df248 Mon Sep 17 00:00:00 2001
From: Robert Kloefkorn <robertk@posteo.net>
Date: Wed, 8 Oct 2014 21:16:10 +0200
Subject: [PATCH] copy poisson implementation from dglagrange. Not working yet.
---
dune/fem-dg/test/poisson/Makefile.am | 74 ++
dune/fem-dg/test/poisson/benchmark.cc | 1160 ++++++++++++++++++
dune/fem-dg/test/poisson/benchmark.hh | 1118 +++++++++++++++++
dune/fem-dg/test/poisson/benchmarkproblem.hh | 266 ++++
dune/fem-dg/test/poisson/models.hh | 500 ++++++++
dune/fem-dg/test/poisson/passtraits.hh | 105 ++
dune/fem-dg/test/poisson/problemcreator.hh | 197 +++
7 files changed, 3420 insertions(+)
create mode 100644 dune/fem-dg/test/poisson/Makefile.am
create mode 100644 dune/fem-dg/test/poisson/benchmark.cc
create mode 100644 dune/fem-dg/test/poisson/benchmark.hh
create mode 100644 dune/fem-dg/test/poisson/benchmarkproblem.hh
create mode 100644 dune/fem-dg/test/poisson/models.hh
create mode 100644 dune/fem-dg/test/poisson/passtraits.hh
create mode 100644 dune/fem-dg/test/poisson/problemcreator.hh
diff --git a/dune/fem-dg/test/poisson/Makefile.am b/dune/fem-dg/test/poisson/Makefile.am
new file mode 100644
index 00000000..31057579
--- /dev/null
+++ b/dune/fem-dg/test/poisson/Makefile.am
@@ -0,0 +1,74 @@
+# 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
diff --git a/dune/fem-dg/test/poisson/benchmark.cc b/dune/fem-dg/test/poisson/benchmark.cc
new file mode 100644
index 00000000..07caba2a
--- /dev/null
+++ b/dune/fem-dg/test/poisson/benchmark.cc
@@ -0,0 +1,1160 @@
+#ifndef FVCA5_BENCHMARK_CC
+#define FVCA5_BENCHMARK_CC
+
+#include <cmath>
+#include <dune/common/exceptions.hh>
+#include "benchmark.hh"
+
+template <int dim, class DomainField, class Field>
+class BenchMark_1 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ Field factor_[dim][dim];
+public:
+ virtual ~BenchMark_1() {}
+ BenchMark_1(Field globalShift, Field factor)
+ : globalShift_( globalShift )
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j)
+ {
+ if( i == j ) factor_[i][j] = 1.5;
+ else if( std::abs( i - j ) == 1 ) factor_[i][j] = 0.5;
+ else factor_[i][j] = 0;
+ }
+ }
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ // copy values
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j)
+ k[i][j] = factor_[i][j];
+ }
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ double x = arg[0];
+ double y = arg[1];
+
+ double uxx = -2.*y*(1-y)*16.;
+ double uxy = (-2.*x+1)*(-2.*y+1)*16.;
+ double uyy = -2.*x*(1-x)*16.;
+
+ return -( factor_[0][0] * uxx + factor_[0][1] * uxy + factor_[1][0] * uxy + factor_[1][1] * uyy);
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ Field val = 16.;
+ for(int i=0; i<dim; ++i)
+ val *= (x[i] - (x[i]*x[i]));
+ val += globalShift_ ;
+ return val;
+ }
+
+ virtual void gradExact(const DomainField arg[dim], Field grad[dim] ) const
+ {
+ double x = arg[0];
+ double y = arg[1];
+
+ grad[0] = (-2.*x+1)*y*(1-y)*16.;
+ grad[1] = (-2.*y+1)*x*(1-x)*16.;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark_1_2 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ Field factor_[dim][dim];
+public:
+ virtual ~BenchMark_1_2() {}
+ BenchMark_1_2(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j)
+ {
+ if( i == j ) factor_[i][j] = 1.5;
+ else factor_[i][j] = 0.5;
+ }
+ }
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_[i][i];
+ for(int j=0; j<i; ++j) k[i][j] = factor_[i][j];
+ for(int j=i+1; j<dim; ++j) k[i][j] = factor_[i][j];
+ }
+ }
+
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ double x1 = 1.-arg[0];
+ double y1 = 1.-arg[1];
+
+ double uxx = -y1*y1*sin(x1*y1) + 6.*x1* y1*y1;
+ double uxy = -x1*y1*sin(x1*y1) + cos(x1*y1) + 6.*y1* x1*x1;
+ double uyy = -x1*x1*sin(x1*y1) + 2.* x1*x1*x1;
+ return -(factor_[0][0] * uxx + factor_[0][1]*uxy + factor_[1][0]* uxy + factor_[1][1]*uyy);
+ }
+
+ virtual Field exact(const DomainField arg[dim]) const
+ {
+ double x1 = 1.-arg[0];
+ double y1 = 1.-arg[1];
+ return sin(x1*y1) + x1*x1*x1 * y1*y1;
+ }
+
+ virtual void gradExact(const DomainField arg[dim], Field grad[dim] ) const
+ {
+ double x1 = 1.-arg[0];
+ double y1 = 1.-arg[1];
+
+ grad[0] = -y1*cos(x1*y1) - 3.* (x1*y1)* (x1*y1);
+ grad[1] = -x1*cos(x1*y1) - 2.*y1* x1*x1*x1;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark_2 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ Field factor_[dim][dim];
+ const Field delta_;
+ const Field sqrtDelta_;
+ const Field x1_;
+ const Field x2_;
+public:
+ virtual ~BenchMark_2() {}
+ BenchMark_2(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , delta_(factor)
+ , sqrtDelta_( sqrt(delta_) )
+ , x1_ ( 8. * atan(1.) )
+ , x2_ ( x1_ / sqrtDelta_ )
+ {
+ factor_[0][0] = 1;
+ factor_[0][1] = factor_[1][0] = 0.0;
+ factor_[1][1] = delta_;
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_[i][i];
+ for(int j=0; j<i; ++j) k[i][j] = factor_[i][j];
+ for(int j=i+1; j<dim; ++j) k[i][j] = factor_[i][j];
+ }
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return sin(x1_* x[0])*exp(-x2_* x[1]);
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = x1_ * cos(x1_* x[0])*exp(-x2_ * x[1]);
+ grad[1] = -x2_ * exact(x);
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact(x);
+ // we have only neumann boundary here (see discretemodel.hh)
+ return true;
+ //return false;
+ }
+};
+
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark_3 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ Field factor_[dim][dim];
+ const Field delta_;
+ const Field pi_;
+ const Field cost_;
+ const Field sint_;
+public:
+ virtual ~BenchMark_3() {}
+ BenchMark_3(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , delta_(1e-3)
+ , pi_ ( 4. * atan(1.) )
+ , cost_ ( cos( 40. * pi_ / 180. ) )
+ , sint_ ( sqrt(1. - cost_*cost_) )
+ {
+ factor_[0][0] = cost_*cost_+delta_*sint_*sint_;
+ factor_[1][0] = factor_[0][1] = cost_*sint_*(1.-delta_);
+ factor_[1][1] = sint_*sint_+delta_*cost_*cost_;
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_[i][i];
+ for(int j=0; j<i; ++j) k[i][j] = factor_[i][j];
+ for(int j=i+1; j<dim; ++j) k[i][j] = factor_[i][j];
+ }
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ return 0.0;
+ }
+
+ double bndFunc (const double x,
+ const double lower,
+ const double upper,
+ const double lowval,
+ const double upval) const
+ {
+ if( x <= lower ) return lowval;
+ if( x >= upper ) return upval;
+
+ const double scale = (x - lower)/(upper - lower);
+ assert( scale >= 0 && scale <= 1 );
+
+ return (1. - scale) * lowval + scale * upval;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ const int bndId = this->getBoundaryId( x , x );
+ if( bndId == 0 )
+ {
+ return bndFunc(x[1],0.2,0.3,1,0.5);
+ }
+ else if( bndId == 1 )
+ {
+ return bndFunc(x[1],0.7,0.8,0.5,0);
+ }
+ else if( bndId == 2 )
+ {
+ return bndFunc(x[0],0.2,0.3,1,0.5);
+ }
+ else if( bndId == 3 )
+ {
+ return bndFunc(x[0],0.7,0.8,0.5,0);
+ }
+ return 0.5;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = 0.0;
+ grad[1] = 0.0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ // only dirichlet for this problem
+ return true;
+ }
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark_4 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+public:
+ virtual ~BenchMark_4() {}
+ BenchMark_4(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ {
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ for(int j=0; j<dim; ++j )
+ k[i][j] = 0;
+
+ if( omega1(x) )
+ {
+ k[0][0] = 1e2;
+ k[1][1] = 1e1;
+ if( dim > 2 )
+ {
+ k[2][2] = 1e3;
+ k[0][1] = 0.5;
+ k[1][0] = 0.5;
+ k[1][2] = 5;
+ k[2][1] = 5;
+ }
+ }
+ else
+ {
+ k[0][0] = 1e-2;
+ k[1][1] = 1e-3;
+ if( dim > 2 )
+ {
+ //k[0][1] = 1e-1;
+ //k[1][0] = 1e-1;
+ //k[1][2] = 1e-2;
+ //k[2][1] = 1e-2;
+ k[2][2] = 1e-2;
+ }
+ }
+ }
+
+ bool omega1(const DomainField x[dim]) const
+ {
+ if( dim == 2 )
+ {
+ if (x[0] <= 0.5)
+ {
+ int inty = int(10.0 * (x[1] + 0.15));
+ // if even then omega1 else omega2
+ return ((inty%2) == 0);
+ }
+ else
+ {
+ int inty = int(10.0 * x[1]);
+ // if even then omega1 else omega2
+ return ((inty%2) == 0);
+ }
+ }
+ else if ( dim == 3 )
+ {
+ if (x[0] <= 0.5)
+ {
+ int inty = int(16.0 * (x[1] + 0.0625));
+ // if even then omega1 else omega2
+ return ((inty%2) == 0);
+ }
+ else
+ {
+ int inty = int(16.0 * x[1]);
+ // if even then omega1 else omega2
+ return ((inty%2) == 0);
+ }
+ }
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ double val = (1.0 - x[0]);
+ if( dim > 2 )
+ val *= (1 - x[2]);
+ return val;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = 0.0;
+ grad[1] = 0.0;
+ grad[dim-1] = 0.0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark_5 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field delta_;
+ const Field pi;
+public:
+ virtual ~BenchMark_5() {}
+ BenchMark_5(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , delta_(1e-3)
+ , pi ( 4. * atan(1.) )
+ {
+ }
+
+ virtual bool constantLocalK () const { return false; }
+
+ virtual void K(const DomainField arg[dim], Field k[dim][dim] ) const
+ {
+ double x = arg[0];
+ double y = arg[1];
+ double rt = x*x+y*y;
+ k[0][0] = (y*y+delta_*x*x)/rt;
+ k[1][1] = (x*x+delta_*y*y)/rt;
+ k[1][0] = k[0][1] = -(1-delta_)*x*y/rt;
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ Field k[dim][dim];
+ K(arg,k);
+ double x = arg[0];
+ double y = arg[1];
+ double rt = x*x+y*y;
+
+ double ux = pi * cos(pi*x)*sin(pi*y);
+ double uy = pi * cos(pi*y)*sin(pi*x);
+
+ double f0 = sin(pi*x)*sin(pi*y)*pi*pi*(1+delta_)*(x*x+y*y)
+ + cos(pi*x)*sin(pi*y)*pi*(1.-3.*delta_)*x
+ + cos(pi*y)*sin(pi*x)*pi*(1.-3.*delta_)*y
+ + cos(pi*y)*cos(pi*x)*2.*pi*pi*(1.-delta_)*x*y;
+ double kxx = k[0][0];
+ double kyy = k[1][1];
+ double kxy = k[0][1];
+ return (f0+2.*(x*(kxx*ux+kxy*uy)+y*(kxy*ux+kyy*uy)))/rt;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return sin(pi*x[0])*sin(pi*x[1]);
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = pi * cos(pi*x[0])*sin(pi*x[1]);
+ grad[1] = pi * cos(pi*x[1])*sin(pi*x[0]);
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark_6 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field delta_;
+ const Field cost_;
+ const Field sint_;
+public:
+ virtual ~BenchMark_6() {}
+ BenchMark_6(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , delta_(0.2)
+ , cost_ ( 1./sqrt(1.+delta_*delta_) )
+ , sint_ ( delta_*cost_ )
+ {
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ double phi1 = x[1] - delta_ * (x[0] - .5) - .475;
+ double phi2 = phi1 - .05;
+
+ double alpha = 0.0;
+ double beta = 0.0;
+ if (phi1<0 || phi2>0)
+ {
+ alpha = 1.0;
+ beta = 0.1;
+ }
+ else
+ {
+ alpha = 100.0;
+ beta = 10.0;
+ }
+
+ k[0][0] = alpha*cost_*cost_+beta*sint_*sint_;
+ k[0][1] = k[1][0] = cost_*sint_*(alpha-beta);
+ k[1][1] = alpha*sint_*sint_+beta*cost_*cost_;
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return - x[0] - x[1] * delta_;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = -1.0;
+ grad[1] = -delta_;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ // we have neumann boundary here
+ return true;
+ }
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark_7 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field delta_;
+public:
+ virtual ~BenchMark_7() {}
+ BenchMark_7(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , delta_(0.2)
+ {
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ double phi1 = phi(x);
+ double phi2 = phi1 - .05;
+
+ int dom = domain( phi1, phi2 );
+ if( dom == 1 || dom == 3 )
+ {
+ k[0][0] = k[1][1] = 1;
+ k[1][0] = k[0][1] = 0;
+ }
+ else
+ {
+ k[0][0] = k[1][1] = 0.01;
+ k[1][0] = k[0][1] = 0;
+ }
+ }
+
+ double phi(const DomainField x[dim]) const
+ {
+ return x[1] - delta_ * (x[0] - .5) - .475;
+ }
+
+ int domain(const double phi1, const double phi2) const
+ {
+ if (phi1<0)
+ return 1;
+ else if (phi2<0)
+ return 2;
+ else
+ return 3;
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ double phi1 = phi(x);
+ double phi2 = phi1 - .05;
+
+ int dom = domain( phi1, phi2 );
+ if( dom == 1 )
+ {
+ return -phi1;
+ }
+ else if( dom == 2 )
+ {
+ return -phi1/.01;
+ }
+ else
+ {
+ return -phi2 - 5.;
+ }
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ double phi1 = phi(x);
+ double phi2 = phi1 - .05;
+ int dom = domain( phi1, phi2 );
+ if( dom == 1 || dom == 3 )
+ {
+ grad[0] = delta_;
+ grad[1] = -1.0;
+ }
+ else // if (dom == 2)
+ {
+ grad[0] = delta_/.01;
+ grad[1] = - 1./.01;
+ }
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+#if PROBLEM_8
+#include <dune/fem/gridpart/gridpart.hh>
+#include <dune/fem/space/fvspace.hh>
+#include <dune/fem/space/lagrangespace.hh>
+#include <dune/fem/function/adaptivefunction.hh>
+
+typedef LeafGridPart<GridType> ParamGridPartType;
+typedef FunctionSpace<double,double,GridType::dimension,1> ParamFunctionSpaceType;
+typedef LagrangeDiscreteFunctionSpace<ParamFunctionSpaceType,ParamGridPartType,2>
+ ParamDiscreteFunctionSpaceType;
+// typedef FiniteVolumeSpace<ParamFunctionSpaceType,ParamGridPartType,0>
+// ParamDiscreteFunctionSpaceType;
+typedef AdaptiveDiscreteFunction<ParamDiscreteFunctionSpaceType>
+ ParamDiscreteFunctionType;
+ParamDiscreteFunctionType* paramDiscreteFunction;
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark_8 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field delta_;
+public:
+ virtual ~BenchMark_8() {}
+ BenchMark_8(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , delta_(0.)
+ {
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ k[0][0] = k[1][1] = 1;
+ k[1][0] = k[0][1] = 0;
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ FieldVector<double,dim> p(0);
+ FieldVector<double,1> ret(0);
+ p[0] = arg[0];
+ p[1] = arg[1];
+
+ const ParamGridPartType& gridPart = paramDiscreteFunction->space().gridPart();
+ const typename ParamGridPartType::IndexSetType& index = gridPart.indexSet();
+ HierarchicSearch<GridType,ParamGridPartType::IndexSetType>
+ search(gridPart.grid(),index);
+ typename ParamDiscreteFunctionType::LocalFunctionType lf =
+ paramDiscreteFunction->localFunction((*(search.findEntity(p))));
+
+ lf.evaluate(search.findEntity(p)->geometry().local(p),ret);
+ /*
+ std::cout << search.findEntity(p)->geometry().local(p)
+ << " " << p
+ << " " << ret[0] << std::endl;
+ */
+ return ret[0];
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = 0.0;
+ grad[1] = 0.0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = 0.0;
+ return true;
+ }
+ virtual int getBoundaryId(const DomainField x[dim],
+ const DomainField n[dim]) const
+ {
+ if (std::abs(n[0]) > 1e-8) {
+ if (n[1] < 0.) return 0;
+ else return 1;
+ }
+ if (std::abs(n[0]) < 1e-8) {
+ if (n[1] < 0.) return 2;
+ else return 3;
+ }
+ abort();
+ return 0;
+ }
+};
+#endif
+
+template <int dim, class DomainField, class Field>
+class BenchMark_9 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field delta_;
+public:
+ virtual ~BenchMark_9() {}
+ BenchMark_9(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , delta_(0.)
+ {
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ k[0][0] = k[1][1] = 1;
+ k[1][0] = k[0][1] = 0;
+ }
+
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ return 0.;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = 0.0;
+ grad[1] = 0.0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = (x[0]<0.6)?1.:0.;
+ if (x[0]<=0 || x[0]>=1)
+ return false;
+ if (x[1]<=0 || x[1]>=1)
+ return false;
+ return true;
+ }
+};
+
+
+/////////////////////////////////////////////////////////////////////
+//
+// 3D Benchmark Problems
+//
+/////////////////////////////////////////////////////////////////////
+template <int dim, class DomainField, class Field>
+class BenchMark3d_1 : public DataFunctionIF<dim,DomainField,Field>
+{
+
+ Field K_[dim][dim];
+public:
+ virtual ~BenchMark3d_1() {}
+ BenchMark3d_1(Field globalShift, Field factor)
+ {
+ if( dim != 3 )
+ {
+ DUNE_THROW(Dune::NotImplemented,"Problem only implemented for dim=3");
+ }
+
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j)
+ {
+ if( i == j ) K_[i][j] = 1.0;
+ else if( std::abs( i - j ) == 1 )
+ {
+ K_[i][j] = 0.5;
+ }
+ else K_[i][j] = 0;
+ }
+ }
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ // just copy tensor
+ for(int i=0; i<dim; ++i)
+ for(int j=0; j<dim; ++j)
+ k[i][j] = K_[i][j];
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ return M_PI*M_PI*(3.0*sin(M_PI*x[0])*sin(M_PI*(x[1]+0.5))*sin(M_PI*(x[2]+(1.0/3.0)))
+ -cos(M_PI*x[0])*cos(M_PI*(x[1]+0.5))*sin(M_PI*(x[2]+(1.0/3.0)))
+ -sin(M_PI*x[0])*cos(M_PI*(x[1]+0.5))*cos(M_PI*(x[2]+(1.0/3.0))));
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return 1.0+sin(M_PI*x[0])*sin(M_PI*(x[1]+0.5))*sin(M_PI*(x[2]+(1.0/3.0)));
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = M_PI*cos(M_PI*x[0])*sin(M_PI*(x[1]+0.5))*sin(M_PI*(x[2]+(1.0/3.0)));
+ grad[1] = M_PI*sin(M_PI*x[0])*cos(M_PI*(x[1]+0.5))*sin(M_PI*(x[2]+(1.0/3.0)));
+ grad[2] = M_PI*sin(M_PI*x[0])*sin(M_PI*(x[1]+0.5))*cos(M_PI*(x[2]+(1.0/3.0)));
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark3d_3 : public DataFunctionIF<dim,DomainField,Field>
+{
+
+ Field K_[dim][dim];
+public:
+ virtual ~BenchMark3d_3() {}
+ BenchMark3d_3(Field globalShift, Field factor)
+ {
+ if( dim != 3 )
+ {
+ DUNE_THROW(Dune::NotImplemented,"Problem only implemented for dim=3");
+ }
+
+ for(int i=0; i<dim; ++i)
+ for(int j=0; j<dim; ++j)
+ K_[ i ][ j ] = 0 ;
+
+ // set diagonal
+ K_[0][0] = 1;
+ K_[1][1] = 1;
+ K_[2][2] = 1e3 ;
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j)
+ {
+ k[i][j] = K_[i][j];
+ }
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ return 1002.0*4.0*M_PI*M_PI*sin(2.0*M_PI*x[0])*sin(2.0*M_PI*x[1])*sin(2.0*M_PI*x[2]);
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return sin(2.0*M_PI*x[0])*sin(2.0*M_PI*x[1])*sin(2.0*M_PI*x[2]);
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ const Field pi2 = 2.0*M_PI;
+ grad[0] = pi2 * cos( pi2 * x[0]) * sin( pi2 * x[1]) * sin( pi2 * x[2]);
+ grad[1] = pi2 * sin( pi2 * x[0]) * cos( pi2 * x[1]) * sin( pi2 * x[2]);
+ grad[2] = pi2 * sin( pi2 * x[0]) * sin( pi2 * x[1]) * cos( pi2 * x[2]);
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark3d_4 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field tau_ ;
+public:
+ virtual ~BenchMark3d_4() {}
+ BenchMark3d_4(Field globalShift, Field factor)
+ : tau_( 0.2 )
+ {
+ if( dim != 3 )
+ {
+ DUNE_THROW(Dune::NotImplemented,"Problem only implemented for dim=3");
+ }
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ for(int j=0; j<dim; ++j )
+ k[i][j] = 0;
+
+ for(int i=0; i<2; ++ i ) k[i][i] = 1;
+ k[2][2] = tau_ ;
+ }
+
+ int getDomain(const DomainField x[dim]) const
+ {
+ return 0;
+ }
+
+ virtual Field rhs (const DomainField arg[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ Field u,u_x,u_y,u_z;
+ calculsolwell(x[0],x[1],x[2],u,u_x,u_y,u_z,1.0,0.0,0.0,1.0,0.0,tau_);
+ return u;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ Field u,u_x,u_y,u_z;
+ calculsolwell(x[0],x[1],x[2],u,u_x,u_y,u_z,1.0,0.0,0.0,1.0,0.0,tau_);
+ grad[0] = u_x;
+ grad[1] = u_y;
+ grad[dim-1] = u_z;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+
+ void calculsolwell (Field x, Field y, Field z,
+ Field& p, Field& px, Field& py, Field& pz,
+ Field lxx, Field lxy, Field lxz, Field lyy, Field lyz, Field lzz) const
+ {
+ Field alpha,mu0,
+ A11,A12,A13,A21,A22,A23,A31,A32,A33,
+ //XX,
+ YY,ZZ,s2,a2,a1,a0,sh,frho,mu,ch,frho_p,
+ GP_X,GP_Y,GP_Z;
+
+ alpha = 1.0;
+ mu0 = 1.85334252045523;
+
+ lxx = 1.;
+ lyy = 1.;
+ lzz = 0.2;
+ lxy = 0.;
+ lxz = 0.;
+ lyz = 0.;
+
+ A11 = 0.122859140982213;
+ A12 = 0.0;
+ A13 = -1.68776357811111;
+ A21 = 0.0;
+ A22 = 0.76472449133173;
+ A23 = 0.0;
+ A31 = 0.754790818120945;
+ A32 = 0.0;
+ A33 = 0.274721390893459;
+
+ a2 = 0.000603774740915493;
+
+
+ //-------------------------------------------------
+ // Eval solution
+ //-------------------------------------------------
+ //XX = A11*x+A12*y+A13*z;
+ YY = A21*x+A22*y+A23*z;
+ ZZ = A31*x+A32*y+A33*z;
+
+ // computation of s2 = sh^2, solution to
+
+ // a2 S^2 + a1 S + a0 = 0
+
+ a1 = a2 - ZZ*ZZ - YY*YY;
+
+ a0 = - YY*YY;
+
+ s2 = (-a1+sqrt(a1*a1 - 4.0*a0*a2))/(2.0*a2);
+
+ sh = sqrt(s2);
+
+ // argsh(w) = log(w+sqrt(w**2+1))
+
+ frho = sh+sqrt(s2+1.0);
+
+ mu = log(frho);
+
+ p = alpha *(mu - mu0);
+
+ //-------------------------------------------------
+ // Eval gradient
+ //-------------------------------------------------
+ ch = (frho+1.0/frho)*0.5;
+
+ frho_p = alpha /( (ZZ*ZZ*sh)/(ch*ch*ch) + (YY*YY*ch)/(sh*sh*sh) );
+
+ GP_X = 0.0;
+ GP_Y = frho_p * YY / ( sh * sh );
+ GP_Z = frho_p * ZZ / ( ch * ch );
+
+ px = A11*GP_X + A21*GP_Y + A31*GP_Z;
+ py = A12*GP_X + A22*GP_Y + A32*GP_Z;
+ pz = A13*GP_X + A23*GP_Y + A33*GP_Z;
+ }
+
+};
+
+
+template <int dim, class DomainField, class Field>
+class BenchMark3d_5 : public DataFunctionIF<dim,DomainField,Field>
+{
+ enum { numDomain = 4 };
+ Field tensor_[ numDomain ][ dim ];
+ Field alpha_[ numDomain ];
+ Field trace_[ numDomain ];
+ const Field pi2_ ;
+public:
+ virtual ~BenchMark3d_5() {}
+ BenchMark3d_5(Field globalShift, Field factor)
+ : pi2_( 2.0 * M_PI )
+ {
+ if( dim != 3 )
+ {
+ DUNE_THROW(Dune::NotImplemented,"Problem only implemented for dim=3");
+ }
+
+ {
+ Field (&tensor)[dim] = tensor_[ 0 ];
+ tensor[0] = 1.0;
+ tensor[1] = 10.0;
+ tensor[2] = 0.01;
+ }
+
+ {
+ Field (&tensor)[dim] = tensor_[ 1 ];
+ tensor[0] = 1.0;
+ tensor[1] = 0.1;
+ tensor[2] = 100.0;
+ }
+
+ {
+ Field (&tensor)[dim] = tensor_[ 2 ];
+ tensor[0] = 1.0;
+ tensor[1] = 0.01;
+ tensor[2] = 10.0;
+ }
+
+ {
+ Field (&tensor)[dim] = tensor_[ 3 ];
+ tensor[0] = 1.0;
+ tensor[1] = 100.0;
+ tensor[2] = 0.1;
+ }
+
+ alpha_[ 0 ] = 0.1;
+ alpha_[ 1 ] = 10.0;
+ alpha_[ 2 ] = 100.0;
+ alpha_[ 3 ] = 0.01;
+
+ for(int i=0; i<numDomain; ++i )
+ {
+ trace_[ i ] = 0;
+ for( int j=0; j<dim; ++ j)
+ trace_[ i ] += tensor_[ i ][ j ];
+ }
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ const int domain = getDomain( x );
+ assert( domain >= 0 && domain < numDomain );
+ for(int i=0; i<dim; ++i)
+ for(int j=0; j<dim; ++j )
+ k[i][j] = 0;
+ // set diagonal
+ for(int j=0; j<dim; ++j )
+ k[j][j] = tensor_[ domain ][ j ];
+ }
+
+ int getDomain(const DomainField x[dim]) const
+ {
+ if (x[1]<=0.5)
+ {
+ if (x[2]<=0.5) return 0; else return 3;
+ }
+ else
+ {
+ if (x[2]<=0.5) return 1; else return 2;
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ const int domain = getDomain(x);
+ assert( domain >= 0 && domain < numDomain );
+ Field result = 4.0*M_PI*M_PI*sin(pi2_*x[0])*sin(pi2_*x[1])*sin(pi2_*x[2]);
+ result *= alpha_[ domain ] * trace_[ domain ];
+ return result;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ const int domain = getDomain(x);
+ assert( domain >= 0 && domain < numDomain );
+ Field val = sin(pi2_*x[0])*sin(pi2_*x[1])*sin(pi2_*x[2]);
+ val *= alpha_[ domain ];
+ return val ;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ const int domain = getDomain(x);
+ assert( domain >= 0 && domain < numDomain );
+
+ const Field x_pi = pi2_*x[0] ;
+ const Field y_pi = pi2_*x[1] ;
+ const Field z_pi = pi2_*x[2] ;
+
+ grad[0] = pi2_ * cos( x_pi ) * sin( y_pi ) * sin( z_pi );
+ grad[1] = pi2_ * sin( x_pi ) * cos( y_pi ) * sin( z_pi );
+ grad[2] = pi2_ * sin( x_pi ) * sin( y_pi ) * cos( z_pi );
+
+ for(int i=0; i<dim; ++i )
+ grad[ i ] *= alpha_[ domain ];
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+#endif
diff --git a/dune/fem-dg/test/poisson/benchmark.hh b/dune/fem-dg/test/poisson/benchmark.hh
new file mode 100644
index 00000000..73670b2d
--- /dev/null
+++ b/dune/fem-dg/test/poisson/benchmark.hh
@@ -0,0 +1,1118 @@
+#ifndef ELLIPTPROBLEM_CC
+#define ELLIPTPROBLEM_CC
+
+#include <cmath>
+#include <set>
+#include <cassert>
+
+template<int dim, class DomainField, class Field>
+class DataFunctionIF
+{
+protected:
+ DataFunctionIF() {}
+public:
+ // destructor
+ virtual ~DataFunctionIF() {}
+
+ // returns true if K is constant on one element
+ virtual bool constantLocalK () const { return true; }
+
+ // diffusion tensor
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const = 0;
+ // right hand side
+ virtual Field rhs (const DomainField x[dim]) const = 0;
+ // right hand side
+ virtual Field rhs (const double time, const DomainField x[dim]) const
+ {
+ return rhs( x );
+ }
+ virtual void velocity(const DomainField x[dim], DomainField v[dim]) const
+ {
+ for(int i=0; i<dim; ++i)
+ v[i] = 0.;
+ }
+
+ // exact solution
+ virtual Field exact(const DomainField x[dim]) const = 0;
+
+ // exact solution (time dependent)
+ virtual Field exact(const double time, const DomainField x[dim]) const
+ {
+ return exact( x );
+ }
+
+ // exact gradient
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const = 0;
+
+ // boundary data
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const = 0;
+
+ // neumann boundary
+ virtual void neumann(const DomainField x[dim], Field grad[dim]) const
+ {
+ Field tmp[dim];
+ gradExact(x,tmp);
+ Field k[dim][dim];
+ K(x, k);
+ for(int i=0; i<dim; ++i)
+ {
+ grad[i] = 0;
+ for(int j=0; j<dim; ++j)
+ grad[i] += tmp[j] * k[i][j];
+ }
+ }
+
+ // boundary data
+ virtual bool boundaryDataFunction(const double time,
+ const DomainField x[dim],
+ Field & val) const
+ {
+ return boundaryDataFunction(x, val);
+ }
+
+ // neumann boundary
+ virtual void neumann(const double time,
+ const DomainField x[dim], Field grad[dim]) const
+ {
+ return neumann(x, grad);
+ }
+
+ virtual int getBoundaryId(const DomainField x[dim]) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ const int id = 2 * i;
+ if( x[i] <= 1e-10 )
+ {
+ return id;
+ }
+ else if ( x[i] >= 0.99999999 )
+ {
+ return id + 1;
+ }
+ }
+
+ //assert( false );
+ //abort();
+ return -1;
+ }
+ virtual int getBoundaryId(const DomainField x[dim],
+ const DomainField n[dim]) const
+ {
+ return getBoundaryId( x );
+ }
+
+ Field SQR( const Field& a ) const
+ {
+ return (a * a);
+ }
+};
+
+template <int dim, class DomainField, class Field>
+class AlbertaProblem : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field factor_;
+public:
+ virtual ~AlbertaProblem() {}
+ AlbertaProblem(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , factor_(1.0)
+ {
+ //assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j) k[i][j] = 0;
+ k[i][i] = factor_;
+ }
+ }
+
+ inline Field xSqr(const DomainField x[dim]) const
+ {
+ Field xsqr = 0.0;
+ for(int i=0; i<dim; ++i) xsqr += x[i] * x[i];
+ return xsqr;
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ const Field xsqr = xSqr( x );
+ return -(400.0 * xsqr - 20.0 * dim) * std :: exp( -10.0 * xsqr );
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return std :: exp( -10.0 * xSqr( x ) );
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ const Field factor = -20.0 * std :: exp( -10.0 * xSqr( x ) );
+ for(int i=0; i<dim; ++i)
+ grad[i] = x[i] * factor ;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ //if(x[0] <= 0.0) return false;
+ //if(x[1] <= 0.0) return false;
+ return true;
+ //return false;
+ }
+};
+
+
+template <int dim, class DomainField, class Field>
+class SinSinSin : public DataFunctionIF<dim,DomainField,Field>
+{
+ Field K_[dim][dim];
+public:
+ virtual ~SinSinSin() {}
+ SinSinSin(Field globalShift, Field factor)
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j)
+ {
+ if( i == j )
+ {
+ if( i == 0 )
+ K_[i][j] = 10.;//factor;
+ else
+ K_[i][j] = 0.1;//factor;
+ }
+ /*
+ else if( std::abs( i - j ) == 1 )
+ {
+ K_[i][j] = 0.5;
+ }
+ */
+ else K_[i][j] = 0;
+ }
+ }
+
+ /*
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j)
+ {
+ std::cout << K_[i][j] << " ";
+ }
+ std::cout << std::endl;
+ }
+ */
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j)
+ {
+ k[i][j] = K_[i][j];
+ }
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ Field sum = 0;
+ int comp[ dim - 1 ] ;
+ for(int i=0; i<dim; ++i)
+ {
+ comp[0] = i;
+ for(int j=0; j<dim; ++j)
+ {
+ comp[ dim - 2 ] = j;
+ sum += K_[j][i] * laplace( x, comp );
+ }
+ }
+ return -sum;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ const Field pi = 2.0 * M_PI;
+ Field val = 1.0;
+ for(int i=0; i<dim; ++i)
+ {
+ val *= sin( pi * x[i] );
+ }
+ return val;
+ }
+
+ double laplace(const DomainField x[dim], const int comp[dim - 1] ) const
+ {
+ const Field pi = 2.0 * M_PI;
+ Field val = pi * pi;
+
+ std::set<int> comps ;
+ for( int j=0; j<dim-1; ++j)
+ {
+ comps.insert( comp[j] );
+ }
+
+ if( comps.size() == 1 )
+ {
+ // add other components
+ for(int i=0; i<dim; ++i)
+ {
+ val *= sin( pi * x[ i ] );
+ }
+
+ // minus because sin'' = -sin
+ return -val;
+ }
+ else
+ {
+ for( int i=0; i<dim-1; ++i)
+ {
+ val *= cos( pi * x[ comp[i] ] );
+ }
+
+ for(int i=0; i<dim; ++i)
+ {
+ if( comps.find( i ) == comps.end() )
+ val *= sin( pi * x[ i ] );
+ }
+ return val;
+ }
+ }
+
+ double gradient(const DomainField x[dim], const int comp) const
+ {
+ const Field pi = 2.0 * M_PI;
+ Field val = pi * cos( pi * x[ comp ] );
+ // add other components
+ for(int j=1; j<dim; ++j)
+ {
+ val *= sin( pi * x[ (comp + j) % dim ] );
+ }
+ return val;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ grad[i] = gradient( x, i );
+ }
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ //if(x[0] <= 0.0) return false;
+ return true;
+ //return false;
+ }
+};
+
+
+template <int dim, class DomainField, class Field>
+class SinSin : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field factor_;
+public:
+ virtual ~SinSin() {}
+ SinSin(Field globalShift, Field factor)
+ : globalShift_(globalShift)
+ , factor_(factor)
+ {
+ //assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ Field sin_x = sin(2.0*M_PI*x[0]);
+ Field sin_y = sin(2.0*M_PI*x[1]);
+
+ Field val = 8.0 * M_PI * M_PI * sin_x * sin_y ;
+ val *= factor_;
+ return val;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ Field val = sin(2.0*M_PI*x[0]) * sin(2.0*M_PI*x[1]);
+ val += globalShift_;
+ return val;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = 2.0*M_PI*cos(2.0*M_PI*x[0])*sin(2.0*M_PI*x[1]);
+ grad[1] = 2.0*M_PI*sin(2.0*M_PI*x[0])*cos(2.0*M_PI*x[1]);
+
+ // initial grad with zero for 3d
+ for( int i=2; i<dim; ++i) grad[i] = 0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ //if(x[0] <= 0.0) return false;
+ return true;
+ //return false;
+ }
+};
+
+template <int dim, class DomainField, class Field>
+class CosCos : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field factor_;
+public:
+ virtual ~CosCos() {}
+ CosCos(Field globalShift, Field factor)
+ : globalShift_(globalShift)
+ , factor_(factor)
+ {
+ //assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ Field cos_x = cos(2.0*M_PI*x[0]);
+ Field cos_y = cos(2.0*M_PI*x[1]);
+
+ Field val = 8.0 * M_PI*M_PI* cos_x * cos_y ;
+ val *= factor_;
+ return val;
+ //return 0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ Field val = cos(2.0*M_PI*x[0]) * cos(2.0*M_PI*x[1]);
+ val += globalShift_;
+ return val;
+ //return x[1];
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = 2.0*M_PI*-sin(2.0*M_PI*x[0])*cos(2.0*M_PI*x[1]);
+ grad[1] = 2.0*M_PI*cos(2.0*M_PI*x[0])*-sin(2.0*M_PI*x[1]);
+ //grad[0] = 0.;
+ //grad[1] = 1.;
+
+ // initial grad with zero for 3d
+ for( int i=2; i<dim; ++i) grad[i] = 0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ //if(x[0] <= 0.0) return false;
+ return true;
+ }
+};
+
+//! problem from Castillo paper
+template <int dim, class DomainField, class Field>
+class CastilloProblem : public DataFunctionIF<dim,DomainField,Field>
+{
+ using DataFunctionIF<dim,DomainField,Field> :: SQR;
+ const Field globalShift_;
+ const Field factor_;
+public:
+ virtual ~CastilloProblem() {}
+ CastilloProblem(Field globalShift, Field factor)
+ : globalShift_(globalShift)
+ , factor_(factor)
+ {
+ assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ Field ret = 0.0;
+ Field tmp =-23+7*SQR(x[1])-24*x[0]+24*x[0]*SQR(x[1])+7*SQR(x[0])+9*SQR(x[0])*SQR(x[1])-24
+ *x[1]+24*x[1]*SQR(x[0]);
+
+ ret=-0.5*exp(0.75*(x[0]+x[1]));
+ ret *= tmp;
+ ret *= factor_;
+ return ret;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ Field val = 4.0 *(1.-SQR(x[0]))*(1.-SQR(x[1]))*exp(0.75*(x[0]+x[1]));
+ val += globalShift_;
+ return val;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ // not implemented yet
+ grad[0] = grad[1] = 0.0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+//! problem Einstringende Ecke
+template <class DomainField, class Field>
+class InSpringingCorner2d : public DataFunctionIF<2,DomainField,Field>
+{
+ enum { dim = 2 };
+ const Field globalShift_;
+ const Field factor_;
+ const Field lambda_;
+public:
+ virtual ~InSpringingCorner2d() {}
+ InSpringingCorner2d(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , factor_(1.0)
+ , lambda_(180./270.)
+ {
+ assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ double r2 = radius(x[0],x[1]);
+ double phi = argphi(x[0],x[1]);
+ return pow(r2,lambda_*0.5)*sin(lambda_*phi);
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ double r2=radius(x[0],x[1]);
+ double phi=argphi(x[0],x[1]);
+ double r2dx=2.*x[0];
+ double r2dy=2.*x[1];
+ double phidx=-x[1]/r2;
+ double phidy=x[0]/r2;
+ double lambdaPow = (lambda_*0.5)*pow(r2,lambda_*0.5-1.);
+ grad[0]= lambdaPow * r2dx * sin(lambda_ * phi)
+ + lambda_ * cos( lambda_ * phi) * phidx * pow(r2,lambda_ * 0.5);
+ grad[1]= lambdaPow * r2dy * sin(lambda_ * phi)
+ + lambda_ * cos( lambda_ * phi) * phidy * pow(r2,lambda_*0.5);
+ assert( grad[0] == grad[0] );
+ assert( grad[1] == grad[1] );
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+
+ private:
+ /** \brief proper implementation of atan(x,y)
+ */
+ inline double argphi(double x,double y) const
+ {
+ double phi=arg(std::complex<double>(x,y));
+ if (y<0) phi+=2.*M_PI;
+ return phi;
+ }
+
+ /** \brief implementation for the radius squared
+ * (^0.5 is part of function implementation)
+ */
+ inline double radius(double x, double y) const
+ {
+ double ret =0;
+ ret = x*x +y*y;
+ return ret;
+ }
+};
+
+//! problem Einstringende Ecke 3d
+template <int dim, class DomainField, class Field>
+class InSpringingCorner : public DataFunctionIF<dim,DomainField,Field>
+{
+ typedef InSpringingCorner2d<DomainField, Field> Corner2dType ;
+ Corner2dType corner2d_;
+ const Field factor_;
+
+public:
+ virtual ~InSpringingCorner() {}
+ InSpringingCorner(Field globalShift, Field factor)
+ : corner2d_(globalShift, factor),
+ factor_( factor )
+ {
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ return 0.0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ DomainField x2d[ 2 ];
+ if ( dim == 2 )
+ {
+ for( int i=0; i<dim; ++i ) x2d[ i ] = x[ i ];
+ return corner2d_.exact( x2d );
+ }
+ else
+ {
+ double sum = 0;
+ // ( x, y)
+ x2d[ 0 ] = -x[ 0 ]; x2d[ 1 ] = x[ 1 ];
+ sum += corner2d_.exact( x2d );
+ // (-y,-z)
+ //x2d[ 0 ] = -x[ 2 ]; x2d[ 1 ] = -x[ 1 ];
+ x2d[ 0 ] = -x[ 1 ]; x2d[ 1 ] = -x[ 2 ];
+ sum += corner2d_.exact( x2d );
+ // ( z,-x)
+ x2d[ 0 ] = -x[ 0 ]; x2d[ 1 ] = -x[ 2 ];
+ sum += corner2d_.exact( x2d );
+
+ //sum *= (x[ 0 ] + 1.0) / 2.0;
+ return sum ;
+ }
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ for( int i=0; i<dim; ++i ) grad[ i ] = 0;
+ if ( dim == 2 )
+ {
+ DomainField x2d[ 2 ];
+ for( int i=0; i<dim; ++i ) x2d[ i ] = x[ i ];
+ corner2d_.gradExact( x2d, &grad[ 0 ] );
+ }
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+//! problem Einstringende Ecke 3d
+template <int dim, class DomainField, class Field>
+class FicheraCorner : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field factor_;
+
+public:
+ virtual ~FicheraCorner() {}
+ FicheraCorner(Field globalShift, Field factor)
+ : globalShift_( globalShift ),
+ factor_( factor )
+ {
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ const double u = exact( x );
+ return -3.0/(4.0* u * u * u );
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ double xAbs = 0;
+ for( int i=0; i<dim; ++i )
+ xAbs += x[ i ] * x[ i ];
+ return std::pow( xAbs, 0.25 );
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ const double u = exact( x );
+ const double u3 = 2.0 * u * u * u ;
+ for( int i=0; i<dim; ++i ) grad[ i ] = x[ i ] / u3 ;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+//! problem single Hump
+template <int dim, class DomainField, class Field>
+class Hump : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field factor_;
+public:
+ virtual ~Hump() {}
+ Hump(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , factor_(1.0)
+ {
+ assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ double w=10.*x[0]*x[0]+10.*x[1];
+ double v=(x[0]-x[0]*x[0])*(x[1]-x[1]*x[1]);
+ double dwx = 20.*x[0];
+ double dwy = 10.;
+ double dwxx = 20.;
+ double dwyy = 0.;
+ double dvx = (1.-2.*x[0])*(x[1]-x[1]*x[1]);
+ double dvy = (1.-2.*x[1])*(x[0]-x[0]*x[0]);
+ double dvxx = -2.*(x[1]-x[1]*x[1]);
+ double dvyy = -2.*(x[0]-x[0]*x[0]);
+ Field grad[dim];
+ grad[0] = exp(w)*(dwx*v*v + 2.*v*dvx);
+ grad[1] = exp(w)*(dwy*v*v + 2.*v*dvy);
+ double dxx = dwx*grad[0] + exp(w)*(dwxx*v*v+dwx*2.*v*dvx + 2.*dvx*dvx+2.*v*dvxx);
+ double dyy = dwy*grad[1] + exp(w)*(dwyy*v*v+dwy*2.*v*dvy + 2.*dvy*dvy+2.*v*dvyy);
+ return -dxx-dyy;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ double w=10.*x[0]*x[0]+10.*x[1];
+ double v=(x[0]-x[0]*x[0])*(x[1]-x[1]*x[1]);
+ return exp(w)*v*v;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ double w=10.*x[0]*x[0]+10.*x[1];
+ double v=(x[0]-x[0]*x[0])*(x[1]-x[1]*x[1]);
+ double dwx = 20.*x[0];
+ double dwy = 10.;
+ double dvx = (1.-2.*x[0])*(x[1]-x[1]*x[1]);
+ double dvy = (1.-2.*x[1])*(x[0]-x[0]*x[0]);
+ grad[0] = exp(w)*(dwx*v*v + 2.*v*dvx);
+ grad[1] = exp(w)*(dwy*v*v + 2.*v*dvy);
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+//! problem Riviere-Bastian
+template <int dim, class DomainField, class Field>
+class RiviereProblem : public DataFunctionIF<dim,DomainField,Field>
+{
+ using DataFunctionIF<dim,DomainField,Field> :: SQR;
+ const Field globalShift_;
+ const Field factor_;
+public:
+ virtual ~RiviereProblem() {}
+ RiviereProblem(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , factor_(1.0)
+ {
+ assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ Field val = exact(x);
+ Field x_part = val * SQR(-2.0 * (x[0] - 0.5)) - 2.0 * val;
+ Field y_part = val * SQR(-2.0 * (x[1] - 0.5)) - 2.0 * val;
+ return -(x_part + y_part);
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ Field power = -( SQR( x[0] - 0.5 ) + SQR( x[1] - 0.5 ) );
+ return pow( M_E , power );
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ Field val = exact( x );
+ grad[0] = val * ( -2.0*x[0] + 1.0 );
+ grad[1] = val * ( -2.0*x[1] + 1.0 );
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+//! problem Riviere-Bastian
+template <int dim, class DomainField, class Field>
+class HeatProblem : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field factor_;
+public:
+ virtual ~HeatProblem() {}
+ HeatProblem(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , factor_(1.0)
+ {
+ assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ double scp(const DomainField x[dim]) const
+ {
+ double r2 = 0.0;
+ for(int i=0; i<dim; ++i) r2 += x[i]*x[i];
+ return r2;
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ return rhs(0.0, x);
+ }
+
+ virtual Field rhs (const double time, const DomainField x[dim]) const
+ {
+ double r2 = scp( x );
+
+ double ux = std::sin(M_PI*time)* std::exp(-10.0*r2);
+ double ut = M_PI*std::cos(M_PI*time)*std::exp(-10.0*r2);
+ return(ut - (400.0*r2 - 20.0*dim)*ux);
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ return exact(0.0, x);
+ }
+
+ virtual Field exact(const double time, const DomainField x[dim]) const
+ {
+ double r2 = scp( x );
+ return(std::sin( M_PI * time) * std::exp( -10.0 * r2));
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = grad[1] = 0.0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ assert( false );
+ val = exact( x );
+ return true;
+ }
+
+ virtual bool boundaryDataFunction(const double time,
+ const DomainField x[dim],
+ Field & val) const
+ {
+ val = exact( time, x );
+ return true;
+ }
+};
+
+template <int dim, class DomainField, class Field>
+class BoundaryLayerProblem : public DataFunctionIF<dim,DomainField,Field>
+{
+public:
+ virtual ~BoundaryLayerProblem() {}
+ BoundaryLayerProblem(Field globalShift, Field factor)
+ : eps_(factor)
+ {
+ assert( eps_ > 0 );
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ for(int j=0; j<dim; ++j) k[i][j] = 0;
+ k[i][i] = 1.; // eps_;
+ }
+ }
+
+ void velocity(const DomainField x[dim], DomainField v[dim]) const
+ {
+ for(int i=0; i<dim; ++i)
+ v[i] = 1./eps_;
+ }
+
+ virtual Field rhs (const DomainField x[dim]) const
+ {
+ return 0;
+ }
+
+ virtual Field exact(const DomainField x[dim]) const
+ {
+ Field ret = 1;
+ for (int i=0;i<dim;++i)
+ ret *= u1( x[i] );
+ return ret;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ for (int i=0;i<dim;++i)
+ {
+ grad[ i ] = 1.;
+ for (int j=0;j<dim;++j)
+ grad[ i ] *= (j==i)? du1( x[i] ):u1( x[i] );
+ }
+ return;
+ }
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+private:
+ double u1(double x) const
+ {
+ return (exp(x/eps_)-1.) / (exp(1./eps_)-1.);
+ }
+ double du1(double x) const
+ {
+ return (exp(x/eps_)/eps_) / (exp(1./eps_)-1.);
+ }
+ const double eps_;
+};
+
+
+//! problem CurvedRidges (deal.II step-14)
+template <int dim, class DomainField, class Field>
+class CurvedRidges : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field factor_;
+public:
+ virtual ~CurvedRidges() {}
+ CurvedRidges(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , factor_(1.0)
+ {
+ assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ virtual Field rhs (const DomainField p[dim]) const
+ {
+ double q = p[ 0 ];
+ for (unsigned int i=1; i<dim; ++i)
+ {
+ q += std::sin(10*p[ i ]+5*p[ 0 ]*p[ 0 ]);
+ }
+ const double u = std::exp(q);
+ double t1 = 1, t2 = 0, t3 = 0;
+ for (unsigned int i=1; i<dim; ++i)
+ {
+ t1 += std::cos(10*p[ i ]+5*p[ 0 ]*p[ 0 ]) * 10 * p[ 0 ];
+ t2 += 10*std::cos(10*p[ i ]+5*p[ 0 ]*p[ 0 ]) -
+ 100*std::sin(10*p[ i ]+5*p[ 0 ]*p[ 0 ]) * p[ 0 ]*p[ 0 ];
+ t3 += 100*std::cos(10*p[ i ]+5*p[ 0 ]*p[ 0 ])*std::cos(10*p[ i ]+5*p[ 0 ]*p[ 0 ]) -
+ 100*std::sin(10*p[ i ]+5*p[ 0 ]*p[ 0 ]);
+ }
+ t1 = t1*t1;
+
+ return -u*(t1+t2+t3);
+ }
+
+ virtual Field exact(const DomainField p[dim]) const
+ {
+ double q = p[ 0 ];
+ for (unsigned int i=1; i<dim; ++i)
+ {
+ q += std::sin(10*p[ i ]+5*p[ 0 ]*p[ 0 ]);
+ }
+ const double exponential = std::exp(q);
+ return exponential;
+ }
+
+ virtual void gradExact(const DomainField p[dim], Field grad[dim] ) const
+ {
+ double u = exact(p);
+ grad[0] = 1.;
+ for (unsigned int i=1; i<dim; ++i)
+ grad[0] += std::cos(10*p[ i ]+5*p[ 0 ]*p[ 0 ]) * 10. * p[0];
+ grad[0] *= u;
+ for (int i=1;i<dim;++i)
+ {
+ grad[i] = std::cos(10*p[ i ]+5*p[ 0 ]*p[ 0 ]) * 10.;
+ grad[i] *= u;
+ }
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+//! problem CurvedRidges (deal.II step-14)
+template <int dim, class DomainField, class Field>
+class Excercise2_3 : public DataFunctionIF<dim,DomainField,Field>
+{
+ const Field globalShift_;
+ const Field factor_;
+ Field point_ [ dim ];
+public:
+ virtual ~Excercise2_3() {}
+ Excercise2_3(Field globalShift, Field factor)
+ : globalShift_(0.0)
+ , factor_(1.0)
+ {
+ for(int i=0; i<dim; ++i) point_[ i ] = 0.75;
+ assert(dim == 2);
+ }
+
+ virtual void K(const DomainField x[dim], Field k[dim][dim] ) const
+ {
+ for(int i=0; i<dim; ++i)
+ {
+ k[i][i] = factor_;
+ for(int j=0; j<i; ++j) k[i][j] = 0;
+ for(int j=i+1; j<dim; ++j) k[i][j] = 0;
+ }
+ }
+
+ virtual Field rhs (const DomainField p[dim]) const
+ {
+ return 1.0;
+ }
+
+ virtual Field exact(const DomainField p[dim]) const
+ {
+ // exact solution not known
+ return 0.0;
+ }
+
+ virtual void gradExact(const DomainField x[dim], Field grad[dim] ) const
+ {
+ grad[0] = 0;
+ grad[1] = 0;
+ }
+
+ virtual bool boundaryDataFunction(const DomainField x[dim], Field & val) const
+ {
+ val = exact( x );
+ return true;
+ }
+};
+
+
+
+
+#include "benchmark.cc"
+#endif
diff --git a/dune/fem-dg/test/poisson/benchmarkproblem.hh b/dune/fem-dg/test/poisson/benchmarkproblem.hh
new file mode 100644
index 00000000..ff3960f4
--- /dev/null
+++ b/dune/fem-dg/test/poisson/benchmarkproblem.hh
@@ -0,0 +1,266 @@
+#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__*/
+
diff --git a/dune/fem-dg/test/poisson/models.hh b/dune/fem-dg/test/poisson/models.hh
new file mode 100644
index 00000000..d290b5ae
--- /dev/null
+++ b/dune/fem-dg/test/poisson/models.hh
@@ -0,0 +1,500 @@
+#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
diff --git a/dune/fem-dg/test/poisson/passtraits.hh b/dune/fem-dg/test/poisson/passtraits.hh
new file mode 100644
index 00000000..1a135782
--- /dev/null
+++ b/dune/fem-dg/test/poisson/passtraits.hh
@@ -0,0 +1,105 @@
+#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
diff --git a/dune/fem-dg/test/poisson/problemcreator.hh b/dune/fem-dg/test/poisson/problemcreator.hh
new file mode 100644
index 00000000..62be0b01
--- /dev/null
+++ b/dune/fem-dg/test/poisson/problemcreator.hh
@@ -0,0 +1,197 @@
+#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
--
GitLab