diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 3a22422262ca80f232c5d45ca749d9a541b0b052..0cb1ef2f826969fc375c2611463e4f23d5bc3fa0 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -12,7 +12,7 @@ from dune.perftool.pdelab.argument import (pymbolic_apply_function,
                                            )
 from dune.perftool.pdelab.basis import (pymbolic_basis,
                                         pymbolic_reference_gradient,
-                                        pymbolic_evaluate_gridfunction,
+                                        pymbolic_gridfunction,
                                         )
 from dune.perftool.pdelab.geometry import (component_iname,
                                            pymbolic_cell_volume,
@@ -99,8 +99,8 @@ class PDELabInterface(object):
     def pymbolic_apply_function(self, element, restriction, index):
         return pymbolic_apply_function(self.visitor, element, restriction, index)
 
-    def pymbolic_evaluate_gridfunction(self, coeff, restriction, grad):
-        return pymbolic_evaluate_gridfunction(self.visitor, coeff, restriction, grad)
+    def pymbolic_gridfunction(self, coeff, restriction, grad):
+        return pymbolic_gridfunction(coeff, restriction, grad)
 
     #
     # Tensor expression related generator functions
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index b2593874aff3294224863a232f6ccb2728f6da89..b786c3be775e158362237e57574dc8d94a97a63e 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -231,20 +231,53 @@ def evaluate_coefficient_gradient(visitor, element, name, container, restriction
                 )
 
 
+def local_gridfunction(coeff, restriction, diffOrder):
+    from dune.perftool.pdelab.localoperator import name_gridfunction_member
+    gf = name_gridfunction_member(coeff, diffOrder)
+    name = "{}_local".format(gf)
+    define_local_gridfunction(name, gf)
+    bind_gridfunction_to_element(name, restriction)
+    return name
+
+
+@preamble
+def define_local_gridfunction(name, gf):
+    return "auto {} = localFunction({});".format(name, gf)
+
+
 @preamble
 def bind_gridfunction_to_element(gf, restriction):
     element = name_cell(restriction)
-    return "localFunction({}).bind({});".format(gf, element)
+    return "{}.bind({});".format(gf, element)
+
+
+def declare_grid_function_range(gridfunction):
+    def _decl(name, *args):
+        return "typename decltype({})::Range {};".format(gridfunction, name)
 
+    return _decl
 
-def pymbolic_evaluate_gridfunction(visitor, coeff, restriction, grad):
+@kernel_cached
+def pymbolic_evaluate_gridfunction(name, coeff, restriction, grad):
     diffOrder = 1 if grad else 0
 
-    from dune.perftool.pdelab.localoperator import name_gridfunction_member
-    gridfunction = name_gridfunction_member(coeff, diffOrder)
+    gridfunction = local_gridfunction(coeff, restriction, diffOrder)
+
+    temporary_variable(name,
+                       shape=(1,) + (world_dimension(),) * diffOrder,
+                       decl_method=declare_grid_function_range(gridfunction),
+                       managed=False,
+                       )
+
+    quadpos = get_backend(interface="qp_in_cell")(restriction)
+    instruction(code="{} = {}({});".format(name, gridfunction, quadpos),
+                assignees=frozenset({name}),
+                within_inames=frozenset(get_backend(interface="quad_inames")()),
+                within_inames_is_final=True,
+                )
 
-    bind_gridfunction_to_element(gridfunction, restriction)
 
-    # Foobar!
-    visitor.indices = None
-    return 1
+def pymbolic_gridfunction(coeff, restriction, grad):
+    name = "coeff{}{}".format(coeff.count(), "_grad" if grad else "")
+    pymbolic_evaluate_gridfunction(name, coeff, restriction, grad)
+    return Subscript(Variable(name), (0,))
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index dffc7c88be19609985234ecddb4cdf1675ff60fa..11c302d9744fdc70eb30395f5dcd0697558b6cc2 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -179,7 +179,7 @@ def define_gridfunction_member(name, coeff, diffOrder):
         newtype = "Dune::PDELab::DiscreteGridViewFunction<typename {0}::GridFunctionSpace, typename {0}::Vector, {1}>".format(_type, diffOrder)
         params = ["{}.gridFunctionSpace()".format(param), "{}.dofs()".format(param)]
         initializer_list(name, params, classtag="operator")
-        return "{} {};".format(_type, name)
+        return "{} {};".format(newtype, name)
     else:
         initializer_list(name, [param], classtag="operator")
         return "const {}& {};".format(_type, name)
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index 153f9713784b3d03d40053d493e34f50d7a645ff..687afd5c5d95538089d030021aade9c9b5f707af 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -164,7 +164,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
             # and exports it through a getter method 'getTime'
             return prim.Call(prim.Variable("getTime"), ())
         else:
-            return self.interface.pymbolic_evaluate_gridfunction(o, restriction, self.reference_grad)
+            return self.interface.pymbolic_gridfunction(o, restriction, self.reference_grad)
 
     #
     # Handlers for all indexing related stuff
diff --git a/test/coeffeval/CMakeLists.txt b/test/coeffeval/CMakeLists.txt
index 5150d01b200fa50bf32b63b7e66e7d24c0ab58ac..fd621382ce4a73468bdfd2548bd8809b4017c8c6 100644
--- a/test/coeffeval/CMakeLists.txt
+++ b/test/coeffeval/CMakeLists.txt
@@ -8,3 +8,5 @@ add_generated_executable(UFLFILE poisson.ufl
 dune_add_test(TARGET coeffeval_poisson
               CMD_ARGS coeffeval_poisson.ini
               )
+
+dune_symlink_to_source_files(FILES coeffeval_poisson.ini)