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 https://github.com/firedrakeproject/firedrake/blob/35f2b3cce725d0d6a3ed277ea8bda3a287e00476/firedrake/mesh.py#L2313 when target mesh cell orientations are requested https://github.com/firedrakeproject/firedrake/blob/35f2b3cce725d0d6a3ed277ea8bda3a287e00476/firedrake/interpolation.py#L960 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.