diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index e1ae4c43b0ab23bcc9a8448bad3b1b1f962f3aeb..571bdef3e3df1a5d10d6d2b60a54d0cab8dcc5ea 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -757,9 +757,11 @@ def local_operator_default_settings(operator, form):
 def generate_residual_kernels(form):
     logger = logging.getLogger(__name__)
     with global_context(form_type='residual'):
+        operator_kernels = {}
+
         # Generate the necessary residual methods
         for measure in set(i.integral_type() for i in form.integrals()):
-            logger.info("generate_localoperator_kernels: measure {}".format(measure))
+            logger.info("generate_residual_kernels: measure {}".format(measure))
             with global_context(integral_type=measure):
                 enum_pattern()
                 pattern_baseclass()
@@ -810,7 +812,9 @@ def generate_residual_kernels(form):
                                              classtag="operator",
                                              )
 
-                return {(measure, 'residual'): kernel}
+            operator_kernels[(measure, 'residual')] =  kernel
+
+        return operator_kernels
 
 
 def generate_jacobian_kernels(form, original_form):
@@ -826,7 +830,7 @@ def generate_jacobian_kernels(form, original_form):
     with global_context(form_type="jacobian"):
         if get_form_option("generate_jacobians"):
             for measure in set(i.integral_type() for i in jacform.integrals()):
-                logger.info("generate_localoperator_kernels: measure {}".format(measure))
+                logger.info("generate_jacobian_kernels: measure {}".format(measure))
                 with global_context(integral_type=measure):
                     from dune.perftool.pdelab.signatures import assembler_routine_name
                     with global_context(kernel=assembler_routine_name()):
@@ -874,6 +878,10 @@ def generate_jacobian_kernels(form, original_form):
     return operator_kernels
 
 
+def generate_control_kernels(forms):
+    pass
+
+
 def generate_localoperator_kernels(operator):
     logger = logging.getLogger(__name__)
 
@@ -890,8 +898,7 @@ def generate_localoperator_kernels(operator):
         # objective function w.r.t the ansatz function to get the
         # final residual formulation of the adjoint.
         #
-        # TODO: This might only be true for linear problems. Adjust
-        # documentation as knowledge improves ;)
+        # Might not be true in all cases but works for the simple ones.
         assert get_form_option("objective_function") is not None
         assert get_form_option("control") is False
 
@@ -937,12 +944,18 @@ def generate_localoperator_kernels(operator):
             return multi_index
 
         forms = []
+        element = original_form.coefficients()[0].ufl_element()
+        coeff = Coefficient(element, count=3)
         for control in controls:
             shape = control.ufl_shape
             flat_length = np.prod(shape)
             for i in range(flat_length):
                 c = control[_unravel(i, shape)]
                 control_form = diff(original_form, control)
+                control_form = action(control_form, coeff)
+                objective = data.object_by_name[get_form_option("objective_function")]
+                objective_gradient = diff(objective, control)
+                control_form = control_form + objective_gradient
                 forms.append(preprocess_form(control_form).preprocessed_form)
 
         # Used to create local operator default settings