diff --git a/python/dune/perftool/ufl/extract_accumulation_terms.py b/python/dune/perftool/ufl/extract_accumulation_terms.py
index 5a8a4d7a1405d1b5298026c9edc845a85110bef9..cefc82751386f521d64ee66b4a3a9ae59353aff9 100644
--- a/python/dune/perftool/ufl/extract_accumulation_terms.py
+++ b/python/dune/perftool/ufl/extract_accumulation_terms.py
@@ -59,7 +59,7 @@ def split_into_accumulation_terms(expr):
 
     # Extract a list of modified terminals for the ansatz function
     # in jacobian forms.
-    all_jacobian_args = extract_modified_arguments(expr, argnumber=1)
+    all_jacobian_args = extract_modified_arguments(expr, argnumber=1, do_index=False)
 
     for test_arg in test_args:
         # Do this as a multi step replacement procedure to avoid UFL nagging
@@ -110,7 +110,6 @@ def split_into_accumulation_terms(expr):
         # Get all appearances of test functions with their indices
         indexed_test_args = extract_modified_arguments(replace_expr, argnumber=0, do_index=True)
         for indexed_test_arg in indexed_test_args:
-            # from pudb import set_trace; set_trace()
             if indexed_test_arg.index:
                 # If the test function is indexed, create a new multiindex of this shape
                 # -> (m,n) in the example above
@@ -133,20 +132,31 @@ def split_into_accumulation_terms(expr):
 
         # 5) Further split according to trial function in jacobian terms
         if all_jacobian_args:
-            # TODO -> Jacobians not yet implemented!
-            assert(False)
-            jac_args = extract_modified_arguments(replace_expr, argnumber=1)
-
-            for restriction in (Restriction.NONE, Restriction.POSITIVE, Restriction.NEGATIVE):
+            jac_args = extract_modified_arguments(replace_expr, argnumber=1, do_index=False)
+            for jac_arg in jac_args:
+                # 5.1) Cut the expression to this ansatz function
                 replacement = {ma.expr: Zero(shape=ma.expr.ufl_shape,
                                              free_indices=ma.expr.ufl_free_indices,
                                              index_dimensions=ma.expr.ufl_index_dimensions)
-                               if ma.restriction != restriction else ma.expr
                                for ma in jac_args}
+                replacement[jac_arg.expr] = jac_arg.expr
                 jac_expr = replace_expression(replace_expr, replacemap=replacement)
 
-                if not isinstance(jac_expr, Zero):
-                    ret.append(AccumulationTerm(jac_expr, test_arg, indexmap, newi))
+                # 5.2) Propagate indexed zeros to simplify expression
+                jac_expr = zero_propagation(jac_expr)
+
+                indexed_jac_args = extract_modified_arguments(jac_expr, argnumber=1, do_index=True)
+                for restriction in (Restriction.NONE, Restriction.POSITIVE, Restriction.NEGATIVE):
+                    replacement = {ma.expr: Zero(shape=ma.expr.ufl_shape,
+                                                 free_indices=ma.expr.ufl_free_indices,
+                                                 index_dimensions=ma.expr.ufl_index_dimensions)
+                                   if ma.restriction != restriction else ma.expr
+                                   for ma in indexed_jac_args}
+
+                    jac_expr = replace_expression(jac_expr, replacemap=replacement)
+
+                    if not isinstance(jac_expr, Zero):
+                        ret.append(AccumulationTerm(jac_expr, test_arg, indexmap, newi))
         else:
             if not isinstance(replace_expr, Zero):
                 ret.append(AccumulationTerm(replace_expr, test_arg, indexmap, newi))