Merge branch 'feature/stdcomplex-test' into 'master'

Feature/stdcomplex test

See merge request !2

dune/pdelab-systemtesting/CMakeLists.txt

 #install headers
 install(FILES pdelab-systemtesting.hh DESTINATION include/dune/pdelab-systemtesting)
-
+add_subdirectory(helmholtz_stdcomplex)

dune/pdelab-systemtesting/helmholtz_stdcomplex/CMakeLists.txt

+add_executable("helmholtz" helmholtz.cc)
+add_dune_ug_flags(helmholtz)
+add_dune_superlu_flags(helmholtz)
+add_dune_umfpack_flags(helmholtz)
+target_link_dune_default_libraries(helmholtz)

dune/pdelab-systemtesting/helmholtz_stdcomplex/helmholtz.cc

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+/** \file
+
+    \brief Solve elliptic problem in unconstrained spaces with conforming finite elements
+*/
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#else
+#warning ERROR NO CONFIG.H
+#endif
+#include<math.h>
+#include<iostream>
+#include<fstream>
+#include<vector>
+#include<map>
+#include<string>
+#include<typeinfo>
+
+#include <random>
+
+
+#include <sys/types.h>  //for mkdir chdir etc
+#include <sys/stat.h>
+#include <unistd.h>
+
+
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/exceptions.hh>
+#include<dune/common/fvector.hh>
+#include<dune/common/static_assert.hh>
+#include<dune/common/timer.hh>
+#include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include<dune/grid/io/file/gmshreader.hh>
+#include<dune/grid/yaspgrid.hh>
+#if HAVE_ALBERTA
+#include<dune/grid/albertagrid.hh>
+#include <dune/grid/albertagrid/dgfparser.hh>
+#endif
+#if HAVE_UG
+#include<dune/grid/uggrid.hh>
+#else
+#error UGERR!
+#endif
+#if HAVE_ALUGRID
+#include<dune/grid/alugrid.hh>
+#include<dune/grid/io/file/dgfparser/dgfalu.hh>
+#include<dune/grid/io/file/dgfparser/dgfparser.hh>
+#endif
+#include<dune/istl/bvector.hh>
+#include<dune/istl/operators.hh>
+#include<dune/istl/solvers.hh>
+#include<dune/istl/preconditioners.hh>
+#include<dune/istl/io.hh>
+
+#include<complex>
+#define SUPERLU_NTYPE 3 //needed for complex superlu, because superlu is in c where is no templates
+#include<dune/istl/superlu.hh>
+// SUPERLU_NTYPE==0 float
+// SUPERLU_NTYPE==1 double
+// SUPERLU_NTYPE==2 std::complex<float>
+// SUPERLU_NTYPE==3 std::complex<double>
+
+
+
+#include <dune/istl/umfpack.hh>
+
+
+#include<dune/pdelab/newton/newton.hh>
+#include<dune/pdelab/finiteelementmap/p0fem.hh>
+#include<dune/pdelab/finiteelementmap/pkfem.hh>
+#include<dune/pdelab/finiteelementmap/qkfem.hh>
+//#include<dune/pdelab/finiteelementmap/rannacher_turek2dfem.hh>
+#include<dune/pdelab/constraints/common/constraints.hh>
+#include<dune/pdelab/constraints/conforming.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
+#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
+#include<dune/pdelab/gridfunctionspace/interpolate.hh>
+#include<dune/pdelab/gridfunctionspace/vtk.hh>
+#include<dune/pdelab/common/function.hh>
+#include<dune/pdelab/common/vtkexport.hh>
+#include<dune/pdelab/backend/istlvectorbackend.hh>
+#include<dune/pdelab/backend/istl/bcrsmatrixbackend.hh>
+#include<dune/pdelab/backend/istlsolverbackend.hh>
+//#include<dune/pdelab/backend/seqistlsolverbackend.hh>
+#include<dune/pdelab/stationary/linearproblem.hh>
+#include<dune/pdelab/gridoperator/gridoperator.hh>
+#include<dune/pdelab/instationary/onestep.hh> //Filename helper
+#include <dune/common/parametertreeparser.hh>
+
+
+
+
+
+#include<dune/geometry/referenceelements.hh>
+#include<dune/geometry/quadraturerules.hh>
+#include<dune/pdelab/common/geometrywrapper.hh>
+#include<dune/pdelab/localoperator/defaultimp.hh>
+#include<dune/pdelab/localoperator/pattern.hh>
+#include<dune/pdelab/localoperator/flags.hh>
+#include <dune/pdelab/boilerplate/pdelab.hh>
+
+
+
+#include"parameters_planewave.hh"
+#include"parameters_sphericalwave.hh"
+#include"helmholtz_bcextension.hh"
+#include"helmholtz_bcanalytic.hh"
+#include"helmholtz_operator.hh"
+#include"helmholtz_Qk.hh"
+
+
+
+
+
+
+
+//===============================================================
+// Main program with grid setup
+//===============================================================
+
+#include"helmholtz_main.hh"

dune/pdelab-systemtesting/helmholtz_stdcomplex/helmholtz.cfg

+[grid]
+level = 8
+[problem_parameter]
+polynomialdegree = 1
+omega = 80.
+[norm]
+errornorm = L2 #H1
+[solvers]
+solver = "UMFPACK" #, "SuperLU", "GMRESILU0", "GMRESILU1", "GMRESSGS", "BiCGSILU0", "BiCGSILU1", "BiCGSSGS"
\ No newline at end of file

