[WIP] Special case for Diagonal Matrix broken

python/dune/perftool/options.py

@@ -32,6 +32,7 @@ def get_form_compiler_arguments():
     parser.add_argument("--interactive", action="store_true", help="whether the optimization process should be guided interactively (also useful for debugging)")
     parser.add_argument("--print-transformations", action="store_true", help="print out dot files after ufl tree transformations")
     parser.add_argument("--print-transformations-dir", type=str, help="place where to put dot files (can be omitted)")
+    parser.add_argument("--diagonal-transformation-matrix", action="store_true", help="set option if the jacoby of the transformation is diagonal (axiparallel grids)")
 
     # These are the options that I deemed necessary in uflpdelab
 #     parser.add_argument("--param-class-file", type=str, help="The filename for the generated parameter class header")

python/dune/perftool/pdelab/basis.py

@@ -290,18 +290,40 @@ def evaluate_basis_gradient(leaf_element, name, restriction):
     jac = name_jacobian_inverse_transposed(restriction)
     index = lfs_iname(leaf_element, restriction, context='transformgrads')
     reference_gradients = name_reference_gradient(leaf_element, restriction)
-    instruction(inames=(index,
-                        quadrature_iname(),
-                        ),
-                code='{}.mv({}[{}][0], {}[{}]);'.format(jac,
-                                                        reference_gradients,
-                                                        index,
-                                                        name,
-                                                        index,
-                                                        ),
-                assignees=frozenset({name}),
-                read_variables=frozenset({reference_gradients}),
-                )
+
+    # Evaluate product of jacobianInvereseTransposed and gradient on reference element with loopy
+    dimindex_0 = dimension_iname(context='gradient_transformation', count=0)
+    from dune.perftool.options import get_option
+    if not get_option('diagonal-transformation-matrix'):
+        # dimindex_1 = dimension_iname(context='gradient_transformation', count=1)
+        reduction_expr = Product((
+            Subscript(Variable(jac), (Variable(dimindex_0), Variable(dimindex_0))),
+            Subscript(Variable(reference_gradients), (Variable(index), 0, Variable(dimindex_0)))))
+
+        assignee = Subscript(Variable(name), (Variable(index), Variable(dimindex_0)))
+
+        # Create exrpession instruction
+        forced_iname_deps = frozenset({quadrature_iname()}).union(frozenset({index}))
+        instruction(
+            expression=Reduction("sum", dimindex_0, reduction_expr, allow_simultaneous=True),
+            assignee=assignee,
+            forced_iname_deps=forced_iname_deps,
+            forced_iname_deps_is_final=True)
+    else:
+        dimindex_1 = dimension_iname(context='gradient_transformation', count=1)
+        reduction_expr = Product((
+            Subscript(Variable(jac), (Variable(dimindex_0), Variable(dimindex_1))),
+            Subscript(Variable(reference_gradients), (Variable(index), 0, Variable(dimindex_1)))))
+
+        assignee = Subscript(Variable(name), (Variable(index), Variable(dimindex_0)))
+
+        # Create exrpession instruction
+        forced_iname_deps = frozenset({quadrature_iname()}).union(frozenset({index})).union(frozenset({dimindex_0}))
+        instruction(
+            expression=Reduction("sum", dimindex_1, reduction_expr, allow_simultaneous=True),
+            assignee=assignee,
+            forced_iname_deps=forced_iname_deps,
+            forced_iname_deps_is_final=True)
 
 
 @symbol