Degree estimation doesn't cater for non-affine cell geometry

[Martin Sandve Alnæs]

For example 1.0*dx(quadratic_mesh) needs higher degree than 1 to be exact because of the |detJ| factor.

[Lawrence Mitchell]

Is the way to fix this to modify the estimated degree when applying integrand scaling?

[Martin Sandve Alnæs]

Maybe, but how? Adding some function of mesh degree?

Does it seem reasonable to say detJ degree is roughly like (x degree - 1)*gdim, which is 0 for affine and gdim for quadratic mesh? Or is that overkill?

[Lawrence Mitchell]

So I would apply geometry lowering to the integrand scaling, estimate the degree of the resulting expression and then add the resulting estimation to the estimated degree of the whole integral.

[Lawrence Mitchell]

Of course the problem with that is that the scaling is non polynomial of the cell is not affine due to the abs

[Miklós Homolya]

I don't think that det(J) changes sign inside a cell. AFAIK that would mean in non-immersed that a cell self-overlaps.

[Prof Garth Wells]

Indeed, the sign of det(J) will not change for a cell that isn't screwed up.

[Martin Sandve Alnæs]

Yes, abs can be ignored. Won't lowering then estimation lead to similar overkill degrees that we saw before?

[Prof Garth Wells]

What if the lower bound on the integration degree is equal to the element mapping degree?

[Lawrence Mitchell]

Martin, I had anticipated something like:

1. estimate degrees for integrand
2. Separately estimate degree for integrand scaling (this should be exact because the detJ and similar terms are polynomial)
3. Sum the estimated degrees from 1 and 2 giving a final estimated degree.

Applying geometry lowering before degree estimation just to the integrand scaling (so that the degree estimator doesn't need to learn about too many new things) should not blow up, I think.

[Lawrence Mitchell]

That would look something like:

diff --git a/ufl/algorithms/apply_integral_scaling.py b/ufl/algorithms/apply_integral_scaling.py
index f36b02e..607b629 100644
--- a/ufl/algorithms/apply_integral_scaling.py
+++ b/ufl/algorithms/apply_integral_scaling.py
@@ -24,6 +24,8 @@ from ufl.log import error, warning
 from ufl.assertions import ufl_assert

 from ufl.classes import JacobianDeterminant, FacetJacobianDeterminant, QuadratureWeight, Form, Integral
+from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering
+from ufl.algorithms.estimate_degrees import estimate_total_polynomial_degree


 def compute_integrand_scaling_factor(integral):
@@ -36,13 +38,21 @@ def compute_integrand_scaling_factor(integral):
     tdim = domain.topological_dimension()
     #gdim = domain.geometric_dimension()

+    # Polynomial degree of integrand scaling
+    degree = 0
     if integral_type == "cell":
-        scale = abs(JacobianDeterminant(domain)) * weight
+        detJ = JacobianDeterminant(domain)
+        degree = estimate_total_polynomial_degree(apply_geometry_lowering(detJ))
+        # Despite the abs, |detJ| is polynomial except for
+        # self-intersecting cells, where we have other problems.
+        scale = abs(detJ) * weight

     elif integral_type.startswith("exterior_facet"):
         if tdim > 1:
             # Scaling integral by facet jacobian determinant and quadrature weight
-            scale = FacetJacobianDeterminant(domain) * weight
+            detFJ = FacetJacobianDeterminant(domain)
+            degree = estimate_total_polynomial_degree(apply_geometry_lowering(detFJ))
+            scale = detFJ * weight
         else:
             # No need to scale 'integral' over a vertex
             scale = 1
@@ -50,7 +60,9 @@ def compute_integrand_scaling_factor(integral):
     elif integral_type.startswith("interior_facet"):
         if tdim > 1:
             # Scaling integral by facet jacobian determinant from one side and quadrature weight
-            scale = FacetJacobianDeterminant(domain)('+') * weight
+            detFJ = FacetJacobianDeterminant(domain)
+            degree = estimate_total_polynomial_degree(apply_geometry_lowering(detFJ))
+            scale = detFJ('+') * weight
         else:
             # No need to scale 'integral' over a vertex
             scale = 1
@@ -66,7 +78,7 @@ def compute_integrand_scaling_factor(integral):
     else:
         error("Unknown integral type {}, don't know how to scale.".format(integral_type))

-    return scale
+    return scale, degree


 def apply_integral_scaling(form):
@@ -81,9 +93,12 @@ def apply_integral_scaling(form):
     elif isinstance(form, Integral):
         integral = form
         # Compute and apply integration scaling factor
-        scale = compute_integrand_scaling_factor(integral)
+        scale, degree = compute_integrand_scaling_factor(integral)
+        md = {}
+        md.update(integral.metadata())
+        md["estimated_polynomial_degree"] += degree
         newintegrand = integral.integrand() * scale
-        return integral.reconstruct(integrand=newintegrand)
+        return integral.reconstruct(integrand=newintegrand, metadata=md)

     else:
         error("Invalid type %s" % (form.__class__.__name__,))

[Martin Sandve Alnæs]

I won't have time to look at this until August or September when I'm back from paternity leave.

Your code looks reasonable to me, feel free to make a PR and ping other FEniCS devs for merging, should be fine as long as it doesn't break FFC regression tests with uflacs.

[Miklós Homolya]

I suppose this is now fixed by pull request #109.