dune/pdelab-systemtesting/helmholtz_stdcomplex/helmholtz_Qk.hh

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+
+
+#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_QK_HH
+#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_QK_HH
+
+
+
+/** \file
+
+    \brief Construct GFS and solve problem defined in parameter
+*/
+
+template<int k, class GV, class PARAM>
+void helmholtz_Qk (const GV& gv, PARAM& param, std::string& errornorm, std::string solver)
+{
+  // <<<1>>> Choose domain and range field type
+  typedef typename GV::Grid::ctype Coord;
+      typedef double R;
+      typedef std::complex<R> C;
+      typedef C RF;
+
+
+
+  // <<<2>>> Make grid function space
+  typedef Dune::PDELab::QkLocalFiniteElementMap<GV,Coord,C,k> FEM;
+  //typedef Dune::PDELab::PkLocalFiniteElementMap<GV,Coord,C,k> FEM;
+  FEM fem(gv);
+  // typedef Dune::PDELab::NoConstraints CON;  // !!!!! NO CONTRAINTS
+  typedef Dune::PDELab::ConformingDirichletConstraints CON; // constraints class
+  typedef Dune::PDELab::ISTLVectorBackend<> VBE;
+  typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE> GFS;
+  GFS gfs(gv,fem);
+  gfs.name("solution");
+  //BCTypeParam bctype; // boundary condition type
+
+  typedef typename GFS::template ConstraintsContainer<C>::Type CC;
+  CC cc;
+  Dune::PDELab::constraints( param, gfs, cc ); // assemble constraints
+  gfs.update();
+  std::cout << "constrained dofs=" << cc.size() << " of " << gfs.globalSize() << std::endl;
+
+
+  // <<<3>>> Make grid operator
+  typedef HelmholtzLocalOperator<PARAM > LOP;
+
+  LOP lop( param, 2*k );
+  typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
+  // Structured 2D grid, Q1 finite elements => 9-point stencil / Q2 => 25
+  MBE mbe(k == 1 ? 9 : 25);
+
+  typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,C,C,C,CC,CC> GO;
+  //GO go(gfs,gfs,lop,mbe);
+  GO go(gfs,cc,gfs,cc,lop,mbe);
+
+  // How well did we estimate the number of entries per matrix row?
+  // => print Jacobian pattern statistics
+  //typename GO::Traits::Jacobian jac(go);
+  //std::cout << jac.patternStatistics() << std::endl;
+
+
+
+  typedef typename GO::Traits::Domain U;
+  C zero(0.);
+  U u(gfs,zero);                                      // initial value
+  //ure.base().bar(); //determine type of u
+  typedef BCExtension<PARAM > G;                      // boundary value = extension
+  G g(gv,param);
+  Dune::PDELab::interpolate(g,gfs,u);                // interpolate coefficient vector
+  // Dune::PDELab::interpolate(param,gfs,u);
+  Dune::PDELab::set_nonconstrained_dofs(cc,0.0,u);
+
+
+
+
+
+  // <<<4>>> Select a linear solver backend
+  // typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_SSOR LS;
+  //typedef Dune::PDELab::ISTLBackend_SEQ_SuperLU LS ;
+  //typedef Dune::PDELab::ISTLBackend_SEQ_CG_SSOR LS;
+  // typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LS;
+  // LS ls(5000,0);
+  // typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,U> SLP;
+  // SLP slp(go,ls,u,1e-10);
+  // slp.apply();
+
+  //-------------------------------------------------------
+  // we do not use the backend, we call solver directly instead
+
+
+  typedef typename Dune::PDELab::BackendVectorSelector<GFS,RF>::Type V;
+  typedef typename GO::Jacobian M;
+
+  typedef typename M::BaseT ISTLM;
+  typedef typename V::BaseT ISTLV;
+
+
+  M m(go,0.);
+  //std::cout << m.patternStatistics() << std::endl;
+  go.jacobian(u,m);
+  V r(gfs);
+  r = 0.0;
+  go.residual(u,r);
+
+
+
+
+
+    if(solver == "UMFPACK") {
+      Dune::UMFPack<ISTLM> solver(m.base(), 0);
+      r *= -1.0; // need -residual
+      //u = r;
+      // u = 0;
+      Dune::InverseOperatorResult stat;
+      solver.apply(u.base(),r.base(),stat);
+    }
+    if(solver == "SuperLU") {
+      Dune::SuperLU<ISTLM> solver(m.base(), 0);
+      r *= -1.0; // need -residual
+      //u = r;
+      // u = 0;
+      Dune::InverseOperatorResult stat;
+      solver.apply(u.base(),r.base(),stat);
+    }
+
+
+
+    Dune::InverseOperatorResult stat;
+
+    if (solver == "GMRESILU0") {
+      Dune::SeqILU0<ISTLM,ISTLV,ISTLV> ilu0(m.base(),1.0);
+      Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(m.base());
+      Dune::RestartedGMResSolver<ISTLV> solver(opa, ilu0, 1E-7, 5000, 5000, 0);
+      r *= -1.0; // need -residual
+      //u = r;
+      // u = 0;
+      solver.apply(u.base(),r.base(),stat);
+      std::cout<<"Iterations: "<< stat.iterations<<std::endl;
+      std::cout<<"Time: "<<  stat.elapsed<< std::endl;
+    }
+    if (solver == "GMRESILU1") {
+      Dune::SeqILUn<ISTLM,ISTLV,ISTLV> ilun(m.base(), 1, 1.0);
+      Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(m.base());
+      Dune::RestartedGMResSolver<ISTLV> solver(opa, ilun, 1E-7, 5000, 5000, 0);
+      r *= -1.0; // need -residual
+      //u = r;
+      // u = 0;
+      solver.apply(u.base(),r.base(),stat);
+      std::cout<<"Iterations: "<< stat.iterations<<std::endl;
+      std::cout<<"Time: "<<  stat.elapsed<< std::endl;
+    }
+    if (solver == "GMRESSGS") {
+      Dune::SeqSSOR<ISTLM,ISTLV,ISTLV> ssor(m.base(), 3, 1.); //ssor with \omega = 1 is SGS (symmetric gauss seidel)
+      Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(m.base());
+      Dune::RestartedGMResSolver<ISTLV> solver(opa, ssor, 1E-7, 5000, 5000, 0);
+      r *= -1.0; // need -residual
+      //u = r;
+      // u = 0;
+      solver.apply(u.base(),r.base(),stat);
+      std::cout<<"Iterations: "<< stat.iterations<<std::endl;
+      std::cout<<"Time: "<<  stat.elapsed<< std::endl;
+    }
+    if (solver == "BiCGSILU0") {
+      Dune::SeqILU0<ISTLM,ISTLV,ISTLV> ilu0(m.base(),1.0);
+      Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(m.base());
+      Dune::BiCGSTABSolver<ISTLV> solver(opa,ilu0,1E-7,20000, 0);
+      // solve the jacobian system
+      r *= -1.0; // need -residual
+      //u = r;
+      // u = 0;
+      solver.apply(u.base(),r.base(),stat);
+      std::cout<<"Iterations: "<< stat.iterations<<std::endl;
+      std::cout<<"Time: "<<  stat.elapsed<< std::endl;
+    }
+
+    //this solver breaks down at 4000 DOFs by solving the helmholtz eq
+    if (solver == "BiCGSILU1") {
+      Dune::SeqILUn<ISTLM,ISTLV,ISTLV> ilun(m.base(), 1, 1.0);
+      Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(m.base());
+      Dune::BiCGSTABSolver<ISTLV> solver(opa,ilun,1E-7,20000, 0);
+      // solve the jacobian system
+      r *= -1.0; // need -residual
+      // //u = r;
+      // u = 0;
+      solver.apply(u.base(),r.base(),stat);
+      std::cout<<"Iterations: "<< stat.iterations<<std::endl;
+      std::cout<<"Time: "<<  stat.elapsed<< std::endl;
+    }
+    if (solver == "BiCGSSGS") {
+      Dune::SeqSSOR<ISTLM,ISTLV,ISTLV> ssor(m.base(), 3, 1.); //ssor with \omega = 1 is SGS (symmetric gauss seidel)
+      Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV> opa(m.base());
+      Dune::BiCGSTABSolver<ISTLV> solver(opa,ssor,1E-7,20000, 0);
+      r *= -1.0; // need -residual
+      // //u = r;
+      // u = 0;
+      solver.apply(u.base(),r.base(),stat);
+      std::cout<<"Iterations: "<< stat.iterations<<std::endl;
+      std::cout<<"Time: "<<  stat.elapsed<< std::endl;
+    }
+
+
+
+
+  // // // compare helmholtz solution to analytic helmholtz solution
+  U ua(gfs, zero);                                      // initial value
+  typedef BCAnalytic<PARAM> Ga;                      // boundary value = extension
+  Ga ga(gv, param);
+  Dune::PDELab::interpolate(ga,gfs,ua);                // interpolate coefficient vector
+
+
+  // // Make a real grid function space
+  typedef Dune::PDELab::QkLocalFiniteElementMap<GV,Coord,R,k> FEMr;
+  typedef Dune::PDELab::GridFunctionSpace<GV,FEMr,CON,VBE> GFSr;
+  FEMr femr(gv);
+  GFSr gfsr(gv,femr);
+
+  //create real analytic solution vectors
+  typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type reu(gfsr, 0.);  // real part u_h
+  typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type imu(gfsr, 0.);  // imag part u_h
+  typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type reua(gfsr, 0.); // real part analytic solution
+  typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type imua(gfsr, 0.); // imega part analytic solution
+  // typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type reue(gfsr, 0.); // real part error
+  // typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type imue(gfsr, 0.); // imag part error
+
+  typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type::iterator reut = reu.begin();
+  typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type::iterator imut = imu.begin();
+  typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type::iterator reuat = reua.begin();
+  typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type::iterator imuat = imua.begin();
+  // typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type::iterator reuet = reue.begin();
+  // typename Dune::PDELab::BackendVectorSelector<GFSr,R>::Type::iterator imuet = imue.begin();
+
+
+  typename U::const_iterator ut = u.begin();
+  typename U::const_iterator uat = ua.begin();
+
+  // //compute error i.e. |ua - u| separated into real and imag part
+  while( ut != u.end() ) {
+
+    *reut = std::real(*ut);
+    *imut = std::imag(*ut);
+    *reuat = std::real(*uat);
+    *imuat = std::imag(*uat);
+    // *reuet = std::abs(*reut - *reuat) ;
+    // *imuet = std::abs(*imut - *imuat) ;
+
+    // std::cout<<h<<"\t"<<*reut - (*reuat)<<std::endl;
+    //   ++h;
+
+    ++ut;
+    ++uat;
+
+    ++reut;
+    ++imut;
+    ++reuat;
+    ++imuat;
+    // ++reuet;
+    // ++imuet;
+
+  }
+
+
+  //<<<7>>> graphical output
+
+  // output u_h
+  //Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,k == 1 ? 0 : 3);
+  Dune::SubsamplingVTKWriter<GV> vtkwriter(gv, 0);
+
+
+  gfsr.name("real");
+  Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfsr,reu);
+
+  gfsr.name("imaginary");
+  Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfsr,imu);
+
+  gfsr.name("real_analytic");
+  Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfsr,reua);
+  gfsr.name("imaginary_analytic");
+  Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfsr,imua);
+
+
+  // gfsr.name("real_error");
+  // Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfsr,reue);
+  // gfsr.name("imaginary_error");
+  // Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfsr,imue);
+
+
+  std::stringstream basename;
+  basename << "solvers01_Q" << k;
+  vtkwriter.write(basename.str(),Dune::VTK::appendedraw);
+
+
+
+
+
+
+ // //  // compute  DISCRETIZATION L2 error
+ // //  /* solution */
+ //  if(errornorm == "L2") {
+ //    // define discrete grid function for solution
+ //    typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
+ //    // DGF udgf(gfs, u);
+
+ //    typedef BCAnalytic<PARAM> ES;
+ //    ES es(gv, param);
+ //    U esh(gfs, zero);
+ //    Dune::PDELab::interpolate(es,gfs,esh);
+ //    DGF esdgf(gfs, esh);
+
+ //    typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquared;
+ //    DifferenceSquared difference(esdgf, es);
+
+ //    typename DifferenceSquared::Traits::RangeType l2normsquared(0.0);
+ //    Dune::PDELab::integrateGridFunction(difference,l2normsquared,10);
+
+ //    std::cout<<"L2DiscrError: "<<std::sqrt(std::abs(l2normsquared[0]))<<std::endl;
+ //  }
+
+
+
+ // //  // compute  POLLUTION L2 error
+ // //  /* solution */
+ //  if(errornorm == "L2") {
+ //    // define discrete grid function for solution
+ //    typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
+ //    DGF udgf(gfs, u);
+
+ //    typedef BCAnalytic<PARAM> ES;
+ //    ES es(gv, param);
+ //    U esh(gfs, zero);
+ //    Dune::PDELab::interpolate(es,gfs,esh);
+ //    DGF esdgf(gfs, esh);
+
+
+ //    typedef DifferenceSquaredAdapter<DGF,DGF> DifferenceSquared;
+ //    DifferenceSquared difference(udgf,esdgf);
+
+ //    typename DifferenceSquared::Traits::RangeType l2normsquared(0.0);
+ //    Dune::PDELab::integrateGridFunction(difference,l2normsquared,10);
+
+ //    std::cout<<"L2PolError: "<<std::sqrt(std::abs(l2normsquared[0]))<<std::endl;
+
+ //  }
+
+
+
+
+
+ //  // compute L2 error
+ //  /* solution */
+  if(errornorm == "L2") {
+    // define discrete grid function for solution
+    typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
+    DGF udgf(gfs, u);
+
+    typedef BCAnalytic<PARAM> ES;
+    ES es(gv, param);
+
+    typedef DifferenceSquaredAdapter<DGF,ES> DifferenceSquared;
+    DifferenceSquared difference(udgf,es);
+
+    typename DifferenceSquared::Traits::RangeType l2normsquared(0.0);
+    Dune::PDELab::integrateGridFunction(difference,l2normsquared,10);
+
+    std::cout<<"L2Error: "<<std::sqrt(std::abs(l2normsquared[0]))<<std::endl;
+  }
+
+
+
+  // compute H1 semi norm error
+  /* solution */
+  if(errornorm == "H1") {
+    //discrete function gradient
+    typedef Dune::PDELab::DiscreteGridFunctionGradient<GFS,U> DGFG;
+    DGFG udgfg(gfs, u);
+
+    //gradient of exact solution
+    typedef BCAnalyticGrad<PARAM> ESG;
+    ESG esg(gv, param);
+
+    // difference of gradient
+    typedef DifferenceSquaredAdapter<DGFG,ESG> DifferenceSquaredAdapterg;
+    DifferenceSquaredAdapterg differenceg(udgfg,esg);
+
+    typename DifferenceSquaredAdapterg::Traits::RangeType l2normsquared(0.0);
+    typename DifferenceSquaredAdapter<DGFG,ESG>::Traits::RangeType normsquared(0.0);
+    Dune::PDELab::integrateGridFunction(differenceg,normsquared,10);
+
+    std::cout<<"H1Error: "<<std::sqrt(std::abs(l2normsquared[0]))<<std::endl;
+
+  }
+
+
+
+
+}
+
+#endif

