diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 91d0e76d33ef6b43f1514d5cbf01590502421ff4..46d54a8db7d4d8fa046ab3676f189750a2bd12fb 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -774,14 +774,17 @@ def name_vtkfile():
 
 @preamble
 def print_matrix():
+    ini = name_initree()
     t_go = type_gridoperator()
     n_go = name_gridoperator()
     v = name_vector()
-    return ["typedef typename {}::Traits::Jacobian M;".format(t_go),
-            "M m({});".format(n_go),
-            "{}.jacobian({},m);".format(n_go, v),
-            "using Dune::PDELab::Backend::native;",
-            "Dune::printmatrix(std::cout, native(m),\"global stiffness matrix\",\"row\",9,1);"]
+    return ["if ({}.get<bool>(\"printmatrix\", false)) {{".format(ini),
+            "  typedef typename {}::Traits::Jacobian M;".format(t_go),
+            "  M m({});".format(n_go),
+            "  {}.jacobian({},m);".format(n_go, v),
+            "  using Dune::PDELab::Backend::native;",
+            "  Dune::printmatrix(std::cout, native(m),\"global stiffness matrix\",\"row\",9,1);",
+            "}"]
 
 
 @preamble
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index abffb46dda7b65578411d0af5b6b83457d7a0b47..d8ac5fd78237fce7ac10d30a2ab2bea89540a521 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -1,6 +1,7 @@
 """ Generators for parameter functions """
 
