diff --git a/python/dune/perftool/blockstructured/__init__.py b/python/dune/perftool/blockstructured/__init__.py
index 17d7025637b13eabfad5dfae37a175d1cd3299b2..9657e81c3336681285d65f0bd22c18c36155ff84 100644
--- a/python/dune/perftool/blockstructured/__init__.py
+++ b/python/dune/perftool/blockstructured/__init__.py
@@ -7,6 +7,8 @@ import dune.perftool.blockstructured.basis
 from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position
 # from dune.perftool.pdelab.geometry import to_global
 from dune.perftool.blockstructured.spaces import lfs_inames
+from dune.perftool.blockstructured.basis import (pymbolic_reference_gradient,
+                                                 pymbolic_basis)
 from dune.perftool.blockstructured.geometry import (pymbolic_jacobian_inverse_transposed,
                                                     pymbolic_jacobian_determinant,
                                                     pymbolic_facet_jacobian_determinant,
@@ -28,6 +30,16 @@ class BlockStructuredInterface(PDELabInterface):
     def lfs_inames(self, element, restriction, number=None, context=''):
         return lfs_inames(element, restriction, number, context) + sub_element_inames()
 
+    #
+    # Test and trial function related generator functions
+    #
+
+    def pymbolic_basis(self, element, restriction, number, context=''):
+        return pymbolic_basis(element, restriction, number, context)
+
+    def pymbolic_reference_gradient(self, element, restriction, number, context=''):
+        return pymbolic_reference_gradient(element, restriction, number, context)
+
     #
     # Geometry related generator functions
     #
diff --git a/python/dune/perftool/blockstructured/basis.py b/python/dune/perftool/blockstructured/basis.py
index 4751b89c033990a6d35d4891ef3b47fd3af07662..c3e42822359eb9869cddd83c47c4815f57bac3ad 100644
--- a/python/dune/perftool/blockstructured/basis.py
+++ b/python/dune/perftool/blockstructured/basis.py
@@ -18,6 +18,10 @@ from dune.perftool.pdelab.driver 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.restriction import restricted_name
+from dune.perftool.blockstructured.spaces import lfs_inames
+
+import pymbolic.primitives as prim
 
 
 # define FE basis explicitly in localoperator
@@ -61,7 +65,6 @@ def name_localbasis(leaf_element):
     return name
 
 
-@backend(interface="evaluate_basis", name="blockstructured")
 @kernel_cached
 def evaluate_basis(leaf_element, name, restriction):
     temporary_variable(name, shape=(leaf_element.degree(),), decl_method=declare_cache_temporary(leaf_element, restriction, 'Function'))
@@ -75,7 +78,16 @@ def evaluate_basis(leaf_element, name, restriction):
                 )
 
 
-@backend(interface="evaluate_grad", name="blockstructured")
+def pymbolic_basis(leaf_element, restriction, number, context=''):
+    assert leaf_element.num_sub_elements() == 0
+    name = "phi_{}".format(FEM_name_mangling(leaf_element))
+    name = restricted_name(name, restriction)
+    evaluate_basis(leaf_element, name, restriction)
+    iname = lfs_inames(leaf_element, restriction, number, context=context)[0]
+
+    return prim.Subscript(prim.Variable(name), (prim.Variable(iname), ))
+
+
 @kernel_cached
 def evaluate_reference_gradient(leaf_element, name, restriction):
     temporary_variable(name, shape=(leaf_element.degree(), 1, world_dimension()), decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian'))
@@ -87,3 +99,13 @@ def evaluate_reference_gradient(leaf_element, name, restriction):
                 assignees=frozenset({name}),
                 read_variables=frozenset({get_pymbolic_basename(qp)}),
                 )
