diff --git a/python/dune/perftool/blockstructured/basis.py b/python/dune/perftool/blockstructured/basis.py
index ced039c976d3dfa440d0909c9fd366ee75892fb6..5fb74b6b46a531ec64fbd1ad5754a807cecabb96 100644
--- a/python/dune/perftool/blockstructured/basis.py
+++ b/python/dune/perftool/blockstructured/basis.py
@@ -8,16 +8,21 @@ from dune.perftool.generation import (backend,
                                       initializer_list,
                                       include_file,)
 from dune.perftool.tools import get_pymbolic_basename
+# TODO clean up some imports
 from dune.perftool.pdelab.basis import (declare_cache_temporary,
                                         name_localbasis_cache,
-                                        type_localbasis
-                                        )
+                                        type_localbasis,
+                                        restricted_name,
+                                        FEM_name_mangling)
 from dune.perftool.pdelab.driver import (basetype_range,
                                          isPk,
                                          isQk)
 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.blockstructured.spaces import lfs_inames
+import pymbolic.primitives as prim
+
 
 
 @backend(interface="typedef_localbasis", name="blockstructured")
@@ -46,6 +51,19 @@ def typedef_localbasis(element, name):
     return "using {} = Dune::{};".format(name, basis_type)
 
 
+@class_member(classtag="operator")
+def define_localbasis(leaf_element, name):
+    localBasis_type = type_localbasis(leaf_element)
+    initializer_list(name, (), classtag="operator")
+    return "const {} {};".format(localBasis_type, name)
+
+
+def name_localbasis(leaf_element):
+    name = "microElementBasis"
+    globalarg(name)
+    define_localbasis(leaf_element, name)
+    return name
+
 
 @backend(interface="evaluate_basis", name="blockstructured")
 @kernel_cached
@@ -75,15 +93,19 @@ def evaluate_reference_gradient(leaf_element, name, restriction):
                 )
 
 
-@class_member(classtag="operator")
-def define_localbasis(leaf_element, name):
-    localBasis_type = type_localbasis(leaf_element)
-    initializer_list(name, (), classtag="operator")
-    return "const {} {};".format(localBasis_type, name)
+def tensor_index_to_sequential_index(inames):
+    index = prim.Sum(tuple(prim.Variable(name)*2**i for i,name in enumerate(inames)))
+    return index
 
 
-def name_localbasis(leaf_element):
-    name = "microElementBasis"
-    globalarg(name)
-    define_localbasis(leaf_element, name)
-    return name
+@backend(interface="pymbolic_reference_gradient", name="blockstructured")
+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)
+    inames = lfs_inames(leaf_element, restriction, number, context=context)
+
+    index = tensor_index_to_sequential_index(inames)
+
+    return prim.Subscript(prim.Variable(name), (index, 0))
\ No newline at end of file
diff --git a/python/dune/perftool/blockstructured/spaces.py b/python/dune/perftool/blockstructured/spaces.py
index 42cd9174a4a167c7a002402600d8fd489089bc42..24f76456f8ee27fa4ef3c587248b4f3e855b25e9 100644
--- a/python/dune/perftool/blockstructured/spaces.py
+++ b/python/dune/perftool/blockstructured/spaces.py
@@ -1,8 +1,6 @@
 from dune.perftool.generation import (backend,
                                       domain)
 from dune.perftool.pdelab.geometry import world_dimension
-from dune.perftool.blockstructured.tools import sub_element_inames
-
 
 #TODO subelem_inames hinzufugen?
 @backend(interface="lfs_inames", name="blockstructured")
@@ -14,7 +12,11 @@ def lfs_inames(element, restriction, count=None, context=''):
         else:
             context = str(count)
 
-    name = "micro_{}_index".format(context)
-    domain(name, pow(2,world_dimension()))
+    inames = tuple()
+    dim = world_dimension()
+    for i in range(dim):
+        name = "micro_{}_index_{}".format(context, i)
+        inames = inames + (name,)
+        domain(name, 2)
 
-    return (name, ) + sub_element_inames()
+    return inames
diff --git a/python/dune/perftool/blockstructured/tools.py b/python/dune/perftool/blockstructured/tools.py
index 38b0d144c38481c4ba8202913ca951b5500b1111..5c5c272d6d538262884b3453d57aac701eb6d9a9 100644
--- a/python/dune/perftool/blockstructured/tools.py
+++ b/python/dune/perftool/blockstructured/tools.py
@@ -73,13 +73,21 @@ def micro_index_to_macro_index(index):
                     within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
                     predicates=predicate("!=")
                     )
+    #
+    # modified_index = prim.Sum((prim.Variable(x), index))
+    # if local_dimension() > 1:
+    #     modified_index = prim.Sum((modified_index, prim.Product(((k+1), prim.Variable(y)))))
+    #     index_div_2 = prim.FloorDiv(index, 2)
+    #     index_div_2 = prim.Product((index_div_2, k-1))
+    #     modified_index = prim.Sum((modified_index, index_div_2))
+    #     if local_dimension() > 2:
+    #         modified_index = prim.Sum((modified_index, prim.Variable(z)))
 
     modified_index = prim.Sum((prim.Variable(x), index))
-    if local_dimension()>1:
-        modified_index = prim.Sum((modified_index, prim.Product(((k+1), prim.Variable(y)))))
-    index_div_2 = prim.FloorDiv(index, 2)
-    index_div_2 = prim.Product((index_div_2, k-1))
-    modified_index = prim.Sum((modified_index, index_div_2))
+    if local_dimension() > 1:
+        modified_index = prim.Sum((modified_index, prim.Product((prim.Variable(y),k+1))))
+        if local_dimension() > 2:
+            modified_index = prim.Sum((modified_index, prim.Product((prim.Variable(z), (k+1)**2))))
 
     return modified_index
 
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index 756a8d711017079a1a4266757c237d928f6e8046..7a107c80de03ad87138875de06ea68525154f40a 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -199,8 +199,9 @@ 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')
-    index, _ = get_pymbolic_indices(basis)
+    basis = get_backend("pymbolic_reference_gradient", selector=option_switch("blockstructured"))(leaf_element, restriction, 0, context='trialgrad')
+    # basis = pymbolic_reference_gradient(leaf_element, restriction, 0, context='trialgrad')
+    index,_ = get_pymbolic_indices(basis)
 
     if isinstance(sub_element, (VectorElement, TensorElement)):
         lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]), symmetry=element.symmetry())