diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 05008cba2c986a74d761df93dfb9924078be1fbe..bb9efa516b971c03007322216aa59e5cff34523b 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -43,7 +43,6 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
         # Reset some state variables that are reinitialized for each accumulation term
         self.argshape = 0
         self.transpose_necessary = False
-        self.redinames = ()
         self.inames = []
         self.substitution_rules = []
 
@@ -282,7 +281,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
         else:
             return Subscript(aggr, use_indices)
 
-    def index_sum(self, o):
+    def index_sum(self, o, additional_inames=()):
         # There is three scenarios here:
         # * This is a normal IndexSum:
         #   We should collect the reduction inames, go into recursion and issue
@@ -294,7 +293,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
 
         # Get the iname for the reduction index
         ind = o.ufl_operands[1][0]
-        self.redinames = self.redinames + (ind,)
+        redinames = additional_inames + (ind,)
         shape = o.ufl_operands[0].ufl_index_dimensions[0]
         from dune.perftool.pdelab import name_index
         domain(name_index(ind), shape)
@@ -302,31 +301,28 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
         # If the left operand is an index sum to, we do it in one reduction
         from ufl.classes import IndexSum
         if isinstance(o.ufl_operands[0], IndexSum):
-            return self.call(o.ufl_operands[0])
+            return self.index_sum(o.ufl_operands[0], additional_inames=redinames)
         else:
             from loopy import Reduction
 
             # Recurse to get the summation expression
             term = self.call(o.ufl_operands[0])
-            self.redinames = tuple(i for i in self.redinames if i not in self.dimension_indices)
-            if len(self.redinames) > 0:
-                ret = Reduction("sum", tuple(name_index(ind) for ind in self.redinames), term)
-                name = get_temporary_name()
+            redinames = tuple(i for i in redinames if i not in self.dimension_indices)
+            if len(redinames) > 0:
+                ret = Reduction("sum", tuple(name_index(ind) for ind in redinames), term)
+#                 name = get_temporary_name()
                 # Generate a substitution rule for this one.
-                from loopy import SubstitutionRule
-                self.substitution_rules.append(SubstitutionRule(name,
-                                                                (),
-                                                                ret
-                                                                )
-                                               )
-
-                ret = Variable(name)
+#                 from loopy import SubstitutionRule
+#                 self.substitution_rules.append(SubstitutionRule(name,
+#                                                                 (),
+#                                                                 ret
+#                                                                 )
+#                                                )
+#
+#                 ret = Variable(name)
             else:
                 ret = term
 
-            # Reset the reduction inames for future indexsums
-            self.redinames = ()
-
             return ret
 
     def _index_or_fixed_index(self, index):
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 611f996d7b21ec8821af5a8ed9aa2b03edcd1fcf..564c2b96c25a249cff525a43bcd70d637337d013 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -25,15 +25,6 @@ from pymbolic.primitives import Call, Subscript, Variable
 import loopy
 
 
-def name_testfunction_gradient(element, restriction):
-    return name_basis_gradient(element, restriction)
-
-
-def name_testfunction(element, restriction):
-    it = get_global_context_value("integral_type")
-    if it == 'exterior_facet':
-        restriction = Restriction.NEGATIVE
-    return name_basis(element, restriction)
 
 
 def name_trialfunction_gradient(element, restriction, component):
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index 912de6ab6165b39810a77e59f9e39a822295fef8..416e1cbb29b000a3c1815445faee7b9e854bec49 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -288,55 +288,6 @@ def name_reference_gradient(leaf_element, restriction):
     return name
 
 
-@cached
-def evaluate_basis_gradient(leaf_element, name, restriction):
-    lfs = name_leaf_lfs(leaf_element, restriction)
-    temporary_variable(name, shape=(name_lfs_bound(lfs), name_dimension()), shape_impl=('vec', 'fv'))
-    jac = name_jacobian_inverse_transposed(restriction)
-    index = lfs_iname(leaf_element, restriction, context='transformgrads')
-    reference_gradients = name_reference_gradient(leaf_element, restriction)
-
-    # Evaluate product of jacobianInvereseTransposed and gradient on reference element with loopy
-    dimindex_0 = dimension_iname(context='gradient_transformation', count=0)
-    from dune.perftool.options import get_option
-    if get_option('diagonal_transformation_matrix'):
-        # dimindex_1 = dimension_iname(context='gradient_transformation', count=1)
-        reduction_expr = Product((
-            Subscript(Variable(jac), (Variable(dimindex_0), Variable(dimindex_0))),
-            Subscript(Variable(reference_gradients), (Variable(index), 0, Variable(dimindex_0)))))
-
-        assignee = Subscript(Variable(name), (Variable(index), Variable(dimindex_0)))
-
-        # Create exrpession instruction
-        forced_iname_deps = frozenset({quadrature_iname()}).union(frozenset({index})).union(frozenset({dimindex_0}))
-        instruction(
-            expression=reduction_expr,
-            assignee=assignee,
-            forced_iname_deps=forced_iname_deps,
-            forced_iname_deps_is_final=True)
-    else:
-        dimindex_1 = dimension_iname(context='gradient_transformation', count=1)
-        reduction_expr = Product((
-            Subscript(Variable(jac), (Variable(dimindex_0), Variable(dimindex_1))),
-            Subscript(Variable(reference_gradients), (Variable(index), 0, Variable(dimindex_1)))))
-
-        assignee = Subscript(Variable(name), (Variable(index), Variable(dimindex_0)))
-
-        # Create exrpession instruction
-        forced_iname_deps = frozenset({quadrature_iname()}).union(frozenset({index})).union(frozenset({dimindex_0}))
-        instruction(
-            expression=Reduction("sum", dimindex_1, reduction_expr, allow_simultaneous=True),
-            assignee=assignee,
-            forced_iname_deps=forced_iname_deps,
-            forced_iname_deps_is_final=True)
-
-
-def name_basis_gradient(leaf_element, restriction):
-    assert leaf_element.num_sub_elements() == 0
-    name = "gradphi_{}".format(FEM_name_mangling(leaf_element))
-    name = restricted_name(name, restriction)
-    evaluate_basis_gradient(leaf_element, name, restriction)
-    return name
 
 
 def shape_as_pymbolic(shape):