diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index e11e659fc0c8a4f2e3045cb13485c492c626bb2b..26de3b67de909ee8156d33aefafcaff3727c718d 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -1783,7 +1783,7 @@ def generate_driver(formdatas, data):
         # In case of operator conunting we only assemble the matrix and evaluate the residual
         #assemble_matrix_timer()
         evaluate_residual_timer()
-#        apply_jacobian_timer()
+        apply_jacobian_timer()
     elif is_stationary():
         # We could also use solve if we are not interested in visualization
         vtkoutput()
diff --git a/python/dune/perftool/sumfact/__init__.py b/python/dune/perftool/sumfact/__init__.py
index 1900e7b51987870a2a5baeadac09f4f01dfd2a22..60b120ce95e9650b9a358632f6d3bad7249a4d25 100644
--- a/python/dune/perftool/sumfact/__init__.py
+++ b/python/dune/perftool/sumfact/__init__.py
@@ -1,3 +1,6 @@
+from dune.perftool.pdelab.argument import (name_applycontainer,
+                                           name_coefficientcontainer,
+                                           )
 from dune.perftool.sumfact.quadrature import (quadrature_inames,
                                               quadrature_weight,
                                               pymbolic_quadrature_position_in_cell,
@@ -6,8 +9,8 @@ from dune.perftool.sumfact.quadrature import (quadrature_inames,
 from dune.perftool.sumfact.basis import (lfs_inames,
                                          pymbolic_basis,
                                          pymbolic_reference_gradient,
-                                         pymbolic_trialfunction,
-                                         pymbolic_trialfunction_gradient,
+                                         pymbolic_coefficient,
+                                         pymbolic_coefficient_gradient,
                                          )
 import dune.perftool.sumfact.switch
 
@@ -25,12 +28,22 @@ class SumFactInterface(PDELabInterface):
         return pymbolic_reference_gradient(element, restriction, number)
 
     def pymbolic_trialfunction_gradient(self, element, restriction, component, visitor=None):
-        ret, indices = pymbolic_trialfunction_gradient(element, restriction, component, visitor)
+        ret, indices = pymbolic_coefficient_gradient(element, restriction, component, name_coefficientcontainer, visitor)
         visitor.indices = indices
         return ret
 
     def pymbolic_trialfunction(self, element, restriction, component, visitor=None):
-        ret, indices = pymbolic_trialfunction(element, restriction, component, visitor)
+        ret, indices = pymbolic_coefficient(element, restriction, component, name_coefficientcontainer, visitor)
+        visitor.indices = indices
+        return ret
+
+    def pymbolic_apply_function_gradient(self, element, restriction, component, visitor=None):
+        ret, indices = pymbolic_coefficient_gradient(element, restriction, component, name_applycontainer, visitor)
+        visitor.indices = indices
+        return ret
+
+    def pymbolic_apply_function(self, element, restriction, component, visitor=None):
+        ret, indices = pymbolic_coefficient(element, restriction, component, name_applycontainer, visitor)
         visitor.indices = indices
         return ret
 
diff --git a/python/dune/perftool/sumfact/basis.py b/python/dune/perftool/sumfact/basis.py
index 3c04cd09b7b4e97005ef90f0eab8e3a640522054..e9b5a6ceb39b06b6757b96c76d7b1d4e7ec49ff8 100644
--- a/python/dune/perftool/sumfact/basis.py
+++ b/python/dune/perftool/sumfact/basis.py
@@ -53,7 +53,7 @@ def name_sumfact_base_buffer():
 
 
 @kernel_cached
-def pymbolic_trialfunction_gradient(element, restriction, component, visitor):
+def pymbolic_coefficient_gradient(element, restriction, component, coeff_func, visitor):
     rawname = "gradu" + "_".join(str(c) for c in component)
     name = restricted_name(rawname, restriction)
 
@@ -92,11 +92,11 @@ def pymbolic_trialfunction_gradient(element, restriction, component, visitor):
 
         if get_option('fastdg'):
             # Name of direct input, shape and globalarg is set in sum_factorization_kernel
-            direct_input = name_coefficientcontainer(restriction)
+            direct_input = coeff_func(restriction)
         else:
             direct_input = None
             # Setup the input!
-            setup_theta(inp, element, restriction, component, index)
+            setup_theta(inp, element, restriction, component, index, coeff_func)
 
         # Add a sum factorization kernel that implements the
         # evaluation of the gradients of basis functions at quadrature
@@ -136,7 +136,7 @@ def pymbolic_trialfunction_gradient(element, restriction, component, visitor):
 
 
 @kernel_cached
-def pymbolic_trialfunction(element, restriction, component, visitor):
+def pymbolic_coefficient(element, restriction, component, coeff_func, visitor):
     # Get geometric dimension
     dim = world_dimension()
 
@@ -161,11 +161,11 @@ def pymbolic_trialfunction(element, restriction, component, visitor):
 
     if get_option('fastdg'):
         # Name of direct input, shape and globalarg is set in sum_factorization_kernel
-        direct_input = name_coefficientcontainer(restriction)
+        direct_input = coeff_func(restriction)
     else:
         direct_input = None
         # Setup the input!
-        setup_theta(inp, element, restriction, component, index)
+        setup_theta(inp, element, restriction, component, index, coeff_func)
 
     # Add a sum factorization kernel that implements the evaluation of
     # the basis functions at quadrature points (stage 1)
diff --git a/python/dune/perftool/sumfact/sumfact.py b/python/dune/perftool/sumfact/sumfact.py
index ccb4b38ea61ecfc0c912bb06cebe4295dadbdd2b..ab5ce9e1500baf2070cc0b8ad69d30ede8ec713e 100644
--- a/python/dune/perftool/sumfact/sumfact.py
+++ b/python/dune/perftool/sumfact/sumfact.py
@@ -95,7 +95,7 @@ def accum_iname(restriction, bound, i):
     return sumfact_iname(bound, "accum")
 
 
-def setup_theta(inp, element, restriction, component, index):
+def setup_theta(inp, element, restriction, component, index, coeff_func):
     if index is None:
         index = ()
     else:
@@ -103,7 +103,7 @@ def setup_theta(inp, element, restriction, component, index):
     # Write initial coefficients into buffer
     lfs = name_lfs(element, restriction, component)
     basisiname = sumfact_iname(name_lfs_bound(lfs), "basis")
-    container = name_coefficientcontainer(restriction)
+    container = coeff_func(restriction)
     coeff = pymbolic_coefficient(container, lfs, basisiname)
     assignee = Subscript(Variable(inp), (Variable(basisiname),) + index)
     return instruction(assignee=assignee,