dune/pdelab-systemtesting/helmholtz_stdcomplex/helmholtz_bcanalytic.hh

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+
+#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_BCANALYTIC_HH
+#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_BCANALYTIC_HH
+
+/** \brief A function that defines analytic solution as Dirichlet boundary conditions AND
+    its extension to the interior */
+template<typename T>
+class BCAnalytic
+  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
+									   typename T::Traits::RangeFieldType,
+									   1,
+									   Dune::FieldVector<typename T::Traits::RangeFieldType,1> >,
+					  BCAnalytic<T> > {
+public:
+ typedef Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
+					  typename T::Traits::RangeFieldType,
+					  1,
+					  Dune::FieldVector<typename T::Traits::RangeFieldType,1> > Traits;
+
+  //! construct from grid view
+  BCAnalytic  (const typename Traits::GridViewType& gv_, T& t_) : gv(gv_), t(t_) {}
+
+  //! evaluate extended function on element
+  inline void evaluate (const typename Traits::ElementType& e,
+                        const typename Traits::DomainType& xlocal,
+                        typename Traits::RangeType& y) const
+  {
+    y = t.ua(e, xlocal);
+  }
+
+  //! get a reference to the grid view
+  inline const typename Traits::GridViewType& getGridView () const
+  {
+    return gv;
+  }
+
+private:
+  const typename Traits::GridViewType gv;
+  T& t;
+};
+
+
+/** \brief A function that defines gradient of analytic solution as Dirichlet boundary conditions AND
+    its extension to the interior */
+template<typename T>
+class BCAnalyticGrad
+  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
+									   typename T::Traits::RangeFieldType,
+									   1,
+									   Dune::FieldVector<typename T::Traits::RangeFieldType,1> >,
+					  BCAnalyticGrad<T> > {
+public:
+ typedef Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
+					  typename T::Traits::RangeFieldType,
+					  1,
+					  Dune::FieldVector<typename T::Traits::RangeFieldType,1> > Traits;
+
+  //! construct from grid view
+  BCAnalyticGrad  (const typename Traits::GridViewType& gv_, T& t_) : gv(gv_), t(t_) {}
+
+  const static int dim = T::Traits::GridViewType::Grid::dimension;
+  // typedef typename T::Traits::GridViewType::Grid::ctype ctype;
+
+  //! evaluate extended function on element
+  inline void evaluate (const typename Traits::ElementType& e,
+                        const typename Traits::DomainType& xlocal,
+                        Dune::FieldVector<typename T::Traits::RangeFieldType,dim>& y) const
+  {
+
+    y = t.ugrada(e, xlocal);
+  }
+
+  //! get a reference to the grid view
+  inline const typename Traits::GridViewType& getGridView () const
+  {
+    return gv;
+  }
+
+private:
+  const typename Traits::GridViewType gv;
+  T& t;
+
+};
+
+
+
+
+/*! \brief Adapter returning ||f1(x)-f2(x)||^2 for two given grid functions
+
+  \tparam T1  a grid function type
+  \tparam T2  a grid function type
+*/
+template<typename T1, typename T2>
+class DifferenceSquaredAdapter
+  : public Dune::PDELab::GridFunctionBase<
+  Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType,
+				   typename T1::Traits::RangeFieldType,
+				   1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> >
+  ,DifferenceSquaredAdapter<T1,T2> >
+{
+public:
+  typedef Dune::PDELab::GridFunctionTraits<typename T1::Traits::GridViewType,
+                                           typename T1::Traits::RangeFieldType,
+                                           1,Dune::FieldVector<typename T1::Traits::RangeFieldType,1> > Traits;
+
+  //! constructor
+  DifferenceSquaredAdapter (const T1& t1_, const T2& t2_) : t1(t1_), t2(t2_) {}
+
+  //! \copydoc GridFunctionBase::evaluate()
+  inline void evaluate (const typename Traits::ElementType& e,
+                        const typename Traits::DomainType& x,
+                        typename Traits::RangeType& y) const
+  {
+    typename T1::Traits::RangeType v1,v2;
+    t1.evaluate( e, x, v1);
+    t2.evaluate( e, x, v2);
+    v1 -= v2;
+    y = v1*v1;
+  }
+
+  inline const typename Traits::GridViewType& getGridView () const
+  {
+    return t1.getGridView();
+  }
+
+private:
+  const T1& t1;
+  const T2& t2;
+};
+#endif

