diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index a56360005533aef4fdadb52ec4da9971cf391c97..3085853afcd434a7f84cee62ae57764510cea67d 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -672,6 +672,8 @@ def define_vector(name):
 
 @preamble
 def define_boundary_lambda(boundary, name):
+    if boundary is None:
+        return "auto {} = [&](const auto& x){{ return 0.0; }};".format(name)
     if boundary.is_global:
         return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr)
     else:
@@ -694,28 +696,53 @@ def define_boundary_function(boundary, name):
                                                                                   )
 
 
+@preamble
+def define_compositegfs_parameterfunction(name, *children):
+    return "Dune::PDELab::CompositeGridFunction<{}> {}({});".format(', '.join('decltype({})'.format(c) for c in children),
+                                                                    name,
+                                                                    ', '.join(children))
+
+
 @symbol
-def name_boundary_function(boundary):
-    define_boundary_function(boundary, "b")
-    # TODO add namespace data here
-    return "b"
+def name_boundary_function(expr):
+    from ufl import MixedElement, VectorElement
+    if isinstance(expr, VectorElement):
+        element, (_, boundary) = get_constraints_metadata(expr)
+        from dune.perftool.generation import get_global_context_value
+        name = get_global_context_value("namedata").get(id(boundary), "zero")
+        define_boundary_function(boundary, name)
+        return name
+    if isinstance(expr, MixedElement):
+        children = tuple(name_boundary_function(el) for el in expr.sub_elements())
+        name = '_'.join(children)
+        define_compositegfs_parameterfunction(name, *children)
+        return name
+
+    # This is a scalar leaf of the ansatz function tree
+    element, (_, boundary) = get_constraints_metadata(expr)
+
+    from dune.perftool.generation import get_global_context_value
+    name = get_global_context_value("namedata").get(id(boundary), "zero")
+
+    define_boundary_function(boundary, name)
+    return name
 
 
 @preamble
 def interpolate_vector(name):
     define_vector(name)
-    element, (_, boundary) = get_constraints_metadata(_form.coefficients()[0].ufl_element())
-    b = name_boundary_function(boundary)
+    element = _form.coefficients()[0].ufl_element()
+    bf = name_boundary_function(element)
     gfs = name_gfs(element)
-    return "Dune::PDELab::interpolate({}, {}, {});".format(b,
+    return "Dune::PDELab::interpolate({}, {}, {});".format(bf,
                                                            gfs,
                                                            name,
                                                            )
 
 
 def maybe_interpolate_vector(name):
-    element, (_, boundary) = get_constraints_metadata(_form.coefficients()[0].ufl_element())
-    if boundary:
+    element = _form.coefficients()[0].ufl_element()
+    if has_constraints(element):
         interpolate_vector(name)
     else:
         define_vector(name)
diff --git a/test/stokes/CMakeLists.txt b/test/stokes/CMakeLists.txt
index 1e8651d27a1e5eed3ce487caec5c081377ea1178..0160ab746066dd84113b5b9a2bfcc920d8d5b7ed 100644
--- a/test/stokes/CMakeLists.txt
+++ b/test/stokes/CMakeLists.txt
@@ -1,3 +1,5 @@
+dune_symlink_to_source_files(FILES hagenpoiseuille_ref.vtu)
+
 add_generated_executable(UFLFILE stokes.ufl
                          TARGET stokes_numdiff
                          FORM_COMPILER_ARGS --numerical-jacobian
@@ -5,4 +7,5 @@ add_generated_executable(UFLFILE stokes.ufl
 
 dune_add_system_test(TARGET stokes_numdiff
                      INIFILE stokes_numdiff.mini
+                     SCRIPT dune_vtkcompare.py
                      )
diff --git a/test/stokes/hagenpoiseuille_ref.vtu b/test/stokes/hagenpoiseuille_ref.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..f6d851aca8bc95b570f13d9b26b9de691f129460
--- /dev/null
+++ b/test/stokes/hagenpoiseuille_ref.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1d887eb3a33f5ad367149e61df7e7f82092d74b6c46b3d649c848775e2b3cbfa
+size 80161
diff --git a/test/stokes/stokes.ufl b/test/stokes/stokes.ufl
index e84868a14111d0bd4e3b587ecdc97a32577de9b7..9252834820f67c72a6b28f08b6b1bd1b164ad04c 100644
--- a/test/stokes/stokes.ufl
+++ b/test/stokes/stokes.ufl
@@ -1,4 +1,4 @@
-v_bctype = Expression("if (x[0] < 1e-8) return 1; else return 0;", on_intersection=True)
+v_bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
 v_dirichlet = Expression("Dune::FieldVector<double, 2> y(0.0); y[0] = 4*x[1]*(1.-x[1]); return y;")
 
 cell = triangle
diff --git a/test/stokes/stokes_numdiff.mini b/test/stokes/stokes_numdiff.mini
index 671667ea46b5a886582c402af34aaf71880c3725..0f52d1358312990c81b151fb4e39fc295927f2ec 100644
--- a/test/stokes/stokes_numdiff.mini
+++ b/test/stokes/stokes_numdiff.mini
@@ -2,6 +2,11 @@ __name = stokes_numdiff
 
 lowerleft = 0.0 0.0
 upperright = 1.0 1.0
-elements = 4 4
+elements = 16 16
 elementType = simplical
-printmatrix = true
+printmatrix = false
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = hagenpoiseuille_ref
+extension = vtu