diff --git a/python/dune/codegen/blockstructured/basis.py b/python/dune/codegen/blockstructured/basis.py
index 5af8bb48f2bb625a0b336c5cae6d42863f0063c5..9a0c67afff4fdac59eb86a0b6e4f332c6a9aa79f 100644
--- a/python/dune/codegen/blockstructured/basis.py
+++ b/python/dune/codegen/blockstructured/basis.py
@@ -7,7 +7,7 @@ from dune.codegen.generation import (basis_mixin,
                                      globalarg,
                                      class_member,
                                      initializer_list,
-                                     include_file,)
+                                     include_file, preamble)
 from dune.codegen.tools import get_pymbolic_basename, get_pymbolic_indices
 from dune.codegen.loopy.target import type_floatingpoint
 from dune.codegen.pdelab.basis import (GenericBasisMixin,
@@ -38,11 +38,22 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
         assert not isinstance(element, MixedElement)
         name = "phi_{}".format(FEM_name_mangling(element))
         name = restricted_name(name, restriction)
+        self.init_basis_cache(element)
         self.evaluate_basis(element, name, restriction)
         inames = self.lfs_inames(element, restriction, number, context=context)
 
         return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, element.degree() + 1), 0))
 
+    @preamble(kernel='operator')
+    def init_basis_cache(self, element):
+        cache = name_localbasis_cache(element)
+        localbasis = name_localbasis(element)
+        qp_name = get_pymbolic_basename(self.quadrature_position_in_micro())
+        return ["for (int i=0; i < {}.size(); ++i)".format(qp_name),
+                "{",
+                "  {}.evaluateFunction({}[i], {});".format(cache, qp_name, localbasis),
+                "}"]
+
     @kernel_cached
     def evaluate_basis(self, element, name, restriction):
         temporary_variable(name, shape=((element.degree() + 1)**world_dimension(), 1),
@@ -60,11 +71,22 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
         assert not isinstance(element, MixedElement)
         name = "js_{}".format(FEM_name_mangling(element))
         name = restricted_name(name, restriction)
+        self.init_gradient_cache(element)
         self.evaluate_reference_gradient(element, name, restriction)
         inames = self.lfs_inames(element, restriction, number, context=context)
 
         return prim.Subscript(prim.Variable(name), (tensor_index_to_sequential_index(inames, element.degree() + 1), 0))
 
+    @preamble(kernel='operator')
+    def init_gradient_cache(self, element):
+        cache = name_localbasis_cache(element)
+        localbasis = name_localbasis(element)
+        qp_name = get_pymbolic_basename(self.quadrature_position_in_micro())
+        return ["for (int i=0; i < {}.size(); ++i)".format(qp_name),
+                "{",
+                "  {}.evaluateJacobian({}[i], {});".format(cache, qp_name, localbasis),
+                "}"]
+
     @kernel_cached
     def evaluate_reference_gradient(self, element, name, restriction):
         temporary_variable(name, shape=((element.degree() + 1)**world_dimension(), 1, world_dimension()),