dune/pdelab-systemtesting/helmholtz_stdcomplex/helmholtz_bcextension.hh

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+
+#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_BCEXTENSION_HH
+#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_BCEXTENSION_HH
+
+/** \brief A function that defines Dirichlet boundary conditions AND
+    its extension to the interior */
+template<typename T>
+class BCExtension
+  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
+									   typename T::Traits::RangeFieldType,
+									   1,
+									   Dune::FieldVector<typename T::Traits::RangeFieldType,1> >,
+					  BCExtension<T> > {
+public:
+ typedef Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
+					  typename T::Traits::RangeFieldType,
+					  1,
+					  Dune::FieldVector<typename T::Traits::RangeFieldType,1> > Traits;
+
+  //! construct from grid view
+  BCExtension  (const typename Traits::GridViewType& gv_, T& t_) : gv(gv_), t(t_) {}
+
+  //! evaluate extended function on element
+  inline void evaluate (const typename Traits::ElementType& e,
+                        const typename Traits::DomainType& xlocal,
+                        typename Traits::RangeType& y) const
+  {
+    y = t.g(e, xlocal);
+  }
+
+  //! get a reference to the grid view
+  inline const typename Traits::GridViewType& getGridView () const
+  {
+    return gv;
+  }
+
+private:
+  const typename Traits::GridViewType gv;
+  T& t;
+
+
+};
+
+#endif

