Add a poisson system test

- 2 Grids
- Polynomial degree 1,2
- dim = 2,3

L2 Error already exported for further use.

dune/pdelab-systemtesting/CMakeLists.txt

 #install headers
-install(FILES pdelab-systemtesting.hh DESTINATION include/dune/pdelab-systemtesting)
-add_subdirectory(helmholtz_stdcomplex)
+#add_subdirectory(helmholtz_stdcomplex)
+
+add_subdirectory(poisson)
\ No newline at end of file

dune/pdelab-systemtesting/poisson/CMakeLists.txt

+add_dune_system_test(BASENAME poisson INIFILE poisson.mini SOURCE poisson_uniform.cc TARGETS output)
+foreach(target ${output})
+  add_dune_ug_flags(${target})
+  add_dune_umfpack_flags(${target})
+  add_dune_mpi_flags(${target})
+  target_link_libraries(${target} ${DUNE_LIBS})
+endforeach()

dune/pdelab-systemtesting/poisson/poisson.mini

+# define some convenience variables
+dim =dim= 2, 3
+grid =grid= yasp, ug
+degree =deg= 1, 2
+naming = {grid}_{dim}d_Q{degree}
+
+# The suffix for the executable
+__exec_suffix = {naming}
+
+outputFile = output_{naming}
+computeError = true
+
+[vtk]
+name = poisson_{naming}
+
+[yaspgrid]
+cells =dim= 1 1, 1 1 1
+extension =dim= 1. 1., 1. 1. 1.
+refinement =dim= 4, 2
+
+[ug]
+cells =dim= 1 1, 1 1 1
+lowerleft =dim= 0. 0., 0. 0. 0.
+upperright =dim= 1. 1., 1. 1. 1.
+refinement =dim= 4, 2
+
+[__STATIC.COMPILE_DEFINITIONS]
+GRIDTYPE =grid= Dune::YaspGrid<{dim}>, Dune::UGGrid<{dim}>
+DEGREE = {degree}
\ No newline at end of file

dune/pdelab-systemtesting/poisson/poisson_uniform.cc

