diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 172f7681771d656f37ec86934c4fdf2482babf1d..0e951fed920aa61ea5e311ebbf5c6efd32f0eaaa 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -23,6 +23,14 @@ def name_trialfunctionspace(*a):
     #TODO
     return "lfsu"
 
+def name_argumentspace(arg, *others):
+    if arg.number() == 0:
+        return name_testfunctionspace(arg, *others)
+    if arg.number() == 1:
+        return name_trialfunctionspace(arg, *others)
+    # We should never encounter an argument other than 0 or 1
+    assert False
+
 @dune_symbol
 def name_residual():
     return "r"
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 99f7ea53d31504119b1d346cd291c2a1cc8ebe26..7af8a8a53df295937dcac66e707c6296aea113cd 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -1,11 +1,24 @@
 from pytools import memoize
 from dune.perftool.options import get_option
 from dune.perftool.generation import generator_factory
+from dune.perftool.pdelab import dune_symbol
 
 # Define the generators used in-here
 operator_include = generator_factory(item_tags=("pdelab", "include", "operator"), on_store=lambda i: "#include<{}>".format(i), no_deco=True)
 base_class = generator_factory(item_tags=("pdelab", "baseclass", "operator"), counted=True, no_deco=True)
-initializer_list = generator_factory(item_tags=("pdelab", "initializer", "operator"), counted=True)
+initializer_list = generator_factory(item_tags=("pdelab", "initializer", "operator"), counted=True, no_deco=True)
+# TODO definition
+#private_member = generator_factory(item_tags=("pdelab", "member", "privatemember"))
+
+@generator_factory(item_tags=("pdelab", "constructor"), counted=True)
+def constructor_parameter(_type, name):
+    return "{} {}".format(_type, name)
+
+@dune_symbol
+def name_initree_constructor():
+    operator_include('dune/common/parametertree.hh')
+    constructor_parameter("const Dune::ParameterTree&", "iniParams")
+    return "iniParams"
 
 @memoize
 def measure_specific_details(measure):
@@ -15,7 +28,14 @@ def measure_specific_details(measure):
     if measure == "cell":
         base_class('Dune::PDELab::FullVolumePattern')
         if get_option("numerical_jacobian"):
-            base_class("Dune::PDELab::NumericalJacobianVolume<{}>".format())
+            # Add a base class
+            from dune.perftool.pdelab.driver import type_localoperator
+            loptype = type_localoperator()
+            base_class("Dune::PDELab::NumericalJacobianVolume<{}>".format(loptype))
+
+            # Add the initializer list for that base class
+            ini = name_initree_constructor()
+            initializer_list("Dune::PDELab::NumericalJacobianVolume<{}>({}.get(\"numerical_epsilon.volume\", 1e-9))".format(loptype, ini))
 
         ret["residual_signature"] = ['template<typename EG, typename LFSV0, typename X, typename LFSV1, typename R>',
                                      'void alpha_volume(const EG& eg, const LFSV0& lfsv0, const X& x, const LFSV1& lfsv1, R& r) const']
@@ -61,42 +81,33 @@ def generate_term(integrand=None, measure=None):
     from dune.perftool.generation import retrieve_cache_items
     from dune.perftool.target import DuneTarget
     domains = [i for i in retrieve_cache_items("domain")]
-    instructions = [i for i in retrieve_cache_items("cinstruction")] \
-                 + [i for i in retrieve_cache_items("exprinstruction")]
+    instructions = [i for i in retrieve_cache_items("instruction")]
     temporaries = {i.name:i for i in retrieve_cache_items("temporary")}
-    preambles = [i[1] for i in retrieve_cache_items("preamble")]
-
-    print "Printing the information that we found:\n\nDomains:"
-    for d in domains:
-        print d
-    print "\nInstructions:"
-    for i in instructions:
-        print i
-    print "\nTemporaries:"
-    for t, v in temporaries.items():
-        print "{}: {}".format(t, v)
-    print "\nPreambles:"
-    for p in preambles:
-        print p
-
-#     import loopy
-#     import numpy
-#     from loopy import make_kernel, preprocess_kernel
-#     from dune.perftool.target import DuneTarget
-#     k = make_kernel(domains, instructions,
-#                 [loopy.GlobalArg("grad_a0", numpy.float64, shape=("argi_n", "dim")),
-#                  loopy.GlobalArg("grad_c0", numpy.float64, shape=("dim",)),
-#                  loopy.ValueArg("dim", numpy.int32),
-#                  loopy.ValueArg("q_n", numpy.int32),
-#                  loopy.ValueArg("argi_n", numpy.int32)],
-#                     temporary_variables=temporaries, target=DuneTarget())
-#     k = preprocess_kernel(k)
-    # Create the kernel
-    from loopy import make_kernel, preprocess_kernel, add_dtypes
-    kernel = make_kernel(domains, instructions, temporary_variables=temporaries, preambles=preambles, target=DuneTarget())
+    preambles = [i for i in retrieve_cache_items("preamble")]
+
+# TODO it might be necessary to list the arguments and their shape manually to avoid automatic shape detection!
+    import loopy
     import numpy