dune/pdelab-systemtesting/helmholtz_stdcomplex/helmholtz_main.hh

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+
+#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_MAIN_HH
+#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_MAIN_HH
+
+/** \file
+
+    \brief main
+*/
+
+int main(int argc, char** argv)
+{
+  try{
+    //Maybe initialize Mpi
+    Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
+    if(Dune::MPIHelper::isFake)
+      std::cout<< "This is a sequential program." << std::endl;
+    else
+      {
+        if(helper.rank()==0)
+          std::cout << "parallel run on " << helper.size() << " process(es)" << std::endl;
+      }
+
+
+    if (argc!=2)
+      {
+        if (helper.rank()==0)
+          std::cout << "usage: ./error <cfg-file>" << std::endl;
+        return 1;
+      }
+
+
+    // Parse configuration file.
+    std::string config_file(argv[1]);
+    Dune::ParameterTree configuration;
+    Dune::ParameterTreeParser parser;
+
+
+    try{
+      parser.readINITree( config_file, configuration );
+    }
+    catch(...){
+      std::cerr << "Could not read config file \""
+                << config_file << "\"!" << std::endl;
+      exit(1);
+    }
+
+    int level = configuration.get<int>("grid.level");
+    int polynomialdegree = configuration.get<int>("problem_parameter.polynomialdegree");
+    double omega = configuration.get<double>("problem_parameter.omega") ;
+    std::string errornorm = configuration.get<std::string>("norm.errornorm");
+    std::string solver = configuration.get<std::string>("solvers.solver");
+
+
+
+
+    // sequential version
+    if (helper.size()==1)
+      {
+
+        std::vector<double> dof;
+        std::vector<double> tsolve; // number of iterations + solver time
+
+
+        typedef double R;
+        typedef std::complex<R> C;
+
+        Dune::FieldVector<double,2> L(1.);
+        Dune::array<int,2> N(Dune::fill_array<int,2>(1));
+        std::bitset<2> periodic(false);
+        int overlap=0;
+        Dune::YaspGrid<2> grid(L,N,periodic,overlap);
+        typedef Dune::YaspGrid<2>::LeafGridView GV;
+
+        //refine grid
+        grid.globalRefine(level);
+        //get leafGridView
+        const GV& gv=grid.leafGridView();
+
+        //define problem
+        typedef ParametersPlaneWave<GV, C, R> PARAM;
+        //typedef ParametersSphericalWave<GV, C, R> PARAM;
+
+        PARAM param(omega);
+
+
+        if(polynomialdegree == 1) {
+          helmholtz_Qk<1,GV,PARAM>(gv, param,  errornorm, solver);
+        }
+
+        if(polynomialdegree == 2) {
+          helmholtz_Qk<2,GV,PARAM>(gv, param, errornorm, solver);
+        }
+
+
+
+
+      }
+  }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+    return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+    return 1;
+  }
+}
+
+#endif

