diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py index cfff4f26d6cbd192571ae96d9d528948a9238b59..4f8da0d7bc715cbcd0a9a2a371ebc15c3d39bea5 100644 --- a/python/dune/perftool/options.py +++ b/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") diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py index d6adc144be19f513a6f27830a9fa5483d9ce84e7..71321cd7d7c6a46792614ff83c20d036ffccef8d 100644 --- a/python/dune/perftool/pdelab/basis.py +++ b/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