BUG: Cross-mesh interpolation from an H(div) space on a manifold

[pbrubeck]

Describe the bug
We are getting a segfault. However it is unclear to me if we can safely implement this operation, as degrees the of freedom mean different things on each side of a facet (the edge normals are not necesarily aligned on a manifold mesh).

Steps to Reproduce
Steps to reproduce the behavior:

from firedrake import *
m0 = UnitIcosahedralSphereMesh(0)
m1 = UnitIcosahedralSphereMesh(1)
m0.init_cell_orientations(SpatialCoordinate(m0))
m1.init_cell_orientations(SpatialCoordinate(m1))

V0 = FunctionSpace(m0, "RT", 1)
V1 = FunctionSpace(m1, "RT", 1)
u0 = Function(V0)

assemble(interpolate(u0, V1))

On main/release this will throw another error telling you that the VertexOnlyMesh does not have cell_orientations, this is fixed by replacing target_mesh -> source_mesh in

firedrake/firedrake/interpolation.py

Line 1029 in f6f7f7d

co = source_mesh.cell_orientations()

[pbrubeck]

Actually, the issue seems to be with interpolation from an H(div) space on a manifold

from firedrake import *
m0 = UnitIcosahedralSphereMesh(0)
m1 = UnitIcosahedralSphereMesh(1)
m0.init_cell_orientations(SpatialCoordinate(m0))
m1.init_cell_orientations(SpatialCoordinate(m1))

V0 = FunctionSpace(m0, "RT", 1)
V1 = VectorFunctionSpace(m1, "DG", 1)
u0 = Function(V0)

assemble(interpolate(u0, V1))

[pbrubeck]

H(curl) on the source and/or target mesh seems to work fine (at least no segfault)

[connorjward]

I'm a bit confused here. Is this supposed to be solved by #5052?

[pbrubeck]

I'm a bit confused here. Is this supposed to be solved by #5052?

The issue is very similar to #5052 which only deals with MeshHiearchy. Here the meshes are unrelated.

[pbrubeck]

The culprit seems to be

firedrake/tsfc/kernel_interface/common.py

Lines 75 to 79 in f8b7a56

	return gem.Conditional(gem.Comparison("==", co_int, gem.Literal(1)),
	gem.Literal(-1),
	gem.Conditional(gem.Comparison("==", co_int, gem.Zero()),
	gem.Literal(1),
	gem.Literal(numpy.nan)))

And the segfault goes away by replacing these lines with

return gem.Literal(1)

So probably we are not expecting the arguments to the interpolation kernel to be of int type somewhere deep down in the code generation stack.