-from dune.perftool.generation import (class_basename,
+from dune.perftool.generation import (cached,
+                                      class_basename,
                                       class_member,
                                       constructor_parameter,
                                       generator_factory,
@@ -105,12 +106,14 @@ def define_parameter_temporary(t):
     return _define
 
 
+@cached
 def cell_parameter_function(name, expr, restriction, t='double'):
     temporary_variable(name, shape=(), decl_method=define_parameter_temporary(t))
     define_parameter_function_class_member(name, expr, t, True)
     evaluate_cell_parameter_function(name, restriction)
 
 
+@cached
 def intersection_parameter_function(name, expr, t='double'):
     temporary_variable(name, shape=(), decl_method=define_parameter_temporary(t))
     define_parameter_function_class_member(name, expr, t, False)
diff --git a/test/laplace/CMakeLists.txt b/test/laplace/CMakeLists.txt
index d430d22135fdc9f0321b0462fecd1803eac2ad7e..d58b9c54d8757c6b71e2c714e3af23f8ae329135 100644
--- a/test/laplace/CMakeLists.txt
+++ b/test/laplace/CMakeLists.txt
@@ -30,4 +30,5 @@ add_generated_executable(UFLFILE laplace_dg.ufl
 dune_add_system_test(TARGET laplace_dg_symdiff
                      INIFILE laplace_dg_symdiff.mini)
 
-add_executable(dgref reference_main.cc)
+add_executable(laplace_dg_ref reference_main.cc)
+set_target_properties(laplace_dg_ref PROPERTIES EXCLUDE_FROM_ALL 1)
diff --git a/test/laplace/laplace_dg_numdiff.mini b/test/laplace/laplace_dg_numdiff.mini
index 269235ec41d246875cd534bb3db4d2d70a6047fc..ba1e0014857fc505f0f0f22e6ea8e0b40a1fb095 100644
--- a/test/laplace/laplace_dg_numdiff.mini
+++ b/test/laplace/laplace_dg_numdiff.mini
@@ -4,3 +4,4 @@ lowerleft = 0.0 0.0
 upperright = 1.0 1.0
 elements = 2 2
 elementType = simplical
+printmatrix = true
diff --git a/test/laplace/laplace_dg_symdiff.mini b/test/laplace/laplace_dg_symdiff.mini
index 6544b4d1192e3bed7ea022576e8ccb841eda7ef4..b03b5ada766d0c520e1afe2ee919a713d80e47ca 100644
--- a/test/laplace/laplace_dg_symdiff.mini
+++ b/test/laplace/laplace_dg_symdiff.mini
@@ -4,3 +4,5 @@ lowerleft = 0.0 0.0
 upperright = 1.0 1.0
 elements = 2 2
 elementType = simplical
+printmatrix = true
+
diff --git a/test/laplace/laplace_numdiff.mini b/test/laplace/laplace_numdiff.mini
index 60d37cd0e416deb3c639dccafe43be3233be39dd..8ee48e24d62dc2a16b875a27a02b0bca7b6bba60 100644
--- a/test/laplace/laplace_numdiff.mini
+++ b/test/laplace/laplace_numdiff.mini
@@ -4,3 +4,4 @@ lowerleft = 0.0 0.0
 upperright = 1.0 1.0
 elements = 4 4
 elementType = simplical
+printmatrix = true
diff --git a/test/laplace/laplace_symdiff.mini b/test/laplace/laplace_symdiff.mini
index ec27540f8ee06a2b7b19519c9272088bab51185d..c0ee0dd33c0a98aeb36cb68199946593ee752d41 100644
--- a/test/laplace/laplace_symdiff.mini
+++ b/test/laplace/laplace_symdiff.mini
@@ -4,3 +4,4 @@ lowerleft = 0.0 0.0
 upperright = 1.0 1.0
 elements = 4 4
 elementType = simplical
+printmatrix = true
diff --git a/test/laplace/reference_driver.hh b/test/laplace/reference_driver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..db1b55946833b34ae8e925737c7c127241b57f2a
--- /dev/null
+++ b/test/laplace/reference_driver.hh
@@ -0,0 +1,144 @@
+#ifndef _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
+#define _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
+
+#include <dune/pdelab/gridoperator/gridoperator.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/gridfunctionspace/vtk.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <string>
+#include <dune/common/parametertree.hh>
+#include <dune/pdelab/finiteelementmap/opbfem.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include <dune/pdelab/stationary/linearproblem.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/testtools/gridconstruction.hh>
+#include <dune/pdelab/localoperator/convectiondiffusiondg.hh>
+#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
+
+template<typename GV, typename RF>
+class CDProb
+{
+  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
+
+  public:
+  typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
+
+  //! tensor diffusion coefficient
+  typename Traits::PermTensorType
+  A (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    typename Traits::PermTensorType I;
+    for (std::size_t i=0; i<Traits::dimDomain; i++)
+      for (std::size_t j=0; j<Traits::dimDomain; j++)
+        I[i][j] = (i==j) ? 1 : 0;
+    return I;
+  }
+
+  //! velocity field
+  typename Traits::RangeType
+  b (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    typename Traits::RangeType v(0.0);
+    return v;
+  }
+
+  //! sink term
+  typename Traits::RangeFieldType
+  c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! source term
+  typename Traits::RangeFieldType
+  f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! boundary condition type function
+  BCType
+  bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
+  }
+
+  //! Dirichlet boundary condition value
+  typename Traits::RangeFieldType
+  g (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! Neumann boundary condition
+  typename Traits::RangeFieldType
+  j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! outflow boundary condition
+  typename Traits::RangeFieldType
+  o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return 0.0;
+  }
+};
+
+
+void driver(int argc, char** argv){  typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none, 1> VectorBackend;
+  static const int dim = 2;
+  typedef Dune::UGGrid<dim> Grid;
+  typedef Grid::LeafGridView GV;
+  typedef Grid::ctype DF;
+  typedef double R;
+  typedef Dune::PDELab::OPBLocalFiniteElementMap<DF, R, 1, dim, Dune::GeometryType::simplex> DG1_FEM;
+  typedef Dune::PDELab::NoConstraints NoConstraintsAssembler;
+  typedef Dune::PDELab::GridFunctionSpace<GV, DG1_FEM, NoConstraintsAssembler, VectorBackend> DG1_DIRICHLET_GFS;
+  Dune::ParameterTree initree;
+  Dune::ParameterTreeParser::readINITree(argv[1], initree);
+  IniGridFactory<Grid> factory(initree);
+  std::shared_ptr<Grid> grid = factory.getGrid();
+  GV gv = grid->leafGridView();
+  DG1_FEM dg1_fem;
+  DG1_DIRICHLET_GFS dg1_dirichlet_gfs(gv, dg1_fem);
+  dg1_dirichlet_gfs.name("bla");
+  typedef Dune::SubsamplingVTKWriter<GV> VTKWriter;
+  int sublevel = initree.get<int>("vtk.subsamplinglevel", 0);
+  VTKWriter vtkwriter(gv, sublevel);
+  using LOP = Dune::PDELab::ConvectionDiffusionDG<CDProb<GV, R>, DG1_FEM>;
+  typedef DG1_DIRICHLET_GFS::ConstraintsContainer<R>::Type DG1_CC;
+  typedef Dune::PDELab::istl::BCRSMatrixBackend<> MatrixBackend;
+  typedef Dune::PDELab::GridOperator<DG1_DIRICHLET_GFS, DG1_DIRICHLET_GFS, LOP, MatrixBackend, DF, R, R, DG1_CC, DG1_CC> GO;
+  typedef GO::Traits::Domain V;
+  V x(dg1_dirichlet_gfs);
+  x = 0.0;
+  std::string vtkfile = initree.get<std::string>("wrapper.vtkcompare.name", "output");
+  typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LinearSolver;
+  typedef Dune::PDELab::StationaryLinearProblemSolver<GO, LinearSolver, V> SLP;
+  DG1_CC dg1_cc;
+  dg1_cc.clear();
+  CDProb<GV, R> params;
+  LOP lop(params, Dune::PDELab::ConvectionDiffusionDGMethod::SIPG);
+  int generic_dof_estimate =  6 * dg1_dirichlet_gfs.maxLocalSize();
+  int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate);
+  MatrixBackend mb(dofestimate);
+  GO go(dg1_dirichlet_gfs, dg1_cc, dg1_dirichlet_gfs, dg1_cc, lop, mb);
+  std::cout << "gfs with " << dg1_dirichlet_gfs.size() << " dofs generated  "<< std::endl;
+  std::cout << "cc with " << dg1_cc.size() << " dofs generated  "<< std::endl;
+  LinearSolver ls(false);
+  double reduction = initree.get<double>("reduction", 1e-12);
+  SLP slp(go, ls, x, reduction);
+  slp.apply();
+  typedef typename GO::Traits::Jacobian M;
+  M m(go);
+  go.jacobian(x,m);
+  using Dune::PDELab::Backend::native;
+  Dune::printmatrix(std::cout, native(m),"global stiffness matrix","row",9,1);
+  Dune::PDELab::addSolutionToVTKWriter(vtkwriter, dg1_dirichlet_gfs, x);
+  vtkwriter.write(vtkfile, Dune::VTK::ascii);
+}
+
+#endif //_HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
diff --git a/test/laplace/reference_main.cc b/test/laplace/reference_main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ff00d5c060de6a23fafc79ffb7aeabb8f5d1ac9d
--- /dev/null
+++ b/test/laplace/reference_main.cc
@@ -0,0 +1,33 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+
+#include"/home/dominic/dune/dune-perftool/test/laplace/reference_driver.hh"
+
+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
+      std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
+        <<" processes!"<<std::endl;
+
+    driver(argc, argv);
+
+    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;
+  }
+}
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index d18236b90f5e2a6e36959f4e152f9de6bdad0225..e00a32f427770b64748a866bfcee13bab7171cf4 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -1,3 +1,9 @@
+# Symlink any reference solutions
+dune_symlink_to_source_files(FILES poisson_ref.vtu
+                                   poisson_dg_ref.vtu
+                             )
+
+
 # 1. Poisson Test Case: source f, bc g + numerical differentiation
 add_generated_executable(UFLFILE poisson.ufl
                          TARGET poisson_numdiff
@@ -9,9 +15,6 @@ dune_add_system_test(TARGET poisson_numdiff
                      SCRIPT dune_vtkcompare.py
                      )
 
-# Symlink any reference solutions
-dune_symlink_to_source_files(FILES poisson_ref.vtu)
-
 
 # 2. Poisson Test Case: source f, bc g + symbolic differentiation
 add_generated_executable(UFLFILE poisson.ufl
@@ -42,3 +45,39 @@ add_generated_executable(UFLFILE poisson_neumann.ufl
 dune_add_system_test(TARGET poisson_neumann_symdiff
                      INIFILE poisson_neumann_symdiff.mini
                      )
+
+# 5. Poisson Test Case: DG, f + pure dirichlet, numeric differentiation
+add_generated_executable(UFLFILE poisson_dg.ufl
+                         TARGET poisson_dg_numdiff
+                         FORM_COMPILER_ARGS --numerical-jacobian
+                         )
+
+dune_add_system_test(TARGET poisson_dg_numdiff
+                     INIFILE poisson_dg_numdiff.mini
+                     SCRIPT dune_vtkcompare.py
+                     )
+
+# 6. Poisson Test Case: DG, f + pure dirichlet, symbolic differentiation
+add_generated_executable(UFLFILE poisson_dg.ufl
+                         TARGET poisson_dg_symdiff
+                         )
+
+dune_add_system_test(TARGET poisson_dg_symdiff
+                     INIFILE poisson_dg_symdiff.mini
+                     SCRIPT dune_vtkcompare.py
+                     )
+
+# 7. Poisson Test Case: DG, mixed bc, numeric differentiation
+add_generated_executable(UFLFILE poisson_dg_neumann.ufl
+                         TARGET poisson_dg_neumann_numdiff
+                         FORM_COMPILER_ARGS --numerical-jacobian
+                         )
+
+dune_add_system_test(TARGET poisson_dg_neumann_numdiff
+                     INIFILE poisson_dg_neumann_numdiff.mini
+                     )
+
+# Add the reference code with the PDELab localoperator that produced
+# the reference vtk file
+add_executable(poisson_dg_ref reference_main.cc)
+set_target_properties(poisson_dg_ref PROPERTIES EXCLUDE_FROM_ALL 1)
diff --git a/test/poisson/poisson_dg.ufl b/test/poisson/poisson_dg.ufl
index 416ccf8fcb50430bdd126cda2408261a5752572e..0ee5a857931d2ceea62b8cf5e52c5aed2b696851 100644
--- a/test/poisson/poisson_dg.ufl
+++ b/test/poisson/poisson_dg.ufl
@@ -1,7 +1,5 @@
 f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
-j = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; double s; if (x[1]>0.5) s=1.; else s=-1.; return -2.*s*(x[1]-0.5)*std::exp(-1.*c.two_norm2());", on_intersection=True)
-bctype = Expression("if ((x[1]<1e-8) || (x[1]>1.-1e-8)) return 1; else return 0;", on_intersection=True)
+g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", on_intersection=True)
 
 cell = triangle
 V = FiniteElement("DG", cell, 1)
@@ -9,20 +7,20 @@ V = FiniteElement("DG", cell, 1)
 u = TrialFunction(V)
 v = TestFunction(V)
 
-n = FacetNormal(cell)
+n = FacetNormal(cell)('+')
 
 gamma = 1.0
 theta = 1.0
 
-ds = ds(subdomain_data=bctype)
-
 r = inner(grad(u), grad(v))*dx \
-  - inner(n, avg(grad(u)))*jump(v)*(dS+ds(0)) \
-  + gamma*jump(u)*jump(v)*(dS+ds(0)) \
-  - theta*jump(u)*inner(grad(v), n)*(dS+ds(0)) \
+  + inner(n, avg(grad(u)))*jump(v)*dS \
+  + gamma*jump(u)*jump(v)*dS \
+  + theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - inner(n, grad(u))*v*ds \
+  + gamma*u*v*ds \
+  - theta*u*inner(grad(v), n)*ds \
   - f*v*dx \
-  + theta*g*inner(grad(v), n)*ds(0) \
-  - gamma*g*v*ds(0) \
-  + j*v*ds(1)
+  + theta*g*inner(grad(v), n)*ds \
+  - gamma*g*v*ds
 
 forms = [r]
diff --git a/test/poisson/poisson_dg_numdiff.mini b/test/poisson/poisson_dg_numdiff.mini
new file mode 100644
index 0000000000000000000000000000000000000000..146955e8a8d3699cfcb0e7a145a1619a1c49ab60
--- /dev/null
+++ b/test/poisson/poisson_dg_numdiff.mini
@@ -0,0 +1,11 @@
+__name = poisson_dg_numdiff
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_dg_ref
+extension = vtu
diff --git a/test/poisson/poisson_dg_ref.vtu b/test/poisson/poisson_dg_ref.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..06396d8a6ecf6158bd7045133fe19c3e4d56d389
--- /dev/null
+++ b/test/poisson/poisson_dg_ref.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:34ba9a3065699851e3bdc68c2ef14f525c0bdfb42fc46c07c4218109e2d30c97
+size 226320
diff --git a/test/poisson/reference_driver.hh b/test/poisson/reference_driver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6c3e45f2f2c59e8c7a08554aa288e86d6dac3cba
--- /dev/null
+++ b/test/poisson/reference_driver.hh
@@ -0,0 +1,145 @@
+#ifndef _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
+#define _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
+
+#include <dune/pdelab/gridoperator/gridoperator.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/gridfunctionspace/vtk.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <string>
+#include <dune/common/parametertree.hh>
+#include <dune/pdelab/finiteelementmap/opbfem.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include <dune/pdelab/stationary/linearproblem.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/testtools/gridconstruction.hh>
+#include <dune/pdelab/localoperator/convectiondiffusiondg.hh>
+#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
+
+template<typename GV, typename RF>
+class CDProb
+{
+  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
+
+  public:
+  typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
+
+  //! tensor diffusion coefficient
+  typename Traits::PermTensorType
+  A (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    typename Traits::PermTensorType I;
+    for (std::size_t i=0; i<Traits::dimDomain; i++)
+      for (std::size_t j=0; j<Traits::dimDomain; j++)
+        I[i][j] = (i==j) ? 1 : 0;
+    return I;
+  }
+
+  //! velocity field
+  typename Traits::RangeType
+  b (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    typename Traits::RangeType v(0.0);
+    return v;
+  }
+
+  //! sink term
+  typename Traits::RangeFieldType
+  c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! source term
+  typename Traits::RangeFieldType
+  f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    auto global = e.geometry().global(x);
+    Dune::FieldVector<double,2> c(0.5);
+    c -= global;
+    return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
+  }
+
+  //! boundary condition type function
+  BCType
+  bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
+  }
+
+  //! Dirichlet boundary condition value
+  typename Traits::RangeFieldType
+  g (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+  {
+    auto global = e.geometry().global(x);
+    Dune::FieldVector<double,2> c(0.5);
+    c-= global;
+    return std::exp(-1.*c.two_norm2());
+  }
+
+  //! Neumann boundary condition
+  typename Traits::RangeFieldType
+  j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return 0.0;
+  }
+
+  //! outflow boundary condition
+  typename Traits::RangeFieldType
+  o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+  {
+    return 0.0;
+  }
+};
+
+
+void driver(int argc, char** argv){  typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none, 1> VectorBackend;
+  static const int dim = 2;
+  typedef Dune::UGGrid<dim> Grid;
+  typedef Grid::LeafGridView GV;
+  typedef Grid::ctype DF;
+  typedef double R;
+  typedef Dune::PDELab::OPBLocalFiniteElementMap<DF, R, 1, dim, Dune::GeometryType::simplex> DG1_FEM;
+  typedef Dune::PDELab::NoConstraints NoConstraintsAssembler;
+  typedef Dune::PDELab::GridFunctionSpace<GV, DG1_FEM, NoConstraintsAssembler, VectorBackend> DG1_DIRICHLET_GFS;
+  Dune::ParameterTree initree;
+  Dune::ParameterTreeParser::readINITree(argv[1], initree);
+  IniGridFactory<Grid> factory(initree);
+  std::shared_ptr<Grid> grid = factory.getGrid();
+  GV gv = grid->leafGridView();
+  DG1_FEM dg1_fem;
+  DG1_DIRICHLET_GFS dg1_dirichlet_gfs(gv, dg1_fem);
+  dg1_dirichlet_gfs.name("bla");
+  typedef Dune::SubsamplingVTKWriter<GV> VTKWriter;
+  int sublevel = initree.get<int>("vtk.subsamplinglevel", 0);
+  VTKWriter vtkwriter(gv, sublevel);
+  using LOP = Dune::PDELab::ConvectionDiffusionDG<CDProb<GV, R>, DG1_FEM>;
+  typedef DG1_DIRICHLET_GFS::ConstraintsContainer<R>::Type DG1_CC;
+  typedef Dune::PDELab::istl::BCRSMatrixBackend<> MatrixBackend;
+  typedef Dune::PDELab::GridOperator<DG1_DIRICHLET_GFS, DG1_DIRICHLET_GFS, LOP, MatrixBackend, DF, R, R, DG1_CC, DG1_CC> GO;
+  typedef GO::Traits::Domain V;
+  V x(dg1_dirichlet_gfs);
+  x = 0.0;
+  std::string vtkfile = initree.get<std::string>("wrapper.vtkcompare.name", "output");
+  typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LinearSolver;
+  typedef Dune::PDELab::StationaryLinearProblemSolver<GO, LinearSolver, V> SLP;
+  DG1_CC dg1_cc;
+  dg1_cc.clear();
+  CDProb<GV, R> params;
+  LOP lop(params, Dune::PDELab::ConvectionDiffusionDGMethod::SIPG);
+  int generic_dof_estimate =  6 * dg1_dirichlet_gfs.maxLocalSize();
+  int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate);
+  MatrixBackend mb(dofestimate);
+  GO go(dg1_dirichlet_gfs, dg1_cc, dg1_dirichlet_gfs, dg1_cc, lop, mb);
+  std::cout << "gfs with " << dg1_dirichlet_gfs.size() << " dofs generated  "<< std::endl;
+  std::cout << "cc with " << dg1_cc.size() << " dofs generated  "<< std::endl;
+  LinearSolver ls(false);
+  double reduction = initree.get<double>("reduction", 1e-12);
+  SLP slp(go, ls, x, reduction);
+  slp.apply();
+  Dune::PDELab::addSolutionToVTKWriter(vtkwriter, dg1_dirichlet_gfs, x);
+  vtkwriter.write(vtkfile, Dune::VTK::ascii);
+}
+
+#endif //_HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
diff --git a/test/poisson/reference_main.cc b/test/poisson/reference_main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4304d9b0c61a7be9979993a5a0d24458d0873f99
--- /dev/null
+++ b/test/poisson/reference_main.cc
@@ -0,0 +1,33 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+
+#include "reference_driver.hh"
+
+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
+      std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
+        <<" processes!"<<std::endl;
+
+    driver(argc, argv);
+
+    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;
+  }
+}