Make turning off diagonal jacobian through option possible

python/dune/perftool/options.py

@@ -48,7 +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.")
+    parser.add_argument("--turn-off-diagonal-jacobian", action="store_true", help="Do not use diagonal_jacobian transformation on the ufl tree and cast result of jacobianInverseTransposed into a FieldMatrix.")
 
     # Modify the positional argument to not be a list
     args = vars(parser.parse_args())

python/dune/perftool/pdelab/geometry.py

@@ -288,8 +288,12 @@ def name_unit_inner_normal():
 
 
 def type_jacobian_inverse_transposed(restriction):
-    geo = type_cell_geometry(restriction)
-    return "typename {}::JacobianInverseTransposed".format(geo)
+    if get_option('turn_off_diagonal_jacobian'):
+        dim = world_dimension()
+        return "typename Dune::FieldMatrix<double,{},{}>".format(dim,dim)
+    else:
+        geo = type_cell_geometry(restriction)
+        return "typename {}::JacobianInverseTransposed".format(geo)
 
 
 @kernel_cached
@@ -316,7 +320,7 @@ def define_constant_jacobian_inveser_transposed(name, restriction):
 
     globalarg(name, dtype=np.float64, shape=(dim, dim), managed=False)
 
-    if get_option('diagonal_jacobian_bug_workaround'):
+    if get_option('turn_off_diagonal_jacobian'):
         jit_type = 'Dune::FieldMatrix<double,{},{}>'.format(dim,dim)
     else:
         jit_type = 'auto'

python/dune/perftool/pdelab/localoperator.py

@@ -376,7 +376,7 @@ def visit_integrals(integrals):
 
         # Maybe make the jacobian inverse diagonal!
         if get_option('diagonal_transformation_matrix'):
-            if not get_option('diagonal_jacobian_bug_workaround'):
+            if not get_option('turn_off_diagonal_jacobian'):
                 from dune.perftool.ufl.transformations.axiparallel import diagonal_jacobian
                 integrand = diagonal_jacobian(integrand)
 

test/stokes/stokes_quadrilateral.mini

@@ -14,4 +14,3 @@ extension = vtu
 numerical_jacobian = 1, 0 | expand num
 exact_solution_expression = g
 compare_l2errorsquared = 1e-11
-diagonal_jacobian_bug_workaround = 1