diff --git a/python/dune/perftool/blockstructured/basis.py b/python/dune/perftool/blockstructured/basis.py
index 23d89469f02536d0345515b2cf84c7e11bac5ff1..9ca4eb08c0c2607b6ba9ecd263e532c80b7ac093 100644
--- a/python/dune/perftool/blockstructured/basis.py
+++ b/python/dune/perftool/blockstructured/basis.py
@@ -17,7 +17,7 @@ from dune.perftool.pdelab.driver import (isPk,
 from dune.perftool.pdelab.driver.gridfunctionspace import basetype_range
 from dune.perftool.pdelab.geometry import world_dimension
 from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell
-from dune.perftool.pdelab.spaces import type_gfs
+from dune.perftool.pdelab.spaces import type_leaf_gfs
 from dune.perftool.pdelab.restriction import restricted_name
 from dune.perftool.blockstructured.spaces import lfs_inames
 
@@ -28,7 +28,7 @@ import pymbolic.primitives as prim
 @backend(interface="typedef_localbasis", name="blockstructured")
 @class_member(classtag="operator")
 def typedef_localbasis(element, name):
-    df = "typename {}::Traits::GridView::ctype".format(type_gfs(element))
+    df = "typename {}::Traits::GridView::ctype".format(type_leaf_gfs(element))
     r = basetype_range()
     dim = world_dimension()
     if isPk(element):
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 9396782c60ede97a6173d1fa1b8d5e8f33f9e233..5ae88fa6b5ff5adc61c8ef2013b3df6bdca94276 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -90,16 +90,16 @@ class PDELabInterface(object):
         return pymbolic_reference_gradient(element, restriction, number, context)
 
     def pymbolic_trialfunction_gradient(self, element, restriction, index):
-        return pymbolic_trialfunction_gradient(element, restriction, index)
+        return pymbolic_trialfunction_gradient(self.visitor, element, restriction, index)
 
     def pymbolic_apply_function_gradient(self, element, restriction, index):
-        return pymbolic_apply_function_gradient(element, restriction, index)
+        return pymbolic_apply_function_gradient(self.visitor, element, restriction, index)
 
     def pymbolic_trialfunction(self, element, restriction, index):
-        return pymbolic_trialfunction(element, restriction, index)
+        return pymbolic_trialfunction(self.visitor, element, restriction, index)
 
     def pymbolic_apply_function(self, element, restriction, index):
-        return pymbolic_apply_function(element, restriction, index)
+        return pymbolic_apply_function(self.visitor, element, restriction, index)
 
     #
     # Parameter function related generator functions
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 5adc0eb48f2ba72a2c29b59ee7a66a046ed78edc..5b745af247fd202d9f261ef6873fbf1d7932e923 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -90,35 +90,35 @@ def accumulation_mangler(target, func, dtypes):
                                   )
 
 
-def pymbolic_trialfunction_gradient(element, restriction, index):
+def pymbolic_trialfunction_gradient(visitor, element, restriction, index):
     rawname = "gradu_{}".format(index)
     name = restricted_name(rawname, restriction)
     container = name_coefficientcontainer(restriction)
-    evaluate_coefficient_gradient(element, name, container, restriction, index)
+    evaluate_coefficient_gradient(visitor, element, name, container, restriction, index)
     return Variable(name)
 
 
-def pymbolic_trialfunction(element, restriction, index):
+def pymbolic_trialfunction(visitor, element, restriction, index):
     rawname = "u_{}".format(index)
     name = restricted_name(rawname, restriction)
     container = name_coefficientcontainer(restriction)
-    evaluate_coefficient(element, name, container, restriction, index)
+    evaluate_coefficient(visitor, element, name, container, restriction, index)
     return Variable(name)
 
 
-def pymbolic_apply_function_gradient(element, restriction, index):
+def pymbolic_apply_function_gradient(visitor, element, restriction, index):
     rawname = "gradz_func_{}".format(index)
     name = restricted_name(rawname, restriction)
     container = name_applycontainer(restriction)
-    evaluate_coefficient_gradient(element, name, container, restriction, index)
+    evaluate_coefficient_gradient(visitor, element, name, container, restriction, index)
     return Variable(name)
 
 
-def pymbolic_apply_function(element, restriction, index):
+def pymbolic_apply_function(visitor, element, restriction, index):
     rawname = "z_func_{}".format(index)
     name = restricted_name(rawname, restriction)
     container = name_applycontainer(restriction)