-    kernel = add_dtypes(kernel, dict(grad_a0=numpy.float64, grad_c0=numpy.float64))
-    kernel = preprocess_kernel(kernel)
+    from loopy import make_kernel, preprocess_kernel
+    from dune.perftool.target import DuneTarget
+    k = make_kernel(domains, instructions,
+                [loopy.GlobalArg("grad_a0", numpy.float64, shape=("argi_n", "dim")),
+                 loopy.GlobalArg("grad_a1", numpy.float64, shape=("argj_n", "dim")),
+                 loopy.GlobalArg("grad_c0", numpy.float64, shape=("dim")),
+                 loopy.ValueArg("dim", numpy.int32),
+                 loopy.ValueArg("q_n", numpy.int32),
+                 loopy.ValueArg("argi_n", numpy.int32),
+                 loopy.ValueArg("argj_n", numpy.int32)],
+                    temporary_variables=temporaries, target=DuneTarget())
+    k = preprocess_kernel(k)
+    kernel = k
+
+#     # Create the kernel
+#     from loopy import make_kernel, preprocess_kernel, add_dtypes
+#     kernel = make_kernel(domains, instructions, temporary_variables=temporaries, preambles=preambles, target=DuneTarget())
+#     import numpy
+#     kernel = add_dtypes(kernel, dict(grad_a0=numpy.float64, grad_c0=numpy.float64))
+#     kernel = preprocess_kernel(kernel)
 
     # Return the actual code (might instead return kernels...)
     from loopy import generate_code
@@ -127,7 +138,15 @@ def generate_localoperator(ufldata, operatorfile):
     if get_option("numerical_jacobian"):
         operator_include("dune/pdelab/localoperator/defaultimp.hh")
     else:
-        pass
+        from ufl import derivative
+        from ufl.algorithms import expand_derivatives
+        jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
+
+        for integral in jacform.integrals():
+            body = generate_term(integrand=integral.integrand(), measure=integral.integral_type())
+            signature = measure_specific_details(integral.integral_type())["jacobian_signature"]
+            operator_methods.append((signature, body))
+
 
     # TODO: JacobianApply for matrix-free computations.
 
@@ -136,9 +155,6 @@ def generate_localoperator(ufldata, operatorfile):
     operator_include('dune/pdelab/localoperator/idefault.hh')
     operator_include('dune/pdelab/localoperator/flags.hh')
     operator_include('dune/pdelab/localoperator/pattern.hh')
-    operator_include('dune/common/parametertree.hh')
     operator_include('dune/geometry/quadraturerules.hh')
 
     base_class('Dune::PDELab::LocalOperatorDefaultFlags')
-
-    print operator_methods[0]
diff --git a/python/dune/perftool/transformer.py b/python/dune/perftool/transformer.py
index 14269bf0a4a913447898891ba0a19164776722fa..ab065db3b20d6681abc20e07a9f42412beda613e 100644
--- a/python/dune/perftool/transformer.py
+++ b/python/dune/perftool/transformer.py
@@ -60,17 +60,14 @@ class UFLVisitor(MultiFunction):
         # Visit the expression. The return value will be a pymbolic expression.
         loopyexpr = self.call(o)
 
-        # TODO no jacobian support yet!
-        assert len(self.arguments) == 1
-
         # The expression is completely processes, an accumulation instruction can be issued!
         # Have some data structures for the analysis of the given arguments.
         accumargs = []
 
         # visit the given arguments and extract necessary information
-        for arg, index, grad, restriction in self.arguments:
-            from dune.perftool.pdelab.argument import name_testfunctionspace
-            accumargs.append(name_testfunctionspace(arg, index, restriction))
+        for arg, index, grad, restriction in sorted(self.arguments, key=lambda x: x[0].number()):
+            from dune.perftool.pdelab.argument import name_argumentspace
+            accumargs.append(name_argumentspace(arg, index, restriction))
             accumargs.append(argument_iname(arg))
 
         # Explicitly state the dependency on the quadrature iname
@@ -110,7 +107,6 @@ class UFLVisitor(MultiFunction):
         self.inames.append(index)
 
         return Subscript(Variable(name), (Variable(index),))
-#         return Subscript(Variable(name), Variable(index))
 
     def coefficient(self, o):
         assert len(o.operands()) == 0
@@ -149,7 +145,6 @@ class UFLVisitor(MultiFunction):
             return ret
         else:
             from dune.perftool.pdelab import name_index
-#            return Subscript(self.call(o.operands()[0]), Variable(name_index(o.operands()[1])))
             ret = self.call(o.operands()[0])
             index = Variable(name_index(o.operands()[1]))
             if isinstance(ret, Subscript):