+
+
+def pymbolic_reference_gradient(leaf_element, restriction, number, context=''):
+    assert leaf_element.num_sub_elements() == 0
+    name = "js_{}".format(FEM_name_mangling(leaf_element))
+    name = restricted_name(name, restriction)
+    evaluate_reference_gradient(leaf_element, name, restriction)
+    iname = lfs_inames(leaf_element, restriction, number, context=context)[0]
+
+    return prim.Subscript(prim.Variable(name), (prim.Variable(iname), 0))
\ No newline at end of file
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 10c944d08723ff4b56617964bfb02bb26e1a990a..7fd14d1f990f62636e5e771586c0d6302315b151 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -64,23 +64,25 @@ class PDELabInterface(object):
     # Test and trial function related generator functions
     #
 
-    def pymbolic_basis(self, element, restriction, number):
-        return pymbolic_basis(element, restriction, number)
+    # TODO: should pymbolic_basis also pass the interface?
+    def pymbolic_basis(self, element, restriction, number, context=''):
+        return pymbolic_basis(element, restriction, number, context)
 
-    def pymbolic_reference_gradient(self, element, restriction, number):
-        return pymbolic_reference_gradient(element, restriction, number)
+    # TODO: should pymbolic_reference_gradient also pass the interface?
+    def pymbolic_reference_gradient(self, element, restriction, number, context=''):
+        return pymbolic_reference_gradient(element, restriction, number, context)
 
     def pymbolic_trialfunction_gradient(self, element, restriction, tree_path):
-        return pymbolic_trialfunction_gradient(element, restriction, tree_path)
+        return pymbolic_trialfunction_gradient(self, element, restriction, tree_path)
 
     def pymbolic_apply_function_gradient(self, element, restriction, tree_path):
-        return pymbolic_apply_function_gradient(element, restriction, tree_path)
+        return pymbolic_apply_function_gradient(self, element, restriction, tree_path)
 
     def pymbolic_trialfunction(self, element, restriction, tree_path):
-        return pymbolic_trialfunction(element, restriction, tree_path)
+        return pymbolic_trialfunction(self, element, restriction, tree_path)
 
     def pymbolic_apply_function(self, element, restriction, tree_path):
-        return pymbolic_apply_function(element, restriction, tree_path)
+        return pymbolic_apply_function(self, element, restriction, tree_path)
 
     #
     # Parameter function related generator functions
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 0d6d4ea7586b8c459cc5569fd176841714031e3f..063e2645e73161cfe77715c40a6e43631cc5b04b 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, tree_path):
+def pymbolic_trialfunction_gradient(interface, element, restriction, tree_path):
     rawname = "gradu" + "_".join(str(c) for c in tree_path)
     name = restricted_name(rawname, restriction)
     container = name_coefficientcontainer(restriction)
-    evaluate_coefficient_gradient(element, name, container, restriction, tree_path)
+    evaluate_coefficient_gradient(interface, element, name, container, restriction, tree_path)
     return Variable(name)
 
 
-def pymbolic_trialfunction(element, restriction, tree_path):
+def pymbolic_trialfunction(interface, element, restriction, tree_path):
     rawname = "u" + "_".join(str(c) for c in tree_path)
     name = restricted_name(rawname, restriction)
     container = name_coefficientcontainer(restriction)
-    evaluate_coefficient(element, name, container, restriction, tree_path)
+    evaluate_coefficient(interface, element, name, container, restriction, tree_path)
     return Variable(name)
 
 
-def pymbolic_apply_function_gradient(element, restriction, tree_path):
+def pymbolic_apply_function_gradient(interface, element, restriction, tree_path):
     rawname = "gradz_func" + "_".join(str(c) for c in tree_path)
     name = restricted_name(rawname, restriction)
     container = name_applycontainer(restriction)
-    evaluate_coefficient_gradient(element, name, container, restriction, tree_path)
+    evaluate_coefficient_gradient(interface, element, name, container, restriction, tree_path)
     return Variable(name)
 
 
