tsfc: cancel Jacobian products and reciprocal factors

.github/actions/install/action.yml

@@ -155,6 +155,8 @@ runs:
           --extra-index-url https://download.pytorch.org/whl/cpu \
           "./firedrake-repo[${{ inputs.deps }}]"
 
+        pip install -v --no-deps --ignore-installed git+https://github.com/firedrakeproject/fiat.git@pbrubeck/cancel-jacobian
+        pip install -v --no-deps --ignore-installed git+https://github.com/firedrakeproject/ufl.git@pbrubeck/cancel-jacobian
         firedrake-clean
         pip list
 

tsfc/ufl_utils.py

@@ -14,6 +14,8 @@
 from ufl.algorithms.apply_derivatives import apply_derivatives
 from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering
 from ufl.algorithms.apply_restrictions import apply_restrictions
+from ufl.algorithms.cancel_jacobian_products import cancel_jacobian_products
+from ufl.algorithms.remove_component_tensors import remove_component_tensors
 from ufl.algorithms.comparison_checker import do_comparison_check
 from ufl.algorithms.remove_complex_nodes import remove_complex_nodes
 from ufl.algorithms.signature import compute_expression_signature
@@ -40,6 +42,7 @@ def compute_form_data(form,
                       do_apply_integral_scaling=True,
                       do_apply_geometry_lowering=True,
                       preserve_geometry_types=preserve_geometry_types,
+                      do_cancel_jacobian_products=True,
                       do_apply_default_restrictions=True,
                       do_apply_restrictions=True,
                       do_estimate_degrees=True,
@@ -62,6 +65,7 @@ def compute_form_data(form,
         do_apply_integral_scaling=do_apply_integral_scaling,
         do_apply_geometry_lowering=do_apply_geometry_lowering,
         preserve_geometry_types=preserve_geometry_types,
+        do_cancel_jacobian_products=do_cancel_jacobian_products,
         do_apply_default_restrictions=do_apply_default_restrictions,
         do_apply_restrictions=do_apply_restrictions,
         do_estimate_degrees=do_estimate_degrees,
@@ -136,11 +140,19 @@ def preprocess_expression(expression, complex_mode=False,
         expression = do_comparison_check(expression)
     else:
         expression = remove_complex_nodes(expression)
+    jacobian_types = (Jacobian, JacobianInverse, JacobianDeterminant)
+    lowering_preserve_types = preserve_geometry_types + jacobian_types
     expression = apply_algebra_lowering(expression)
     expression = apply_derivatives(expression)
     expression = apply_function_pullbacks(expression)
-    expression = apply_geometry_lowering(expression, preserve_geometry_types)
+    expression = apply_geometry_lowering(expression, lowering_preserve_types)
+    expression = apply_derivatives(expression)
+    expression = apply_geometry_lowering(expression, lowering_preserve_types)
     expression = apply_derivatives(expression)
+    # Cancel contractions of the Jacobian with its inverse before
+    # expanding the inverse into individual matrix entries
+    expression = remove_component_tensors(expression)
+    expression = cancel_jacobian_products(expression)
     expression = apply_geometry_lowering(expression, preserve_geometry_types)
     expression = apply_derivatives(expression)
     if not complex_mode: