diff --git a/python/dune/perftool/sumfact/amatrix.py b/python/dune/perftool/sumfact/amatrix.py
index 1b56523a2ce8f950e4ef5c6bcd821eb2971855db..f01d3bf36f2533e831eedea60d223421c2a5a80e 100644
--- a/python/dune/perftool/sumfact/amatrix.py
+++ b/python/dune/perftool/sumfact/amatrix.py
@@ -193,7 +193,7 @@ def sort_quadrature_points_weights():
 
 
 @constructor_block("operator")
-def construct_theta(name, transpose):
+def construct_theta(name, transpose, derivative):
     # Make sure that the quadrature points are sorted
     sort_quadrature_points_weights()
 
@@ -205,10 +205,11 @@ def construct_theta(name, transpose):
     qp = name_oned_quadrature_points()
 
     access = "j,i" if transpose else "i,j"
+    basispol = "dp" if derivative else "p"
 
     return ["for (std::size_t i=0; i<{}; i++){{".format(shape[0]),
             "  for (std::size_t j=0; j<{}; j++){{".format(shape[1]),
-            "    {}.colmajoraccess({}) = {}.p(j,{}[i]);".format(name, access, polynomials, qp),
+            "    {}.colmajoraccess({}) = {}.{}(j,{}[i]);".format(name, access, polynomials, basispol, qp),
             "  }",
             "}"]
 
@@ -228,10 +229,10 @@ def type_theta():
 
 
 @class_member("operator")
-def define_theta(name, shape, transpose):
+def define_theta(name, shape, transpose, derivative):
     theta_type = type_theta()
     initializer_list(name, [str(axis) for axis in shape], classtag="operator")
-    construct_theta(name, transpose)
+    construct_theta(name, transpose, derivative)
     return "{} {};".format(theta_type, name)
 
 
@@ -239,7 +240,7 @@ def name_theta():
     name = "Theta"
     shape = (quadrature_points_per_direction(), basis_functions_per_direction())
     globalarg(name, shape=shape, dtype=numpy.float64, dim_tags="f,f")
-    define_theta(name, shape, False)
+    define_theta(name, shape, False, False)
     return name
 
 
@@ -247,5 +248,13 @@ def name_theta_transposed():
     name = "ThetaT"
     shape = (basis_functions_per_direction(), quadrature_points_per_direction())
     globalarg(name, shape=shape, dtype=numpy.float64, dim_tags="f,f")
-    define_theta(name, shape, True)
+    define_theta(name, shape, True, False)
+    return name
+
+
+def name_dtheta():
+    name = "dTheta"
+    shape = (quadrature_points_per_direction(), basis_functions_per_direction())
+    globalarg(name, shape=shape, dtype=numpy.float64, dim_tags="f,f")
+    define_theta(name, shape, False, True)
     return name
diff --git a/python/dune/perftool/sumfact/basis.py b/python/dune/perftool/sumfact/basis.py
index fb4b411c9d0b097a334710cc427514850dad4094..455d806a251ebb2a48cf4c3630bc87b19b0173a1 100644
--- a/python/dune/perftool/sumfact/basis.py
+++ b/python/dune/perftool/sumfact/basis.py
@@ -14,6 +14,7 @@ from dune.perftool.generation import (backend,
 from dune.perftool.sumfact.amatrix import (AMatrix,
                                            basis_functions_per_direction,
                                            ColMajorAccess,
+                                           name_dtheta,
                                            name_theta,
                                            quadrature_points_per_direction,
                                            )
@@ -33,6 +34,11 @@ import pymbolic.primitives as prim
 
 @cached
 def pymbolic_trialfunction_gradient(element, restriction, component):
+
+    theta = name_theta()
+    dtheta = name_dtheta()
+
+
     # palpo TODO -> copied from argument.py
     rawname = "gradu" + "_".join(str(c) for c in component)
     name = restricted_name(rawname, restriction)