+// -*- tab-width: 4; indent-tabs-mode: nil -*-
+/** \file
+    \brief Solve Poisson problem on various grids (sequential).
+    HANGING_NODES_REFINEMENT is macro used to switch on hanging nodes tests.
+    It is set in "Makefile.am" to generate the executable 'poisson_HN'.
+*/
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include<iostream>
+#include<vector>
+#include<map>
+#include<string>
+#include<dune/common/parallel/mpihelper.hh>
+#include<dune/common/parametertree.hh>
+#include<dune/common/parametertreeparser.hh>
+#include<dune/common/exceptions.hh>
+#include<dune/common/fvector.hh>
+#include<dune/common/float_cmp.hh>
+#include<dune/common/static_assert.hh>
+#include<dune/grid/yaspgrid.hh>
+#include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#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<dune/pdelab/finiteelementmap/pkfem.hh>
+#include<dune/pdelab/finiteelementmap/qkfem.hh>
+#include<dune/pdelab/function/minus.hh>
+#include<dune/pdelab/function/sqr.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
+#include<dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
+#include<dune/pdelab/gridfunctionspace/interpolate.hh>
+#include<dune/pdelab/constraints/common/constraints.hh>
+#include<dune/pdelab/constraints/common/constraintsparameters.hh>
+#include<dune/pdelab/constraints/conforming.hh>
+#include<dune/pdelab/constraints/hangingnode.hh>
+#include<dune/pdelab/common/function.hh>
+#include<dune/pdelab/common/functionutilities.hh>
+#include<dune/pdelab/common/vtkexport.hh>
+#include<dune/pdelab/backend/istlvectorbackend.hh>
+#include<dune/pdelab/backend/istl/bcrsmatrixbackend.hh>
+#include<dune/pdelab/backend/istlmatrixbackend.hh>
+#include<dune/pdelab/backend/istlsolverbackend.hh>
+#include<dune/pdelab/backend/seqistlsolverbackend.hh>
+#include<dune/pdelab/localoperator/laplacedirichletp12d.hh>
+#include<dune/pdelab/localoperator/poisson.hh>
+#include<dune/pdelab/gridoperator/gridoperator.hh>
+#include<dune/pdelab/stationary/linearproblem.hh>
+#include<dune/pdelab/gridfunctionspace/vtk.hh>
+
+#include<dune/testtools/gridconstruction.hh>
+#include<dune/testtools/outputtree.hh>
+
+/*
+  HANGING_NODES_REFINEMENT is macro used to switch on hanging nodes tests.
+  It is set in "Makefile.am" to generate the executable 'poisson_HN'.
+*/
+
+//===============================================================
+//===============================================================
+// Solve the Poisson equation
+//           - \Delta u = f in \Omega,
+//                    u = g on \partial\Omega_D
+//  -\nabla u \cdot \nu = j on \partial\Omega_N
+//===============================================================
+//===============================================================
+
+//===============================================================
+// Define parameter functions f,g,j and \partial\Omega_D/N
+//===============================================================
+
+// function for defining the source term
+template<typename GV, typename RF>
+class F
+  : public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>,
+                                                  F<GV,RF> >
+{
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,F<GV,RF> > BaseT;
+
+  F (const GV& gv) : BaseT(gv) {}
+  inline void evaluateGlobal (const typename Traits::DomainType& x,
+                typename Traits::RangeType& y) const
+  {
+    typename Traits::DomainType center;
+    for (int i=0; i<GV::dimension; i++)
+      center[i] = 0.5;
+    center -= x;
+    typename Traits::RangeType norm2 = center.two_norm2();
+    y = (2*GV::dimension - 4 * norm2) * exp(-1*norm2);
+ }
+};
+
+
+
+// constraints parameter class for selecting boundary condition type
+class BCTypeParam
+  : public Dune::PDELab::DirichletConstraintsParameters /*@\label{bcp:base}@*/
+{
+public:
+
+  template<typename I>
+  bool isDirichlet(
+           const I & intersection,   /*@\label{bcp:name}@*/
+           const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord
+           ) const
+  {
+    return true;
+  }
+
+  template<typename I>
+  bool isNeumann(
+                 const I & intersection,   /*@\label{bcp:name}@*/
+                 const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord
+                 ) const
+  {
+    return !isDirichlet(intersection,coord);
+  }
+
+};
+
+// function for Dirichlet boundary conditions and initialization
+template<typename GV, typename RF>
+class G
+  : public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>,
+                                                  G<GV,RF> >
+{
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,G<GV,RF> > BaseT;
+
+  G (const GV& gv) : BaseT(gv) {}
+  inline void evaluateGlobal (const typename Traits::DomainType& x,
+                typename Traits::RangeType& y) const
+  {
+    typename Traits::DomainType center;
+    for (int i=0; i<GV::dimension; i++)
+      center[i] = 0.5;
+    center -= x;
+    y = exp(-center.two_norm2());
+  }
+};
+
+// function for defining the flux boundary condition
+template<typename GV, typename RF>
+class J
+  : public Dune::PDELab::AnalyticGridFunctionBase<Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1>,
+                                                  J<GV,RF> >
+{
+public:
+  typedef Dune::PDELab::AnalyticGridFunctionTraits<GV,RF,1> Traits;
+  typedef Dune::PDELab::AnalyticGridFunctionBase<Traits,J<GV,RF> > BaseT;
+
+  J (const GV& gv) : BaseT(gv) {}
+  inline void evaluateGlobal (const typename Traits::DomainType& x,
+                typename Traits::RangeType& y) const
+  {
+    y = 0;
+  }
+};
+
+int main(int argc, char** argv)
+{
+  try{
+    //Maybe initialize Mpi
+    Dune::MPIHelper::instance(argc, argv);
+
+    // load the parameter file
+    Dune::ParameterTree params;
+    Dune::ParameterTreeParser::readINITree(argv[1], params);
+
+    typedef GRIDTYPE Grid;
+    static const int dim = Grid::dimension;
+    IniGridFactory<Grid> factory(params);
+    auto grid = factory.getGrid();
+
+    // get view
+    typedef Grid::LeafGridView GV;
+    GV gv = grid->leafGridView();
+
+    const int k=DEGREE;
+    const int intorder=2*k;
+    // make finite element map
+    typedef Grid::ctype DF;
+    typedef Dune::PDELab::QkLocalFiniteElementMap<GV,DF,double,k> FEM;
+    FEM fem(gv);
+
+    BCTypeParam bctype;
+    // solve problem
+    typedef Dune::PDELab::ConformingDirichletConstraints Constraints;
+    Constraints constraints;
+
+    // constants and types
+    typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType R;
+
+    // make grid function space
+    typedef Dune::PDELab::ISTLVectorBackend<> VBE;
+    typedef Dune::PDELab::GridFunctionSpace<GV,FEM,Constraints,VBE> GFS;
+    GFS gfs(gv, fem, constraints);
+    gfs.name("poisson solution");
+
+    // make constraints map and initialize it from a function
+    typedef typename GFS::template ConstraintsContainer<R>::Type C;
+    C cg;
+    cg.clear();
+
+    Dune::PDELab::constraints(bctype, gfs, cg);
+
+    // make grid operator
+    typedef F<GV,R> FType;
+    FType f(gv);
+    typedef J<GV,R> JType;
+    JType j(gv);
+    typedef Dune::PDELab::Poisson<FType,BCTypeParam,JType> LOP;
+    LOP lop(f,bctype,j,intorder);
+
+    typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
+    MBE mbe(45); // Maximal number of nonzeroes per row can be cross-checked by patternStatistics().
+
+    typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,R,R,R,C,C> GO;
+    GO go(gfs,cg,gfs,cg,lop,mbe);
+
+    // make coefficent Vector and initialize it from a function
+    typedef typename GO::Traits::Domain V;
+    V x0(gfs);
+    x0 = 0.0;
+    typedef G<GV,R> GType;
+    GType g(gv);
+    Dune::PDELab::interpolate(g,gfs,x0);
+
+    // Choose ISTL Solver Backend
+    typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LS;
+    LS ls(5000,2);
+
+    typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,V> SLP;
+    SLP slp(go,ls,x0,1e-12);
+    slp.apply();
+
+    // make discrete function object
+    Dune::SubsamplingVTKWriter<GV> vtkwriter( gv, 1 );
+    //Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
+    Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfs,x0);
+    vtkwriter.write(params.get<std::string>("vtk.name"),Dune::VTK::ascii);
+
+    if (params.get<bool>("computeError", false))
+    {
+      // calculate L2 error
+      typedef Dune::PDELab::DiscreteGridFunction<GFS, V> DGF;
+      DGF dgf(gfs, x0);
+      Dune::PDELab::MinusGridFunctionAdapter<DGF, GType> diff(dgf, g);
+      typedef Dune::PDELab::SqrGridFunctionAdapter<Dune::PDELab::MinusGridFunctionAdapter<DGF, GType> > L2Error;
+      L2Error l2(diff);
+      L2Error::Traits::RangeType currentL2Error = static_cast<L2Error::Traits::RangeType>(0.0);
+
+      Dune::PDELab::integrateGridFunction(l2,currentL2Error,intorder);
+      Dune::OutputTree output(params.get<std::string>("outputFile"));
+      output.set("l2error", std::sqrt(currentL2Error));
+    }
+    return 0;
+  }
+  catch (Dune::Exception &e){
+    std::cerr << "Dune reported error: " << e << std::endl;
+  return 1;
+  }
+  catch (...){
+    std::cerr << "Unknown exception thrown!" << std::endl;
+  return 1;
+  }
+}