dune/pdelab-systemtesting/helmholtz_stdcomplex/helmholtz_operator.hh

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+
+
+#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_OPERATOR_HH
+#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_OPERATOR_HH
+
+/** a local operator for solving the equation
+ *
+ *   - \Delta u - \omega^2 u   = f   in \Omega
+ *    \nabla u \cdot n - i \omega u  = 0   on \partial\Omega
+ *
+ * with conforming finite elements on all types of grids in any dimension
+ *
+ */
+template<class PARAM>
+class HelmholtzLocalOperator :
+  public Dune::PDELab::NumericalJacobianApplyVolume<HelmholtzLocalOperator<PARAM> >,
+  public Dune::PDELab::NumericalJacobianVolume<HelmholtzLocalOperator<PARAM> >,
+  public Dune::PDELab::NumericalJacobianApplyBoundary<HelmholtzLocalOperator<PARAM> >,
+  public Dune::PDELab::NumericalJacobianBoundary<HelmholtzLocalOperator<PARAM> >,
+  public Dune::PDELab::FullVolumePattern,
+  public Dune::PDELab::LocalOperatorDefaultFlags
+
+
+{
+public:
+  // pattern assembly flags
+  enum { doPatternVolume = true };
+
+  // residual assembly flags
+  enum { doAlphaVolume = true };
+  enum { doAlphaBoundary = true };                                // assemble boundary
+
+
+
+  HelmholtzLocalOperator (const PARAM& param_, unsigned int intorder_=2)
+     : param(param_ ), intorder(intorder_)
+  {
+  }
+
+
+  // volume integral depending on test and ansatz functions
+  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
+  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
+  {
+    // extract some types
+    typedef typename LFSU::Traits::FiniteElementType::
+      Traits::LocalBasisType::Traits::DomainFieldType DF;
+    typedef typename LFSU::Traits::FiniteElementType::
+      Traits::LocalBasisType::Traits::RangeFieldType RF;
+    typedef typename LFSU::Traits::FiniteElementType::
+      Traits::LocalBasisType::Traits::JacobianType JacobianType;
+    typedef typename LFSU::Traits::FiniteElementType::
+      Traits::LocalBasisType::Traits::RangeType RangeType;
+    typedef typename LFSU::Traits::SizeType size_type;
+
+    // dimensions
+    const int dim = EG::Geometry::dimension;
+
+    // select quadrature rule
+    Dune::GeometryType gt = eg.geometry().type();
+    const Dune::QuadratureRule<DF,dim>&
+      rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
+
+    // loop over quadrature points
+    for (typename Dune::QuadratureRule<DF,dim>::const_iterator
+           it=rule.begin(); it!=rule.end(); ++it)
+      {
+        // evaluate basis functions on reference element
+        std::vector<RangeType> phi(lfsu.size());
+        lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
+
+        // compute u at integration point
+        RF u=0.0;
+        for (size_type i=0; i<lfsu.size(); i++)
+          u = u + x(lfsu,i)*phi[i];
+
+        // evaluate gradient of basis functions on reference element
+        std::vector<JacobianType> js(lfsu.size());
+        lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
+
+        // transform gradients from reference element to real element
+        const typename EG::Geometry::JacobianInverseTransposed
+          jac = eg.geometry().jacobianInverseTransposed(it->position());
+        std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
+        for (size_type i=0; i<lfsu.size(); i++)
+          jac.mv(js[i][0],gradphi[i]);
+
+        // compute gradient of u
+        Dune::FieldVector<RF,dim> gradu(0.0);
+        for (size_type i=0; i<lfsu.size(); i++)
+          gradu.axpy(x(lfsu,i),gradphi[i]);
+
+        // evaluate parameters
+        RF f( param.f(eg, it->position()) );
+
+        // integrate grad u * grad phi_i - omega*omega*u*phi_i - f phi_i
+        RF factor = it->weight()*eg.geometry().integrationElement(it->position());
+        for (size_type i=0; i<lfsu.size(); i++)
+          r.accumulate(lfsu,i,( gradu*gradphi[i] - param.omega*param.omega*u*phi[i] - f*phi[i] )*factor);
+      }
+  }
+
+  // - n \nabla u  = j
+  // n \nabla u − i \omega u = 0
+  // boundary integral
+  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
+  void alpha_boundary (const IG& ig, const LFSU& lfsu_s, const X& x_s,
+                       const LFSV& lfsv_s, R& r_s) const
+  {
+    // assume Galerkin: lfsu_s == lfsv_s
+    // This yields more efficient code since the local functionspace only
+    // needs to be evaluated once, but would be incorrect for a finite volume
+    // method
+
+    // some types
+    typedef typename LFSU::Traits::FiniteElementType::
+      Traits::LocalBasisType::Traits::DomainFieldType DF;
+    typedef typename LFSU::Traits::FiniteElementType::
+      Traits::LocalBasisType::Traits::RangeFieldType RF;
+    typedef typename LFSU::Traits::FiniteElementType::
+      Traits::LocalBasisType::Traits::RangeType Range;
+    typedef typename LFSU::Traits::SizeType size_type;
+
+    // dimensions
+    const int dim = IG::dimension;
+
+
+    Dune::GeometryType gtface = ig.geometryInInside().type();
+
+
+    // select quadrature rule for face
+    const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
+
+    // loop over quadrature points and integrate normal flux
+    for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin();
+         it!=rule.end(); ++it)
+      {
+        // skip rest if we are on Dirichlet boundary
+        if ( param.isDirichlet( ig, it->position() ) )
+          continue;
+
+        // position of quadrature point in local coordinates of element
+        Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
+
+        // evaluate basis functions at integration point
+        std::vector<Range> phi(lfsu_s.size());
+        lfsu_s.finiteElement().localBasis().evaluateFunction(local,phi);
+
+        // evaluate u (e.g. flux may depend on u)
+        RF u=0.0;
+        for (size_type i=0; i<lfsu_s.size(); ++i)
+          u = u + x_s(lfsu_s,i)*phi[i];
+
+
+        if ( param.isNeumann( ig.intersection(), it->position() ) )
+          {
+            // evaluate flux boundary condition
+            RF j = param.j(ig, it->position(), u);
+
+            // integrate j
+            RF factor = it->weight()*ig.geometry().integrationElement(it->position());
+            for (size_type i=0; i<lfsu_s.size(); ++i)
+              r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
+          }
+      }
+  }
+
+
+
+
+
+private:
+
+  const PARAM& param;
+  unsigned int intorder;
+
+};
+
+#endif

