BUG: Immersed manifold point evaluation gets wrong result for N1curl, fails for N1div

[ReubenHill]

Describe the bug
Point evaluation fails to get the expected answer when using an immersed sphere with certain function spaces:

from firedrake import *
import numpy as np

m = UnitIcosahedralSphereMesh(refinement_level=2, name='immersedsphere')
m.init_cell_orientations(SpatialCoordinate(m))
V = FunctionSpace(m, 'N1curl', 2)  # also fails with N2curl
expr = 2 * SpatialCoordinate(m)
f = Function(V).interpolate(expr)
point = m.coordinates.dat.data_ro[0]  # array([-0.52573111, 0.85065081, 0.0]) - definitely on the manifold!
f.at(point)  # array([ 2.36328118e-01,  2.03143488e-01, -6.67673661e-16])
2*point  # array([-1.05146222,  1.70130162,  0.0])
assert np.allclose(f.at(point), 2*point)  # fails

I get the same result if I interpolate onto a vertex-only mesh (see #3068).

If I switch to N1div or N2div I get loopy.diagnostic.TypeInferenceFailure: name not known in type inference: cell_orientations when using f.at(point).

With vertex-only mesh interpolation, I raise this runtime error

firedrake/firedrake/mesh.py

Line 2313 in 35f2b3c

raise RuntimeError("No cell orientations found, did you forget to call init_cell_orientations?")

when target mesh cell orientations are requested

firedrake/firedrake/interpolation.py

Line 960 in 35f2b3c

co = target_mesh.cell_orientations()

despite them being set in the script.

Expected behavior
Either the correct answer or a NotImplementedError

Error message
See above.

Environment:

• OS: MacOS
• Python version: 3.11.3
• Output of firedrake-status

Firedrake Configuration:
    package_manager: True
    minimal_petsc: False
    mpicc: None
    mpicxx: None
    mpif90: None
    mpiexec: None
    disable_ssh: False
    honour_petsc_dir: False
    with_parmetis: False
    slepc: False
    packages: []
    honour_pythonpath: False
    opencascade: False
    tinyasm: False
    torch: False
    petsc_int_type: int32
    cache_dir: /Users/rwh10/firedrake/.cache
    complex: False
    remove_build_files: False
    with_blas: None
    netgen: False
Additions:
    None
Environment:
    PYTHONPATH: None
    PETSC_ARCH: None
    PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT               |master                        |47f6c37   |False     |
|PyOP2               |master                        |8e4ad930  |False     |
|fiat                |master                        |8c66270   |False     |
|firedrake           |master                        |35f2b3cce |False     |
|h5py                |firedrake                     |6cc4c912  |False     |
|libspatialindex     |master                        |4768bf3   |True      |
|libsupermesh        |master                        |b145b65   |False     |
|loopy               |main                          |8158afdb  |False     |
|petsc               |firedrake                     |9364cb008b|False     |
|pyadjoint           |master                        |42959ef   |False     |
|pytest-mpi          |main                          |a478bc8   |False     |
|tsfc                |master                        |6f72c9c   |False     |
|ufl                 |master                        |55c281d7  |False     |
---------------------------------------------------------------------------

• Output of firedrake-status for VertexOnlyMesh investigations:

Firedrake Configuration:
    package_manager: True
    minimal_petsc: False
    mpicc: None
    mpicxx: None
    mpif90: None
    mpiexec: None
    disable_ssh: False
    honour_petsc_dir: False
    with_parmetis: False
    slepc: False
    packages: []
    honour_pythonpath: False
    opencascade: False
    tinyasm: False
    torch: False
    petsc_int_type: int32
    cache_dir: /Users/rwh10/firedrake/.cache
    complex: False
    remove_build_files: False
    with_blas: None
    netgen: False
Additions:
    None
Environment:
    PYTHONPATH: None
    PETSC_ARCH: None
    PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT               |master                        |47f6c37   |False     |
|PyOP2               |master                        |8e4ad930  |False     |
|fiat                |master                        |8c66270   |False     |
|firedrake           |ReubenHill/vom-immersed-manifold|455129ff2 |False     |
|h5py                |firedrake                     |6cc4c912  |False     |
|libspatialindex     |master                        |4768bf3   |True      |
|libsupermesh        |master                        |b145b65   |False     |
|loopy               |main                          |8158afdb  |False     |
|petsc               |firedrake                     |9364cb008b|False     |
|pyadjoint           |master                        |42959ef   |False     |
|pytest-mpi          |main                          |a478bc8   |False     |
|tsfc                |master                        |6f72c9c   |False     |
|ufl                 |master                        |55c281d7  |False     |
---------------------------------------------------------------------------

[ReubenHill]

Also get the loopy error for

V = FunctionSpace(m, 'BDM', 2)

and get the wrong answer for Regge:

from firedrake import *
import numpy as np

m = UnitIcosahedralSphereMesh(refinement_level=2, name='immersedsphere')
m.init_cell_orientations(SpatialCoordinate(m))
V = FunctionSpace(m, 'Regge', 2)
x = SpatialCoordinate(m)
expr = outer(x, x)
f = Function(V).interpolate(expr)
point = m.coordinates.dat.data_ro[0]  # array([-0.52573111, 0.85065081, 0.0])
f.at(point)
# array([[ 1.39627448e-02,  1.20021295e-02, -1.68957223e-17],
#        [ 1.20021295e-02,  1.03168192e-02, -1.40488013e-17],
#        [-1.55534977e-17, -1.33695127e-17, -1.52575501e-17]])
np.outer(point, point)
# array([[ 0.2763932, -0.4472136, -0.       ],
#        [-0.4472136,  0.7236068,  0.       ],
#        [-0.       ,  0.       ,  0.       ]])
assert np.allclose(f.at(point), np.outer(point, point))  # fails

[ReubenHill]

Update: N1curl/N2curl result is not necessarily wrong. We need a complete polynomial space for point evaluation to give an the naively expected result. A plane in R^3 (at z=0) with N1curl (degree 2), N2curl (degree 1) BDM (degree 1), RT (degree 2) should all give us the answer we expect in the test above.

Other issues need to be debugged.