Interpolating Constant onto Real on ExtrudedMesh fails

[ReubenHill]

MFE:

from firedrake import *
m = ExtrudedMesh(UnitSquareMesh(1, 1), 1)
V = FunctionSpace(m, "Real", 0)
v = interpolate(Constant(1.0), V) # Failure

fails with traceback

Traceback (most recent call last):
  File "/Users/rwh10/firedrake/src/firedrake/failing", line 4, in <module>
    v = interpolate(Constant(1.0), V)
  File "PETSc/Log.pyx", line 115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "PETSc/Log.pyx", line 116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/Users/rwh10/firedrake/src/firedrake/firedrake/interpolation.py", line 61, in interpolate
    return Interpolator(expr, V, subset=subset, access=access).interpolate(ad_block_tag=ad_block_tag)
  File "/Users/rwh10/firedrake/src/firedrake/firedrake/interpolation.py", line 87, in __init__
    self.callable, arguments = make_interpolator(expr, V, subset, access)
  File "PETSc/Log.pyx", line 115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "PETSc/Log.pyx", line 116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/Users/rwh10/firedrake/src/firedrake/firedrake/interpolation.py", line 218, in make_interpolator
    loops.extend(_interpolator(V, tensor, expr, subset, arguments, access))
  File "<decorator-gen-28>", line 2, in _interpolator
  File "/Users/rwh10/firedrake/src/firedrake/firedrake/utils.py", line 74, in wrapper
    return f(*args, **kwargs)
  File "/Users/rwh10/firedrake/src/firedrake/firedrake/interpolation.py", line 235, in _interpolator
    to_element = create_element(V.ufl_element())
  File "/Users/rwh10/firedrake/src/tsfc/tsfc/finatinterface.py", line 284, in create_element
    finat_element, deps = _create_element(ufl_element,
  File "/Users/rwh10/firedrake/src/tsfc/tsfc/finatinterface.py", line 314, in _create_element
    finat_element, deps = convert(ufl_element, **kwargs)
  File "/usr/local/Cellar/python@3.10/3.10.8/Frameworks/Python.framework/Versions/3.10/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/Users/rwh10/firedrake/src/tsfc/tsfc/finatinterface.py", line 188, in convert_finiteelement
    return lmbda(cell, element.degree()), set()
  File "/Users/rwh10/firedrake/src/FInAT/finat/fiat_elements.py", line 432, in __init__
    super(DiscontinuousLagrange, self).__init__(FIAT.DiscontinuousLagrange(cell, degree))
  File "/Users/rwh10/firedrake/src/fiat/FIAT/discontinuous_lagrange.py", line 188, in DiscontinuousLagrange
    return P0.P0(ref_el)
  File "/Users/rwh10/firedrake/src/fiat/FIAT/P0.py", line 54, in __init__
    poly_set = polynomial_set.ONPolynomialSet(ref_el, 0)
  File "/Users/rwh10/firedrake/src/fiat/FIAT/polynomial_set.py", line 136, in __init__
    num_exp_functions = expansions.polynomial_dimension(ref_el, degree)
  File "/Users/rwh10/firedrake/src/fiat/FIAT/expansions.py", line 435, in polynomial_dimension
    raise ValueError("Unknown reference element type.")
ValueError: Unknown reference element type.

In this function

 420     def polynomial_dimension(ref_el, degree):
 421         """Returns the dimension of the space of polynomials of degree no
 422         greater than degree on the reference element."""
 423         if ref_el.get_shape() == reference_element.POINT:
 424             if degree > 0:
 425                 raise ValueError("Only degree zero polynomials supported on point elements.")
 426             return 1
 427         elif ref_el.get_shape() == reference_element.LINE:
 428             return max(0, degree + 1)
 429         elif ref_el.get_shape() == reference_element.TRIANGLE:
 430             return max((degree + 1) * (degree + 2) // 2, 0)
 431         elif ref_el.get_shape() == reference_element.TETRAHEDRON:
 432             return max(0, (degree + 1) * (degree + 2) * (degree + 3) // 6)
 433         else:
 434  ->         raise ValueError("Unknown reference element type.")

For some reason ref_el.get_shape() reports 99 but I'm not sure where that number is coming from. The number 99 does not change if I change to ExtrudedMesh(UnitSquareMesh(2, 3), 15) nor if I switch to Constant(2.0, domain=m). It looks like the conversion of the UFL element to a FInAT one is not working as expected. I wonder if a further MFE can be found which specifically does that?

[wence-]

All elements on extruded meshes need to be tensor product elements (or else special). You don't realise this because UFL does some magic for P, DQ, etc...:

In [11]: FunctionSpace(m, "Real", 0).ufl_element()
Out[11]: FiniteElement('Real', TensorProductCell(triangle, interval), 0)

In [12]: FunctionSpace(m, "P", 1).ufl_element()
Out[12]: TensorProductElement(FiniteElement('Lagrange', triangle, 1), FiniteElement('Lagrange', interval, 1), TensorProductCell(triangle, interval)) # Not Finiteelement('P', TensorProductCell(triangle, interval), 1)

However it doesn't do the magic for Real (see __new__ in ufl.FiniteElement).

So what happens is that tsfc just falls through to normal fiat conversions and produces a P0 elements with a tensor product cell (ref_el.shape == 99), and things break.

You might think the workaround should be:

from firedrake import *
m = ExtrudedMesh(UnitSquareMesh(1, 1), 1)
V = FunctionSpace(m, TensorProductElement(FiniteElement("Real", triangle, 0), FiniteElement("Real", interval, 0)))
v = interpolate(Constant(1.0), V) 

But this doesn't work because it defeat's Firedrake's detection of the element as a real element (so it's just treated as DG0). (This should probably be fixed, if so, then you can update the UFL magic), otherwise you'll need to handle this case in tsfc's finatinterface, perhaps like this?

diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py
index 0c8b211..b1b072d 100644
--- a/tsfc/finatinterface.py
+++ b/tsfc/finatinterface.py
@@ -122,6 +122,12 @@ def convert_finiteelement(element, **kwargs):
     if element.family() == "Real" and element.cell().cellname() in {"quadrilateral", "hexahedron"}:
         lmbda = None
         element = ufl.FiniteElement("DQ", element.cell(), 0)
+    elif element.family() == "Real" and isinstance(element.cell(), ufl.TensorProductCell):
+        element = ufl.TensorProductElement(
+            *(ufl.FiniteElement("Real", cell, 0) for cell in element.cell().sub_cells()),
+            cell=element.cell(),
+        )
+        return _create_element(element, **kwargs)
     if lmbda is None:
         if element.cell().cellname() == "quadrilateral":
             # Handle quadrilateral short names like RTCF and RTCE.