Have `Function.at` make a vertex-only mesh at the specified points

[ReubenHill]

What we want

VertexOnlyMesh has a much better point location algorithm than Function.at. Without this algorithm, .at regularly produces a runtime error

firedrake/firedrake/function.py

Line 622 in 2a388c3

raise RuntimeError("Point evaluation gave different results across processes.")

which can only be fixed by reducing the mesh tolerance parameter.

The internals of .at should be replaced by code which creates a VertexOnlyMesh at the specified points and returns the values found by interpolating onto a P0DG function space at the point locations. A bunch of old point location code can then be removed.

What's blocking this

• VertexOnlyMesh cannot yet be created in an immersed manifold but as of Point evaluation on immersed manifold  #2879 , .at works for them. I'm working on fixing this in VertexOnlyMesh in immersed manifold and cross mesh interpolation from immersed manifold #3068 .
• Similarly, VertexOnlyMesh cannot yet be created in a high order mesh but .at will happily produce a value at a given point. Fixing this needs a more advanced voting algorithm. From the VertexOnlyMesh code:

  firedrake/firedrake/mesh.py

  Lines 2822 to 2826 in 2a388c3

  	# Bendy meshes require a smarter bounding box algorithm at partition and
  	# (especially) cell level. Projecting coordinates to Bernstein may be
  	# sufficient.
  	if np.any(np.asarray(mesh.coordinates.function_space().ufl_element().degree())) > 1:
  	raise NotImplementedError("Only straight edged meshes are supported")

Possible intermediate solution

We can either hold off on doing this until the above features have been added or we can implement this for cases where this does work and keep the old code for where it doesn't. Eventually we will have feature parity and can remove the old code.

[leo-collins]

Solved by #4543