diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 9f75b3cebf075ab3dc6280f0d563bf9e039b9597..62da51b6024fcf8d4d28b5fc6bfce700fe27ffc2 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -48,6 +48,7 @@ def get_form_compiler_arguments():
     parser.add_argument("--sumfact", action="store_true", help="Use sumfactorization")
     parser.add_argument("--vectorize-quad", action="store_true", help="whether to generate code with explicit vectorization")
     parser.add_argument("--vectorize-grads", action="store_true", help="whether to generate code with explicit vectorization")
+    parser.add_argument("--diagonal_jacobian_bug_workaround", action="store_true", help="Do not use diagonal_jacobian transformation and cast result of jacobianInverseTransposed to FieldMatrx. This is a temporary workaround to avoid a nasty splitting bug.")
 
     # Modify the positional argument to not be a list
     args = vars(parser.parse_args())
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 9e98245c3410528780dcf31702dade76751086d2..1b9669f48b989b2563530b5972ca19392a1f4075 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -316,10 +316,16 @@ def define_constant_jacobian_inveser_transposed(name, restriction):
 
     globalarg(name, dtype=np.float64, shape=(dim, dim), managed=False)
 
-    return 'auto {} = {}.jacobianInverseTransposed({});'.format(name,
-                                                                geo,
-                                                                pos,
-                                                                )
+    if get_option('diagonal_jacobian_bug_workaround'):
+        jit_type = 'Dune::FieldMatrix<double,{},{}>'.format(dim,dim)
+    else:
+        jit_type = 'auto'
+
+    return '{} {} = {}.jacobianInverseTransposed({});'.format(jit_type,
+                                                              name,
+                                                              geo,
+                                                              pos,
+                                                              )
 
 
 @backend(interface="define_jit", name="default")
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 088fa7425ffadffbf549ed8caf0fb79336de7bab..a925c2134a2b9a10d17346e7e6493b4e7dae4dff 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -376,8 +376,9 @@ def visit_integrals(integrals):
 
         # Maybe make the jacobian inverse diagonal!
         if get_option('diagonal_transformation_matrix'):
-            from dune.perftool.ufl.transformations.axiparallel import diagonal_jacobian
-            integrand = diagonal_jacobian(integrand)
+            if not get_option('diagonal_jacobian_bug_workaround'):
+                from dune.perftool.ufl.transformations.axiparallel import diagonal_jacobian
+                integrand = diagonal_jacobian(integrand)
 
         # Gather dimension indices
         from dune.perftool.ufl.dimensionindex import dimension_index_mapping
diff --git a/test/stokes/stokes_quadrilateral.mini b/test/stokes/stokes_quadrilateral.mini
index 0c2a459d2ebdec9279635e47c679d13aff548c89..a35a389f3fae9642065498c140ef0ec581339968 100644
--- a/test/stokes/stokes_quadrilateral.mini
+++ b/test/stokes/stokes_quadrilateral.mini
@@ -14,3 +14,4 @@ extension = vtu
 numerical_jacobian = 1, 0 | expand num
 exact_solution_expression = g
 compare_l2errorsquared = 1e-11
+diagonal_jacobian_bug_workaround = 1