dune/pdelab-systemtesting/helmholtz_stdcomplex/parameters_planewave.hh

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+
+#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_PARAMETERS_PLANEWAVE_HH
+#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_PARAMETERS_PLANEWAVE_HH
+
+
+/** \file
+
+    \brief Plane Wave Parameter = Problem description
+*/
+
+/** \brief Defines problem parameter analytical solution
+ */
+template<typename GV, typename RF, typename DF>
+class ParametersPlaneWave  : public Dune::PDELab::DirichletConstraintsParameters
+{
+public:
+
+  typedef Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> > Traits;
+  // typedef DiffusionParameterTraits<GV,RF> Traits;
+
+  ParametersPlaneWave (double omega_)
+     : omega(omega_)
+  {
+  }
+
+
+  //! Dirichlet boundary condition type function
+  template<typename IG>
+  inline bool isDirichlet(const IG& intersection, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal) const
+  {
+    return false;
+  }
+
+  //! Neumann boundary condition type function
+  template<typename IG>
+  inline bool isNeumann(const IG& intersection, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal) const
+  {
+    return true;
+  }
+
+  //! Dirichlet boundary condition value
+  RF  g (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const
+  {
+      return RF(0.,0.);
+  }
+
+  //! Neumann boundary condition value
+  template <typename IG>
+  RF  j (const IG& ig, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal, RF u) const
+  {
+    const int dim = Traits::GridViewType::Grid::dimension;
+    typedef typename Traits::GridViewType::Grid::ctype ctype;
+    Dune::FieldVector<ctype,dim> xglobal = ig.geometry().global(xlocal);
+    RF i(0., 1.);
+    RF x,y;
+    x = xglobal[0];
+    y = xglobal[1];
+    // get the outer normal vector at the face center
+    const Dune::FieldVector<RF, dim> normal = ig.centerUnitOuterNormal();
+    Dune::FieldVector<RF, dim> gradu;
+    gradu[0] = i*omega*cos(theta)*std::exp( i*omega*(x*cos(theta) + y *sin(theta)) );
+    gradu[1] = i*omega*sin(theta)*std::exp( i*omega*(x*cos(theta) + y *sin(theta)) );
+    RF res = gradu * normal;
+    res *= -1;
+    return res;
+  }
+
+
+
+  //! source term
+  template <typename EG>
+  RF f (const EG& eg,  const Dune::FieldVector<typename EG::Geometry::ctype, EG::Geometry::dimension>& position) const
+     {
+       RF res(0.,0.);
+       return res;
+     }
+
+
+  //! exact solution / analytic solution
+  RF  ua (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const
+  {
+
+    const int dim = Traits::GridViewType::Grid::dimension;
+    typedef typename Traits::GridViewType::Grid::ctype ctype;
+    Dune::FieldVector<ctype,dim> xglobal = e.geometry().global(xlocal);
+
+    RF i(0., 1.);
+    RF x,y;
+    x = xglobal[0];
+    y = xglobal[1];
+
+    return std::exp( i*omega*(x*cos(theta) + y *sin(theta)) );
+
+  }
+
+    const static int dim = Traits::GridViewType::Grid::dimension;
+
+    //! exact grad of solution / analytic grad of solution
+   Dune::FieldVector<RF,dim> ugrada (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const
+  {
+
+    typedef typename Traits::GridViewType::Grid::ctype ctype;
+    Dune::FieldVector<ctype,dim> xglobal = e.geometry().global(xlocal);
+    RF i(0., 1.);
+    RF x,y;
+    x = xglobal[0];
+    y = xglobal[1];
+    Dune::FieldVector<RF,dim> res;
+    res[0] = i*omega*cos(theta)*std::exp( i*omega*(x*cos(theta) + y *sin(theta)) );
+    res[1] = i*omega*sin(theta)*std::exp( i*omega*(x*cos(theta) + y *sin(theta)) );
+
+    return res;
+
+  }
+
+
+  double theta = M_PI/2.;
+  double omega;
+
+};
+
+#endif

dune/pdelab-systemtesting/helmholtz_stdcomplex/parameters_sphericalwave.hh

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+
+#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_PARAMETERS_SPHERICALWAVE_HH
+#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_PARAMETERS_SPHERICALWAVE_HH
+
+
+/** \file
+
+    \brief Spherical Wave Parameter = Problem description
+*/
+/** \brief Defines problem parameter: Dirichlet boundary conditions AND its
+ *         extension to the interior, sinks and sources and the reaction rate.
+ */
+template<typename GV, typename RF, typename DF>
+class ParametersSphericalWave  : public Dune::PDELab::DirichletConstraintsParameters
+{
+public:
+
+  typedef Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> > Traits;
+  // typedef DiffusionParameterTraits<GV,RF> Traits;
+
+  ParametersSphericalWave (double omega_)
+     : omega(omega_)
+  {
+  }
+
+
+  //! Dirichlet boundary condition type function
+  template<typename IG>
+  inline bool isDirichlet(const IG& intersection, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal) const
+  {
+    return false;
+  }
+
+  //! Neumann boundary condition type function
+  template<typename IG>
+  inline bool isNeumann(const IG& intersection, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal) const
+  {
+    return true;
+  }
+
+  //! Dirichlet boundary condition value
+  RF  g (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const
+  {
+      return RF(0.,0.);
+  }
+
+  //! Neumann boundary condition value
+  template <typename IG>
+  RF  j (const IG& ig, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal, RF u) const
+  {
+    RF i(0., 1.);
+    return i*omega*u;
+  }
+
+
+
+    //! source term
+    template <typename EG>
+    RF f (const EG& eg,  const Dune::FieldVector<typename EG::Geometry::ctype, EG::Geometry::dimension>& position) const
+    {
+
+
+        Dune::FieldVector<typename EG::Geometry::ctype, EG::Geometry::dimension> xg = eg.geometry().global( position );
+        Dune::FieldVector<DF,EG::Geometry::dimension> midpoint;
+
+        for(int i=0; i<EG::Geometry::dimension; ++i) {
+            midpoint[i] = 0.5;
+        }
+
+        xg -= midpoint;
+
+        RF res(0.,0.);
+
+        double r=0.2;
+
+        double tmp=0.;
+        for(int i=0; i<EG::Geometry::dimension; ++i) {
+            tmp += xg[i]*xg[i];
+        }
+
+        if( std::sqrt( tmp ) < r )
+            res = RF(25,0.);
+
+        return res;
+
+    }
+
+
+    //! exact solution / analytic solution
+    RF  ua (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const
+    {
+
+        // const int dim = Traits::GridViewType::Grid::dimension;
+        // typedef typename Traits::GridViewType::Grid::ctype ctype;
+        // Dune::FieldVector<ctype,dim> xglobal = e.geometry().global(xlocal);
+
+        std::cout<<"Error: Analytical Solution has not been implemented"<<std::endl;
+
+        return RF(0.,0.);
+
+    }
+
+    const static int dim = Traits::GridViewType::Grid::dimension;
+
+    //! exact grad of solution / analytic grad of solution
+    Dune::FieldVector<RF,dim> ugrada (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const
+    {
+
+
+        // typedef typename Traits::GridViewType::Grid::ctype ctype;
+        // Dune::FieldVector<ctype,dim> xglobal = e.geometry().global(xlocal);
+
+        std::cout<<"Error: Analytical Solution has not been implemented"<<std::endl;
+        Dune::FieldVector<RF,dim> res;
+        res[0] = 0.;
+        res[1] = 0.;
+        return res;
+
+    }
+
+
+  double omega;
+
+};
+
+#endif