-    evaluate_coefficient(element, name, container, restriction, index)
+    evaluate_coefficient(visitor, element, name, container, restriction, index)
     return Variable(name)
 
 
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index 4dc813e7e7ee509b2c8858ca9956083a65dbe3ca..af3c78132b6a8cb6f6a3e20f19a7396517d767ed 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -154,7 +154,7 @@ def shape_as_pymbolic(shape):
 
 
 @kernel_cached
-def evaluate_coefficient(element, name, container, restriction, index):
+def evaluate_coefficient(visitor, element, name, container, restriction, index):
     sub_element = element
     if element.num_sub_elements() > 0:
         sub_element = element.extract_component(index)[1]
@@ -165,7 +165,7 @@ def evaluate_coefficient(element, name, container, restriction, index):
     temporary_variable(name, shape=(),)
 
     lfs = name_lfs(element, restriction, index)
-    basis = pymbolic_basis(sub_element, restriction, 0, context='trial')
+    basis = visitor.interface.pymbolic_basis(sub_element, restriction, 0, context='trial')
     basisindex, = get_pymbolic_indices(basis)
 
 
@@ -188,7 +188,7 @@ def evaluate_coefficient(element, name, container, restriction, index):
 
 
 @kernel_cached
-def evaluate_coefficient_gradient(element, name, container, restriction, index):
+def evaluate_coefficient_gradient(visitor, element, name, container, restriction, index):
     sub_element = element
     if element.num_sub_elements() > 0:
         sub_element = element.extract_component(index)[1]
@@ -203,7 +203,7 @@ def evaluate_coefficient_gradient(element, name, container, restriction, index):
     dimindex = component_iname(count=0)
 
     lfs = name_lfs(element, restriction, index)
-    basis = pymbolic_reference_gradient(sub_element, restriction, 0, context='trialgrad')
+    basis = visitor.interface.pymbolic_reference_gradient(sub_element, restriction, 0, context='trialgrad')
     basisindex, _ = get_pymbolic_indices(basis)
     from dune.perftool.tools import maybe_wrap_subscript
     basis = maybe_wrap_subscript(basis, Variable(dimindex))
diff --git a/python/dune/perftool/pdelab/driver/gridfunctionspace.py b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
index 9248d8f8a74f0d2040d79ed7ee0c0a47905a7b4b..8509b30fd99f0fdb2de9a29e504d88975b74712b 100644
--- a/python/dune/perftool/pdelab/driver/gridfunctionspace.py
+++ b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
@@ -124,6 +124,9 @@ def typedef_fem(element, name):
     df = type_domainfield()
     r = type_range()
     dim = get_dimension()
+    if get_option("blockstructured"):
+        include_file("dune/perftool/blockstructured/blockstructuredqkfem.hh", filetag="driver")
+        return "using {} = Dune::PDELab::BlockstructuredQkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
     if isQk(element):
         include_file("dune/pdelab/finiteelementmap/qkfem.hh", filetag="driver")
         return "using {} = Dune::PDELab::QkLocalFiniteElementMap<{}, {}, {}, {}>;".format(name, gv, df, r, element.degree())
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 08e36c709c729d5563843428e9b92e5052f28556..cbcb93d3ea18cbafaca5df2bc12e168d09f6839b 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -234,7 +234,7 @@ def determine_accumulation_space(info, number):
     if get_option("blockstructured"):
         from dune.perftool.blockstructured.tools import micro_index_to_macro_index
         from dune.perftool.blockstructured.spaces import lfs_inames
-        lfsi = micro_index_to_macro_index(subel, lfs_inames(subel, ma.restriction, count=number)[0])
+        lfsi = micro_index_to_macro_index(subel, lfs_inames(subel, info.restriction, count=number)[0])
     else:
         from dune.perftool.pdelab.spaces import lfs_inames
         lfsi = Variable(lfs_iname(subel, info.restriction, count=number))
@@ -426,6 +426,9 @@ def visit_integral(integral):
     if get_option('sumfact'):
         from dune.perftool.sumfact import SumFactInterface
         interface = SumFactInterface()
+    elif get_option('blockstructured'):
+        from dune.perftool.blockstructured import BlockStructuredInterface
+        interface = BlockStructuredInterface()
     else:
         from dune.perftool.pdelab import PDELabInterface
         interface = PDELabInterface()
diff --git a/test/blockstructured/stokes/hagenpoiseuille_ref.vtu b/test/blockstructured/stokes/hagenpoiseuille_ref.vtu
deleted file mode 100644
index f6d851aca8bc95b570f13d9b26b9de691f129460..0000000000000000000000000000000000000000
--- a/test/blockstructured/stokes/hagenpoiseuille_ref.vtu
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:1d887eb3a33f5ad367149e61df7e7f82092d74b6c46b3d649c848775e2b3cbfa
-size 80161