-def pymbolic_apply_function(element, restriction, tree_path):
+def pymbolic_apply_function(interface, element, restriction, tree_path):
     rawname = "z_func" + "_".join(str(c) for c in tree_path)
     name = restricted_name(rawname, restriction)
     container = name_applycontainer(restriction)
-    evaluate_coefficient(element, name, container, restriction, tree_path)
+    evaluate_coefficient(interface, element, name, container, restriction, tree_path)
     return Variable(name)
 
 
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index eb89b04c5c32ed0f82f5dbfd8a1305f1e43810d9..f4b92cf26830ecabbab304513000576eb99f520f 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -103,7 +103,7 @@ def pymbolic_basis(leaf_element, restriction, number, context=''):
     assert leaf_element.num_sub_elements() == 0
     name = "phi_{}".format(FEM_name_mangling(leaf_element))
     name = restricted_name(name, restriction)
-    get_backend("evaluate_basis", selector=option_switch("blockstructured"))(leaf_element, name, restriction)
+    evaluate_basis(leaf_element, name, restriction)
     iname = get_backend("lfs_inames", selector=option_switch("blockstructured"))(leaf_element, restriction, number, context=context)[0]
 
     return Subscript(Variable(name), (Variable(iname), ))
@@ -131,7 +131,7 @@ def pymbolic_reference_gradient(leaf_element, restriction, number, context=''):
     assert leaf_element.num_sub_elements() == 0
     name = "js_{}".format(FEM_name_mangling(leaf_element))
     name = restricted_name(name, restriction)
-    get_backend("evaluate_grad", selector=option_switch("blockstructured"))(leaf_element, name, restriction)
+    evaluate_reference_gradient(leaf_element, name, restriction)
     iname = get_backend("lfs_inames", selector=option_switch("blockstructured"))(leaf_element, restriction, number, context=context)[0]
 
     return Subscript(Variable(name), (Variable(iname), 0))
@@ -147,7 +147,7 @@ def shape_as_pymbolic(shape):
 
 
 @kernel_cached
-def evaluate_coefficient(element, name, container, restriction, tree_path):
+def evaluate_coefficient(interface, element, name, container, restriction, tree_path):
     from ufl.functionview import select_subelement
     sub_element = select_subelement(element, tree_path)
 
@@ -164,7 +164,7 @@ def evaluate_coefficient(element, name, container, restriction, tree_path):
 
     temporary_variable(name, shape=shape, shape_impl=shape_impl)
     lfs = name_lfs(element, restriction, tree_path)
-    basis = pymbolic_basis(leaf_element, restriction, 0, context='trial')
+    basis = interface.pymbolic_basis(leaf_element, restriction, 0, context='trial')
     index, = get_pymbolic_indices(basis)
 
     if isinstance(sub_element, (VectorElement, TensorElement)):
@@ -188,7 +188,7 @@ def evaluate_coefficient(element, name, container, restriction, tree_path):
 
 
 @kernel_cached
-def evaluate_coefficient_gradient(element, name, container, restriction, tree_path):
+def evaluate_coefficient_gradient(interface, element, name, container, restriction, tree_path):
     # First we determine the rank of the tensor we are talking about
     from ufl.functionview import select_subelement
     sub_element = select_subelement(element, tree_path)
@@ -207,7 +207,8 @@ def evaluate_coefficient_gradient(element, name, container, restriction, tree_pa
     # and proceed to call the necessary generator functions
     temporary_variable(name, shape=shape, shape_impl=shape_impl)
     lfs = name_lfs(element, restriction, tree_path)
-    basis = pymbolic_reference_gradient(leaf_element, restriction, 0, context='trialgrad')
+    from pudb import set_trace; set_trace()
+    basis = interface.pymbolic_reference_gradient(leaf_element, restriction, 0, context='trialgrad')
     index, _ = get_pymbolic_indices(basis)
 
     if isinstance(sub_element, (VectorElement, TensorElement)):