From 702449a09467936c50ddbe9fb5163a60546b2d3f Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 13:04:36 +0100
Subject: [PATCH 01/35] first attempt

---
 firedrake/function.py | 39 ++++++++++++++++++++++++++++++++++++++-
 1 file changed, 38 insertions(+), 1 deletion(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 4cce6b9edb..74cb7cae43 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -18,7 +18,7 @@
 from finat.ufl import MixedElement
 from firedrake.utils import ScalarType, IntType, as_ctypes
 
-from firedrake import functionspaceimpl
+from firedrake import functionspaceimpl, VertexOnlyMesh, FunctionSpace, assemble, interpolate
 from firedrake.cofunction import Cofunction, RieszMap
 from firedrake import utils
 from firedrake.adjoint_utils import FunctionMixin
@@ -698,6 +698,43 @@ def __str__(self):
         return "domain %s does not contain point %s" % (self.domain, self.point)
 
 
+class PointEvaluator:
+    r"""Convenience class for evaluating a :class:`Function` at a set of points."""
+
+    def __init__(self, mesh, points, tolerance=None, ignore_missing_points=False, input_ordered=True):
+        r"""
+        :arg mesh: The :class:`Mesh` on which the :class:`Function` is defined.
+        :arg points: Array of points to evaluate the function at.
+        :kwarg tolerance: Tolerance to use when checking if a point is
+            in a cell. Default is the ``tolerance`` of the :class:`Mesh`.
+        :kwarg ignore_missing_points: If ``True``, do not raise an error if a point is not found.
+        :kwarg input_ordered: If ``True``, return results in the order of the input points.
+        """
+        self.mesh = mesh
+        self.points = np.asarray(points, dtype=utils.ScalarType)
+        self.tolerance = tolerance or mesh.tolerance
+        self.ignore_missing_points = ignore_missing_points
+        self.input_ordered = input_ordered
+        self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour="warn", redundant=False, tolerance=tolerance)
+
+    def evaluate(self, function):
+        r"""Evaluate the given :class:`Function` at the points provided to this
+        :class:`PointEvaluator`.
+
+        :arg function: The :class:`Function` to evaluate.
+        :returns: A NumPy array of values at the points.
+        """
+        if not isinstance(function, Function):
+            raise TypeError(f"Expected a Function, got f{type(function).__name__}")
+
+        P0DG = FunctionSpace(self.vom, "DG", 0)
+        f_at_points = assemble(interpolate(function, P0DG))
+        if self.input_ordered:
+            P0DG_io = FunctionSpace(self.vom.input_ordering, "DG", 0)
+            f_at_points = assemble(interpolate(f_at_points, P0DG_io))
+        return f_at_points.dat.data_ro
+        
+
 @PETSc.Log.EventDecorator()
 def make_c_evaluate(function, c_name="evaluate", ldargs=None, tolerance=None):
     r"""Generates, compiles and loads a C function to evaluate the

From a54a56d354fcdd110ba17d5627532ec590591063 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 13:14:36 +0100
Subject: [PATCH 02/35] fix imports

---
 firedrake/function.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 74cb7cae43..94ea8c97cf 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -18,14 +18,14 @@
 from finat.ufl import MixedElement
 from firedrake.utils import ScalarType, IntType, as_ctypes
 
-from firedrake import functionspaceimpl, VertexOnlyMesh, FunctionSpace, assemble, interpolate
+from firedrake import functionspaceimpl
 from firedrake.cofunction import Cofunction, RieszMap
 from firedrake import utils
 from firedrake.adjoint_utils import FunctionMixin
 from firedrake.petsc import PETSc
 
 
-__all__ = ['Function', 'PointNotInDomainError', 'CoordinatelessFunction']
+__all__ = ['Function', 'PointNotInDomainError', 'CoordinatelessFunction', 'PointEvaluator']
 
 
 class _CFunction(ctypes.Structure):
@@ -710,6 +710,7 @@ def __init__(self, mesh, points, tolerance=None, ignore_missing_points=False, in
         :kwarg ignore_missing_points: If ``True``, do not raise an error if a point is not found.
         :kwarg input_ordered: If ``True``, return results in the order of the input points.
         """
+        from firedrake.mesh import VertexOnlyMesh
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
         self.tolerance = tolerance or mesh.tolerance
@@ -726,6 +727,7 @@ def evaluate(self, function):
         """
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got f{type(function).__name__}")
+        from firedrake import interpolate, assemble, FunctionSpace
 
         P0DG = FunctionSpace(self.vom, "DG", 0)
         f_at_points = assemble(interpolate(function, P0DG))
@@ -733,7 +735,7 @@ def evaluate(self, function):
             P0DG_io = FunctionSpace(self.vom.input_ordering, "DG", 0)
             f_at_points = assemble(interpolate(f_at_points, P0DG_io))
         return f_at_points.dat.data_ro
-        
+
 
 @PETSc.Log.EventDecorator()
 def make_c_evaluate(function, c_name="evaluate", ldargs=None, tolerance=None):

From 171700844d6d2a25e9ddcfb5c51c462019006b4b Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 13:21:14 +0100
Subject: [PATCH 03/35] fix missing_points_behaviour

---
 firedrake/function.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 94ea8c97cf..dbbea062a6 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -701,22 +701,24 @@ def __str__(self):
 class PointEvaluator:
     r"""Convenience class for evaluating a :class:`Function` at a set of points."""
 
-    def __init__(self, mesh, points, tolerance=None, ignore_missing_points=False, input_ordered=True):
+    def __init__(self, mesh, points, tolerance=None, missing_points_behaviour: str = "error", input_ordered=True):
         r"""
         :arg mesh: The :class:`Mesh` on which the :class:`Function` is defined.
         :arg points: Array of points to evaluate the function at.
         :kwarg tolerance: Tolerance to use when checking if a point is
             in a cell. Default is the ``tolerance`` of the :class:`Mesh`.
-        :kwarg ignore_missing_points: If ``True``, do not raise an error if a point is not found.
+        :kwarg missing_points_behaviour: Behaviour when a point is not found in the mesh. Options are:
+        "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh. 
+        "warn": warn if a point is not found in the mesh, but continue.
+        "ignore": ignore points not found in the mesh.
         :kwarg input_ordered: If ``True``, return results in the order of the input points.
         """
         from firedrake.mesh import VertexOnlyMesh
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
         self.tolerance = tolerance or mesh.tolerance
-        self.ignore_missing_points = ignore_missing_points
         self.input_ordered = input_ordered
-        self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour="warn", redundant=False, tolerance=tolerance)
+        self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour=missing_points_behaviour, redundant=False, tolerance=tolerance)
 
     def evaluate(self, function):
         r"""Evaluate the given :class:`Function` at the points provided to this

From aaf0f9387f9ecac08b00ff8c33b8aea84fbf8bb1 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 13:22:22 +0100
Subject: [PATCH 04/35] move input_ordered to evaluate method

---
 firedrake/function.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index dbbea062a6..8f072cd315 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -701,7 +701,7 @@ def __str__(self):
 class PointEvaluator:
     r"""Convenience class for evaluating a :class:`Function` at a set of points."""
 
-    def __init__(self, mesh, points, tolerance=None, missing_points_behaviour: str = "error", input_ordered=True):
+    def __init__(self, mesh, points, tolerance=None, missing_points_behaviour: str = "error"):
         r"""
         :arg mesh: The :class:`Mesh` on which the :class:`Function` is defined.
         :arg points: Array of points to evaluate the function at.
@@ -711,20 +711,19 @@ def __init__(self, mesh, points, tolerance=None, missing_points_behaviour: str =
         "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh. 
         "warn": warn if a point is not found in the mesh, but continue.
         "ignore": ignore points not found in the mesh.
-        :kwarg input_ordered: If ``True``, return results in the order of the input points.
         """
         from firedrake.mesh import VertexOnlyMesh
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
         self.tolerance = tolerance or mesh.tolerance
-        self.input_ordered = input_ordered
         self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour=missing_points_behaviour, redundant=False, tolerance=tolerance)
 
-    def evaluate(self, function):
+    def evaluate(self, function: Function, input_ordered: bool = True):
         r"""Evaluate the given :class:`Function` at the points provided to this
         :class:`PointEvaluator`.
 
         :arg function: The :class:`Function` to evaluate.
+        :kwarg input_ordered: If ``True``, return results in the order of the input points.
         :returns: A NumPy array of values at the points.
         """
         if not isinstance(function, Function):
@@ -733,7 +732,7 @@ def evaluate(self, function):
 
         P0DG = FunctionSpace(self.vom, "DG", 0)
         f_at_points = assemble(interpolate(function, P0DG))
-        if self.input_ordered:
+        if input_ordered:
             P0DG_io = FunctionSpace(self.vom.input_ordering, "DG", 0)
             f_at_points = assemble(interpolate(f_at_points, P0DG_io))
         return f_at_points.dat.data_ro

From 741265e462d8ddcf65fb3b8cd0382bb5408ea5e5 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 13:58:30 +0100
Subject: [PATCH 05/35] fix missing points behaviour

---
 firedrake/function.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 8f072cd315..46fcb371aa 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -723,7 +723,8 @@ def evaluate(self, function: Function, input_ordered: bool = True):
         :class:`PointEvaluator`.
 
         :arg function: The :class:`Function` to evaluate.
-        :kwarg input_ordered: If ``True``, return results in the order of the input points.
+        :kwarg input_ordered: If ``True``, return results in the order of the input points. If any 
+        points were not found in the mesh, they will be return as np.nan.
         :returns: A NumPy array of values at the points.
         """
         if not isinstance(function, Function):
@@ -734,8 +735,11 @@ def evaluate(self, function: Function, input_ordered: bool = True):
         f_at_points = assemble(interpolate(function, P0DG))
         if input_ordered:
             P0DG_io = FunctionSpace(self.vom.input_ordering, "DG", 0)
-            f_at_points = assemble(interpolate(f_at_points, P0DG_io))
-        return f_at_points.dat.data_ro
+            f_at_points_io = Function(P0DG_io).assign(np.nan)
+            f_at_points_io.interpolate(f_at_points)
+            return f_at_points_io.dat.data_ro
+        else:
+            return f_at_points.dat.data_ro
 
 
 @PETSc.Log.EventDecorator()

From add7af003c9481e709cb53661606f109ff6a07f8 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 15:01:17 +0100
Subject: [PATCH 06/35] deal with vector and tensor and mixed functions

---
 firedrake/function.py | 27 ++++++++++++++++++++++-----
 1 file changed, 22 insertions(+), 5 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 46fcb371aa..7241cb10b0 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -11,6 +11,7 @@
 from collections.abc import Collection
 from numbers import Number
 from pathlib import Path
+from functools import partial
 
 from pyop2 import op2, mpi
 from pyop2.exceptions import DataTypeError, DataValueError
@@ -725,16 +726,32 @@ def evaluate(self, function: Function, input_ordered: bool = True):
         :arg function: The :class:`Function` to evaluate.
         :kwarg input_ordered: If ``True``, return results in the order of the input points. If any 
         points were not found in the mesh, they will be return as np.nan.
-        :returns: A NumPy array of values at the points.
+        :returns: A Numpy array of values at the points. If the function is a mixed function, a list
+        of Numpy arrays is returned, one for each subfunction.
         """
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got f{type(function).__name__}")
-        from firedrake import interpolate, assemble, FunctionSpace
-
-        P0DG = FunctionSpace(self.vom, "DG", 0)
+        from firedrake import interpolate, assemble, FunctionSpace, VectorFunctionSpace, TensorFunctionSpace
+
+        subfunctions = function.subfunctions
+        if len(subfunctions) > 1:
+            result = []
+            for subfunction in subfunctions:
+                result.append(self.evaluate(subfunction, input_ordered=input_ordered))
+            return result
+
+        shape = function.ufl_function_space().value_shape
+        if len(shape) == 0:
+            fs = FunctionSpace
+        elif len(shape) == 1:
+            fs = VectorFunctionSpace
+        else:
+            fs = partial(TensorFunctionSpace, shape=shape)
+        P0DG = fs(self.vom, "DG", 0)
         f_at_points = assemble(interpolate(function, P0DG))
+
         if input_ordered:
-            P0DG_io = FunctionSpace(self.vom.input_ordering, "DG", 0)
+            P0DG_io = fs(self.vom.input_ordering, "DG", 0)
             f_at_points_io = Function(P0DG_io).assign(np.nan)
             f_at_points_io.interpolate(f_at_points)
             return f_at_points_io.dat.data_ro

From e554eeadd9ab408acbd7172d9cd13a6cc01b4638 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 15:19:30 +0100
Subject: [PATCH 07/35] docstrings and type hints

fix

fix

fix
---
 firedrake/function.py | 49 ++++++++++++++++++++++++++++---------------
 1 file changed, 32 insertions(+), 17 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 7241cb10b0..e6488874df 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -24,7 +24,7 @@
 from firedrake import utils
 from firedrake.adjoint_utils import FunctionMixin
 from firedrake.petsc import PETSc
-
+from firedrake.mesh import MeshGeometry
 
 __all__ = ['Function', 'PointNotInDomainError', 'CoordinatelessFunction', 'PointEvaluator']
 
@@ -702,32 +702,47 @@ def __str__(self):
 class PointEvaluator:
     r"""Convenience class for evaluating a :class:`Function` at a set of points."""
 
-    def __init__(self, mesh, points, tolerance=None, missing_points_behaviour: str = "error"):
+    def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None, 
+                 missing_points_behaviour: str | None = "error") -> None:
         r"""
-        :arg mesh: The :class:`Mesh` on which the :class:`Function` is defined.
-        :arg points: Array of points to evaluate the function at.
-        :kwarg tolerance: Tolerance to use when checking if a point is
-            in a cell. Default is the ``tolerance`` of the :class:`Mesh`.
-        :kwarg missing_points_behaviour: Behaviour when a point is not found in the mesh. Options are:
-        "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh. 
-        "warn": warn if a point is not found in the mesh, but continue.
-        "ignore": ignore points not found in the mesh.
+        Parameters
+        ----------
+        mesh : :class:`MeshGeometry`
+            The :class:`Mesh` on which the :class:`Function` is defined.
+        points : np.ndarray | list
+            Array of points to evaluate the function at.
+        tolerance : float | None
+            Tolerance to use when checking if a point is in a cell. Default is the ``tolerance`` of the :class:`Mesh`.
+        missing_points_behaviour : str | None
+            Behaviour when a point is not found in the mesh. Options are:
+            "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh. 
+            "warn": warn if a point is not found in the mesh, but continue.
+            None: ignore points not found in the mesh.
         """
         from firedrake.mesh import VertexOnlyMesh
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
-        self.tolerance = tolerance or mesh.tolerance
         self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour=missing_points_behaviour, redundant=False, tolerance=tolerance)
 
-    def evaluate(self, function: Function, input_ordered: bool = True):
+    def evaluate(self, function: Function, input_ordered: bool = True) -> np.ndarray | list[np.ndarray]:
         r"""Evaluate the given :class:`Function` at the points provided to this
         :class:`PointEvaluator`.
 
-        :arg function: The :class:`Function` to evaluate.
-        :kwarg input_ordered: If ``True``, return results in the order of the input points. If any 
-        points were not found in the mesh, they will be return as np.nan.
-        :returns: A Numpy array of values at the points. If the function is a mixed function, a list
-        of Numpy arrays is returned, one for each subfunction.
+        Parameters
+        ----------
+        function : :class:`Function`
+            The :class:`Function` to evaluate.
+        input_ordered : bool
+            If ``True``, return results in the order of the input points. If any 
+            points were not found in the mesh, they will be return as np.nan.
+
+        Returns
+        -------
+        np.ndarray | list[np.ndarray]
+            A Numpy array of values at the points. If the function is scalar-valued, the Numpy array
+            has shape (len(points),). If the function is vector-valued with shape (n,), the Numpy array has shape
+            (len(points), n). If the function is tensor-valued with shape (n, m), the Numpy array has shape
+            (len(points), n, m). If the function is a mixed function, a list of Numpy arrays is returned, one for each subfunction.
         """
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got f{type(function).__name__}")

From 7c603c652c9576173e1ca15b9e794271b0cf0810 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 16:11:33 +0100
Subject: [PATCH 08/35] add tests

---
 .../regression/test_point_eval_api.py         | 59 +++++++++++++++++++
 1 file changed, 59 insertions(+)

diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index e953a423db..c85397712b 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -161,3 +161,62 @@ def test_tolerance():
     with pytest.raises(PointNotInDomainError):
         f.at((-0.1, 0.4))
     assert np.allclose([-1e-12, 0.4], f.at((-1e-12, 0.4)))
+
+
+@pytest.fixture(scope="module")
+def mesh_and_points():
+    points = [[0.1, 0.1], [0.2, 0.2], [0.3, 0.3]]
+    mesh = UnitSquareMesh(10, 10)
+    evaluator = PointEvaluator(mesh, points)
+    return mesh, evaluator
+
+
+@pytest.mark.parallel([1, 3])
+def test_point_evaluator_scalar(mesh_and_points):
+    mesh, evaluator = mesh_and_points
+    points = evaluator.points
+    V = FunctionSpace(mesh, "CG", 1)
+    f = Function(V)
+    x, y = SpatialCoordinate(mesh)
+    f.interpolate(x + y)
+
+    # Test standard scalar function evaluation at points
+    eval = PointEvaluator(mesh, points)
+    f_at_points = eval.evaluate(f)
+    assert np.allclose(f_at_points, [0.2, 0.4, 0.6])
+
+    # Test standard scalar function with missing points
+    eval_missing = PointEvaluator(mesh, np.append(points, [[1.5, 1.5]], axis=0), missing_points_behaviour=None)
+    f_at_points_missing = eval_missing.evaluate(f, input_ordered=False)
+    assert not np.isnan(f_at_points_missing).any()
+    f_at_points_missing_io = eval_missing.evaluate(f)
+    assert np.isnan(f_at_points_missing_io[-1])
+
+@pytest.mark.parallel([1, 3])
+def test_point_evaluator_vector_tensor_mixed(mesh_and_points):
+    mesh, evaluator = mesh_and_points
+    V_vec = VectorFunctionSpace(mesh, "CG", 1)
+    f_vec = Function(V_vec)
+    x, y = SpatialCoordinate(mesh)
+    f_vec.interpolate(as_vector([x, y]))
+    f_vec_at_points = evaluator.evaluate(f_vec)
+    vec_expected = np.array([[0.1, 0.1], [0.2, 0.2], [0.3, 0.3]])
+    assert np.allclose(f_vec_at_points, vec_expected)
+
+    V_tensor = TensorFunctionSpace(mesh, "CG", 1, shape=(2, 3))
+    f_tensor = Function(V_tensor)
+    f_tensor.interpolate(as_matrix([[x, y, x*y], [y, x, x*y]]))
+    f_tensor_at_points = evaluator.evaluate(f_tensor)
+    tensor_expected = np.array([[[0.1, 0.1, 0.01], [0.1, 0.1, 0.01]], 
+                                [[0.2, 0.2, 0.04], [0.2, 0.2, 0.04]], 
+                                [[0.3, 0.3, 0.09], [0.3, 0.3, 0.09]]])
+    assert np.allclose(f_tensor_at_points, tensor_expected)
+
+    V_mixed = V_vec * V_tensor
+    f_mixed = Function(V_mixed)
+    f_vec, f_tensor = f_mixed.subfunctions
+    f_vec.interpolate(as_vector([x, y]))
+    f_tensor.interpolate(as_matrix([[x, y, x*y], [y, x, x*y]]))
+    f_mixed_at_points = evaluator.evaluate(f_mixed)
+    assert np.allclose(f_mixed_at_points[0], vec_expected)
+    assert np.allclose(f_mixed_at_points[1], tensor_expected)
\ No newline at end of file

From d04c04cc584d9499f7cfa7b0b03a3a1a6df72ddc Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 16:11:41 +0100
Subject: [PATCH 09/35] add redundant kwarg

---
 firedrake/function.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index e6488874df..d2baa645dc 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -703,7 +703,7 @@ class PointEvaluator:
     r"""Convenience class for evaluating a :class:`Function` at a set of points."""
 
     def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None, 
-                 missing_points_behaviour: str | None = "error") -> None:
+                 missing_points_behaviour: str | None = "error", redundant: bool = True) -> None:
         r"""
         Parameters
         ----------
@@ -722,7 +722,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
         from firedrake.mesh import VertexOnlyMesh
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
-        self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour=missing_points_behaviour, redundant=False, tolerance=tolerance)
+        self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour=missing_points_behaviour, redundant=redundant, tolerance=tolerance)
 
     def evaluate(self, function: Function, input_ordered: bool = True) -> np.ndarray | list[np.ndarray]:
         r"""Evaluate the given :class:`Function` at the points provided to this

From 2b77180dbec0506deb64417992bcd8003a232b47 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 21 Aug 2025 16:27:12 +0100
Subject: [PATCH 10/35] broadcast result to other ranks

fix
---
 firedrake/function.py | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index d2baa645dc..7e6ffbfcc9 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -722,9 +722,10 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
         from firedrake.mesh import VertexOnlyMesh
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
+        self.redundant = redundant
         self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour=missing_points_behaviour, redundant=redundant, tolerance=tolerance)
 
-    def evaluate(self, function: Function, input_ordered: bool = True) -> np.ndarray | list[np.ndarray]:
+    def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
         r"""Evaluate the given :class:`Function` at the points provided to this
         :class:`PointEvaluator`.
 
@@ -752,7 +753,7 @@ def evaluate(self, function: Function, input_ordered: bool = True) -> np.ndarray
         if len(subfunctions) > 1:
             result = []
             for subfunction in subfunctions:
-                result.append(self.evaluate(subfunction, input_ordered=input_ordered))
+                result.append(self.evaluate(subfunction))
             return result
 
         shape = function.ufl_function_space().value_shape
@@ -765,14 +766,16 @@ def evaluate(self, function: Function, input_ordered: bool = True) -> np.ndarray
         P0DG = fs(self.vom, "DG", 0)
         f_at_points = assemble(interpolate(function, P0DG))
 
-        if input_ordered:
-            P0DG_io = fs(self.vom.input_ordering, "DG", 0)
-            f_at_points_io = Function(P0DG_io).assign(np.nan)
-            f_at_points_io.interpolate(f_at_points)
-            return f_at_points_io.dat.data_ro
-        else:
-            return f_at_points.dat.data_ro
-
+        P0DG_io = fs(self.vom.input_ordering, "DG", 0)
+        f_at_points_io = Function(P0DG_io).assign(np.nan)
+        f_at_points_io.interpolate(f_at_points)
+        result = f_at_points_io.dat.data_ro
+        # If redundant, only rank 0 did the work. Broadcast ordered results to all ranks.
+        if self.redundant and self.mesh.comm.size > 1:
+            if self.mesh.comm.rank != 0:
+                result = np.empty((len(self.points),) + shape, dtype=utils.ScalarType)
+            self.mesh.comm.Bcast(result, root=0)
+        return result
 
 @PETSc.Log.EventDecorator()
 def make_c_evaluate(function, c_name="evaluate", ldargs=None, tolerance=None):

From 085c2254542e1c26172c7cb199678d9dd97cc50c Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 10:34:36 +0100
Subject: [PATCH 11/35] fix tests

---
 tests/firedrake/regression/test_point_eval_api.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index c85397712b..9366ecc41a 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -187,10 +187,8 @@ def test_point_evaluator_scalar(mesh_and_points):
 
     # Test standard scalar function with missing points
     eval_missing = PointEvaluator(mesh, np.append(points, [[1.5, 1.5]], axis=0), missing_points_behaviour=None)
-    f_at_points_missing = eval_missing.evaluate(f, input_ordered=False)
-    assert not np.isnan(f_at_points_missing).any()
-    f_at_points_missing_io = eval_missing.evaluate(f)
-    assert np.isnan(f_at_points_missing_io[-1])
+    f_at_points_missing = eval_missing.evaluate(f)
+    assert np.isnan(f_at_points_missing[-1])
 
 @pytest.mark.parallel([1, 3])
 def test_point_evaluator_vector_tensor_mixed(mesh_and_points):

From bfb45adb8b99254eefa4e6d7af5e772ba8610280 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 13:33:50 +0100
Subject: [PATCH 12/35] update docstrings; add check for mesh

---
 firedrake/function.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 7e6ffbfcc9..a397503bee 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -718,6 +718,10 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh. 
             "warn": warn if a point is not found in the mesh, but continue.
             None: ignore points not found in the mesh.
+        redundant : bool
+            If True, only the points on rank 0 are evaluated, and the result is broadcast to all ranks.
+            If False, each rank evaluates the points it has been given.
+            Default is True.
         """
         from firedrake.mesh import VertexOnlyMesh
         self.mesh = mesh
@@ -733,9 +737,6 @@ def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
         ----------
         function : :class:`Function`
             The :class:`Function` to evaluate.
-        input_ordered : bool
-            If ``True``, return results in the order of the input points. If any 
-            points were not found in the mesh, they will be return as np.nan.
 
         Returns
         -------
@@ -747,6 +748,8 @@ def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
         """
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got f{type(function).__name__}")
+        if function.function_space().mesh() is not self.mesh:
+            raise ValueError("Function mesh must be the same Mesh object as the PointEvaluator mesh.")
         from firedrake import interpolate, assemble, FunctionSpace, VectorFunctionSpace, TensorFunctionSpace
 
         subfunctions = function.subfunctions

From 7f1aa0916fcb943e08b1bbb97a6bce743cffdfef Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 13:34:03 +0100
Subject: [PATCH 13/35] add parallel test redundant=False

---
 .../regression/test_point_eval_api.py         | 24 ++++++++++++++++++-
 1 file changed, 23 insertions(+), 1 deletion(-)

diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index 9366ecc41a..9a5d489394 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -217,4 +217,26 @@ def test_point_evaluator_vector_tensor_mixed(mesh_and_points):
     f_tensor.interpolate(as_matrix([[x, y, x*y], [y, x, x*y]]))
     f_mixed_at_points = evaluator.evaluate(f_mixed)
     assert np.allclose(f_mixed_at_points[0], vec_expected)
-    assert np.allclose(f_mixed_at_points[1], tensor_expected)
\ No newline at end of file
+    assert np.allclose(f_mixed_at_points[1], tensor_expected)
+
+@pytest.mark.parallel(3)
+def test_point_evaluator_nonredundant(mesh_and_points):
+    mesh = mesh_and_points[0]
+    if mesh.comm.rank == 0:
+        points = [[0.1, 0.1]]
+    elif mesh.comm.rank == 1:
+        points = [[0.4, 0.4]]
+    else:
+        points = [[0.8, 0.8]]
+    V = FunctionSpace(mesh, "CG", 1)
+    f = Function(V)
+    x, y = SpatialCoordinate(mesh)
+    f.interpolate(x + y)
+    evaluator = PointEvaluator(mesh, points, redundant=False)
+    f_at_points = evaluator.evaluate(f)
+    if mesh.comm.rank == 0:
+        assert np.allclose(f_at_points, [0.2])
+    elif mesh.comm.rank == 1:
+        assert np.allclose(f_at_points, [0.8])
+    else:
+        assert np.allclose(f_at_points, [1.6])

From e9de19c7cb5adddd541a0512407ff1b7c59a8ce6 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 13:35:33 +0100
Subject: [PATCH 14/35] lint

---
 firedrake/function.py                             | 5 +++--
 tests/firedrake/regression/test_point_eval_api.py | 6 ++++--
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index a397503bee..2ab79023fd 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -702,7 +702,7 @@ def __str__(self):
 class PointEvaluator:
     r"""Convenience class for evaluating a :class:`Function` at a set of points."""
 
-    def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None, 
+    def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None,
                  missing_points_behaviour: str | None = "error", redundant: bool = True) -> None:
         r"""
         Parameters
@@ -715,7 +715,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             Tolerance to use when checking if a point is in a cell. Default is the ``tolerance`` of the :class:`Mesh`.
         missing_points_behaviour : str | None
             Behaviour when a point is not found in the mesh. Options are:
-            "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh. 
+            "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh.
             "warn": warn if a point is not found in the mesh, but continue.
             None: ignore points not found in the mesh.
         redundant : bool
@@ -780,6 +780,7 @@ def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
             self.mesh.comm.Bcast(result, root=0)
         return result
 
+
 @PETSc.Log.EventDecorator()
 def make_c_evaluate(function, c_name="evaluate", ldargs=None, tolerance=None):
     r"""Generates, compiles and loads a C function to evaluate the
diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index 9a5d489394..a6a8ad7b96 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -190,6 +190,7 @@ def test_point_evaluator_scalar(mesh_and_points):
     f_at_points_missing = eval_missing.evaluate(f)
     assert np.isnan(f_at_points_missing[-1])
 
+
 @pytest.mark.parallel([1, 3])
 def test_point_evaluator_vector_tensor_mixed(mesh_and_points):
     mesh, evaluator = mesh_and_points
@@ -205,8 +206,8 @@ def test_point_evaluator_vector_tensor_mixed(mesh_and_points):
     f_tensor = Function(V_tensor)
     f_tensor.interpolate(as_matrix([[x, y, x*y], [y, x, x*y]]))
     f_tensor_at_points = evaluator.evaluate(f_tensor)
-    tensor_expected = np.array([[[0.1, 0.1, 0.01], [0.1, 0.1, 0.01]], 
-                                [[0.2, 0.2, 0.04], [0.2, 0.2, 0.04]], 
+    tensor_expected = np.array([[[0.1, 0.1, 0.01], [0.1, 0.1, 0.01]],
+                                [[0.2, 0.2, 0.04], [0.2, 0.2, 0.04]],
                                 [[0.3, 0.3, 0.09], [0.3, 0.3, 0.09]]])
     assert np.allclose(f_tensor_at_points, tensor_expected)
 
@@ -219,6 +220,7 @@ def test_point_evaluator_vector_tensor_mixed(mesh_and_points):
     assert np.allclose(f_mixed_at_points[0], vec_expected)
     assert np.allclose(f_mixed_at_points[1], tensor_expected)
 
+
 @pytest.mark.parallel(3)
 def test_point_evaluator_nonredundant(mesh_and_points):
     mesh = mesh_and_points[0]

From 99229bb2e47be5ceed2ed5bc8c48f4f4c2ba0f00 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 14:03:54 +0100
Subject: [PATCH 15/35] tidy

---
 firedrake/function.py                             | 7 ++++---
 tests/firedrake/regression/test_point_eval_api.py | 4 ++--
 2 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 2ab79023fd..a602612d57 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -731,7 +731,8 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
 
     def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
         r"""Evaluate the given :class:`Function` at the points provided to this
-        :class:`PointEvaluator`.
+        :class:`PointEvaluator`. Points that were not found in the mesh will be
+        evaluated to np.nan.
 
         Parameters
         ----------
@@ -763,7 +764,7 @@ def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
         if len(shape) == 0:
             fs = FunctionSpace
         elif len(shape) == 1:
-            fs = VectorFunctionSpace
+            fs = partial(VectorFunctionSpace, dim=shape[0])
         else:
             fs = partial(TensorFunctionSpace, shape=shape)
         P0DG = fs(self.vom, "DG", 0)
@@ -777,7 +778,7 @@ def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
         if self.redundant and self.mesh.comm.size > 1:
             if self.mesh.comm.rank != 0:
                 result = np.empty((len(self.points),) + shape, dtype=utils.ScalarType)
-            self.mesh.comm.Bcast(result, root=0)
+            self.mesh.comm.Bcast(result)
         return result
 
 
diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index a6a8ad7b96..77f4101533 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -227,7 +227,7 @@ def test_point_evaluator_nonredundant(mesh_and_points):
     if mesh.comm.rank == 0:
         points = [[0.1, 0.1]]
     elif mesh.comm.rank == 1:
-        points = [[0.4, 0.4]]
+        points = [[0.4, 0.4], [0.5, 0.5]]
     else:
         points = [[0.8, 0.8]]
     V = FunctionSpace(mesh, "CG", 1)
@@ -239,6 +239,6 @@ def test_point_evaluator_nonredundant(mesh_and_points):
     if mesh.comm.rank == 0:
         assert np.allclose(f_at_points, [0.2])
     elif mesh.comm.rank == 1:
-        assert np.allclose(f_at_points, [0.8])
+        assert np.allclose(f_at_points, [0.8, 1.0])
     else:
         assert np.allclose(f_at_points, [1.6])

From 5ff6e34210017da8e73373b4572305c26deb5281 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 14:34:18 +0100
Subject: [PATCH 16/35] tidy

---
 firedrake/function.py                             | 8 ++++++--
 tests/firedrake/regression/test_point_eval_api.py | 3 +--
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index a602612d57..eb4418f179 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -727,7 +727,10 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
         self.redundant = redundant
-        self.vom = VertexOnlyMesh(mesh, points, missing_points_behaviour=missing_points_behaviour, redundant=redundant, tolerance=tolerance)
+        self.vom = VertexOnlyMesh(
+            mesh, points, missing_points_behaviour=missing_points_behaviour, 
+            redundant=redundant, tolerance=tolerance
+        )
 
     def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
         r"""Evaluate the given :class:`Function` at the points provided to this
@@ -745,7 +748,8 @@ def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
             A Numpy array of values at the points. If the function is scalar-valued, the Numpy array
             has shape (len(points),). If the function is vector-valued with shape (n,), the Numpy array has shape
             (len(points), n). If the function is tensor-valued with shape (n, m), the Numpy array has shape
-            (len(points), n, m). If the function is a mixed function, a list of Numpy arrays is returned, one for each subfunction.
+            (len(points), n, m). If the function is a mixed function, a list of Numpy arrays is returned, 
+            one for each subfunction.
         """
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got f{type(function).__name__}")
diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index 77f4101533..d7aba3c92f 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -181,8 +181,7 @@ def test_point_evaluator_scalar(mesh_and_points):
     f.interpolate(x + y)
 
     # Test standard scalar function evaluation at points
-    eval = PointEvaluator(mesh, points)
-    f_at_points = eval.evaluate(f)
+    f_at_points = evaluator.evaluate(f)
     assert np.allclose(f_at_points, [0.2, 0.4, 0.6])
 
     # Test standard scalar function with missing points

From f06f311e031b974fb3642f6d85e114a1095da6bb Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 14:34:59 +0100
Subject: [PATCH 17/35] lint

---
 firedrake/function.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index eb4418f179..8c479b77cd 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -728,7 +728,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
         self.points = np.asarray(points, dtype=utils.ScalarType)
         self.redundant = redundant
         self.vom = VertexOnlyMesh(
-            mesh, points, missing_points_behaviour=missing_points_behaviour, 
+            mesh, points, missing_points_behaviour=missing_points_behaviour,
             redundant=redundant, tolerance=tolerance
         )
 
@@ -748,7 +748,7 @@ def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
             A Numpy array of values at the points. If the function is scalar-valued, the Numpy array
             has shape (len(points),). If the function is vector-valued with shape (n,), the Numpy array has shape
             (len(points), n). If the function is tensor-valued with shape (n, m), the Numpy array has shape
-            (len(points), n, m). If the function is a mixed function, a list of Numpy arrays is returned, 
+            (len(points), n, m). If the function is a mixed function, a list of Numpy arrays is returned,
             one for each subfunction.
         """
         if not isinstance(function, Function):

From 8c0ddf0b40938d328a1c51e7d799ae7c7eeb9640 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 15:57:34 +0100
Subject: [PATCH 18/35] review suggestions

---
 firedrake/function.py | 26 ++++++++++++--------------
 1 file changed, 12 insertions(+), 14 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 8c479b77cd..ac63fdfa0f 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -12,6 +12,7 @@
 from numbers import Number
 from pathlib import Path
 from functools import partial
+from typing import Tuple
 
 from pyop2 import op2, mpi
 from pyop2.exceptions import DataTypeError, DataValueError
@@ -712,7 +713,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
         points : np.ndarray | list
             Array of points to evaluate the function at.
         tolerance : float | None
-            Tolerance to use when checking if a point is in a cell. Default is the ``tolerance`` of the :class:`Mesh`.
+            Tolerance to use when checking if a point is in a cell. Default is the `tolerance` of the :class:`MeshGeometry`.
         missing_points_behaviour : str | None
             Behaviour when a point is not found in the mesh. Options are:
             "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh.
@@ -720,7 +721,8 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             None: ignore points not found in the mesh.
         redundant : bool
             If True, only the points on rank 0 are evaluated, and the result is broadcast to all ranks.
-            If False, each rank evaluates the points it has been given.
+            If False, each rank evaluates the points it has been given. False is useful if you are inputting
+            external data that is already distributed across ranks.
             Default is True.
         """
         from firedrake.mesh import VertexOnlyMesh
@@ -732,23 +734,22 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             redundant=redundant, tolerance=tolerance
         )
 
-    def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
-        r"""Evaluate the given :class:`Function` at the points provided to this
-        :class:`PointEvaluator`. Points that were not found in the mesh will be
-        evaluated to np.nan.
+    def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray]:
+        r"""Evaluate the :class:`Function`.
+        Points that were not found in the mesh will be evaluated to np.nan.
 
         Parameters
         ----------
-        function : :class:`Function`
-            The :class:`Function` to evaluate.
+        function :
+            The `Function` to evaluate.
 
         Returns
         -------
-        np.ndarray | list[np.ndarray]
+        np.ndarray | Tuple[np.ndarray]
             A Numpy array of values at the points. If the function is scalar-valued, the Numpy array
             has shape (len(points),). If the function is vector-valued with shape (n,), the Numpy array has shape
             (len(points), n). If the function is tensor-valued with shape (n, m), the Numpy array has shape
-            (len(points), n, m). If the function is a mixed function, a list of Numpy arrays is returned,
+            (len(points), n, m). If the function is a mixed function, a tuple of Numpy arrays is returned,
             one for each subfunction.
         """
         if not isinstance(function, Function):
@@ -759,10 +760,7 @@ def evaluate(self, function: Function) -> np.ndarray | list[np.ndarray]:
 
         subfunctions = function.subfunctions
         if len(subfunctions) > 1:
-            result = []
-            for subfunction in subfunctions:
-                result.append(self.evaluate(subfunction))
-            return result
+            return tuple(self.evaluate(subfunction) for subfunction in subfunctions)
 
         shape = function.ufl_function_space().value_shape
         if len(shape) == 0:

From 5194ea582b04df5df970b9a6a061b0fcb677d447 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 16:06:34 +0100
Subject: [PATCH 19/35] move imports

---
 firedrake/function.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index ac63fdfa0f..ab37b22a6a 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -25,7 +25,8 @@
 from firedrake import utils
 from firedrake.adjoint_utils import FunctionMixin
 from firedrake.petsc import PETSc
-from firedrake.mesh import MeshGeometry
+from firedrake.mesh import MeshGeometry, VertexOnlyMesh
+from firedrake.functionspace import FunctionSpace, VectorFunctionSpace, TensorFunctionSpace
 
 __all__ = ['Function', 'PointNotInDomainError', 'CoordinatelessFunction', 'PointEvaluator']
 
@@ -725,7 +726,6 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             external data that is already distributed across ranks.
             Default is True.
         """
-        from firedrake.mesh import VertexOnlyMesh
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
         self.redundant = redundant
@@ -752,12 +752,11 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray]:
             (len(points), n, m). If the function is a mixed function, a tuple of Numpy arrays is returned,
             one for each subfunction.
         """
+        from firedrake import assemble, interpolate
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got f{type(function).__name__}")
         if function.function_space().mesh() is not self.mesh:
             raise ValueError("Function mesh must be the same Mesh object as the PointEvaluator mesh.")
-        from firedrake import interpolate, assemble, FunctionSpace, VectorFunctionSpace, TensorFunctionSpace
-
         subfunctions = function.subfunctions
         if len(subfunctions) > 1:
             return tuple(self.evaluate(subfunction) for subfunction in subfunctions)

From 317dd1eb3a342f391938b9794c25285d415b6e7c Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 16:07:33 +0100
Subject: [PATCH 20/35] Change missing points ignore

fix test
---
 firedrake/function.py                             | 6 +++---
 tests/firedrake/regression/test_point_eval_api.py | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index ab37b22a6a..0c58101d9d 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -705,7 +705,7 @@ class PointEvaluator:
     r"""Convenience class for evaluating a :class:`Function` at a set of points."""
 
     def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None,
-                 missing_points_behaviour: str | None = "error", redundant: bool = True) -> None:
+                 missing_points_behaviour: str = "error", redundant: bool = True) -> None:
         r"""
         Parameters
         ----------
@@ -715,11 +715,11 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             Array of points to evaluate the function at.
         tolerance : float | None
             Tolerance to use when checking if a point is in a cell. Default is the `tolerance` of the :class:`MeshGeometry`.
-        missing_points_behaviour : str | None
+        missing_points_behaviour : str
             Behaviour when a point is not found in the mesh. Options are:
             "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh.
             "warn": warn if a point is not found in the mesh, but continue.
-            None: ignore points not found in the mesh.
+            "ignore": ignore points not found in the mesh.
         redundant : bool
             If True, only the points on rank 0 are evaluated, and the result is broadcast to all ranks.
             If False, each rank evaluates the points it has been given. False is useful if you are inputting
diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index d7aba3c92f..69c46c4a32 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -185,7 +185,7 @@ def test_point_evaluator_scalar(mesh_and_points):
     assert np.allclose(f_at_points, [0.2, 0.4, 0.6])
 
     # Test standard scalar function with missing points
-    eval_missing = PointEvaluator(mesh, np.append(points, [[1.5, 1.5]], axis=0), missing_points_behaviour=None)
+    eval_missing = PointEvaluator(mesh, np.append(points, [[1.5, 1.5]], axis=0), missing_points_behaviour="ignore")
     f_at_points_missing = eval_missing.evaluate(f)
     assert np.isnan(f_at_points_missing[-1])
 

From ac6d0b10fc18c945496c3ff1d8cc56ef2f389540 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 16:12:59 +0100
Subject: [PATCH 21/35] fix tuple type hint

---
 firedrake/function.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 0c58101d9d..4fcba20be0 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -734,7 +734,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             redundant=redundant, tolerance=tolerance
         )
 
-    def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray]:
+    def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         r"""Evaluate the :class:`Function`.
         Points that were not found in the mesh will be evaluated to np.nan.
 
@@ -745,7 +745,7 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray]:
 
         Returns
         -------
-        np.ndarray | Tuple[np.ndarray]
+        np.ndarray | Tuple[np.ndarray, ...]
             A Numpy array of values at the points. If the function is scalar-valued, the Numpy array
             has shape (len(points),). If the function is vector-valued with shape (n,), the Numpy array has shape
             (len(points), n). If the function is tensor-valued with shape (n, m), the Numpy array has shape

From 393ccb37294622a3398a32f17c9c08a50e108e55 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Fri, 22 Aug 2025 23:07:30 +0100
Subject: [PATCH 22/35] revert missing points behaviour change

---
 firedrake/function.py                             | 6 +++---
 firedrake/mesh.py                                 | 2 +-
 tests/firedrake/regression/test_point_eval_api.py | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 4fcba20be0..2f476b84c8 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -705,7 +705,7 @@ class PointEvaluator:
     r"""Convenience class for evaluating a :class:`Function` at a set of points."""
 
     def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None,
-                 missing_points_behaviour: str = "error", redundant: bool = True) -> None:
+                 missing_points_behaviour: str | None = "error", redundant: bool = True) -> None:
         r"""
         Parameters
         ----------
@@ -715,11 +715,11 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             Array of points to evaluate the function at.
         tolerance : float | None
             Tolerance to use when checking if a point is in a cell. Default is the `tolerance` of the :class:`MeshGeometry`.
-        missing_points_behaviour : str
+        missing_points_behaviour : str | None
             Behaviour when a point is not found in the mesh. Options are:
             "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh.
             "warn": warn if a point is not found in the mesh, but continue.
-            "ignore": ignore points not found in the mesh.
+            None: ignore points not found in the mesh.
         redundant : bool
             If True, only the points on rank 0 are evaluated, and the result is broadcast to all ranks.
             If False, each rank evaluates the points it has been given. False is useful if you are inputting
diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 7358702adc..d9220879bc 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3319,7 +3319,7 @@ def ExtrudedMesh(mesh, layers, layer_height=None, extrusion_type='uniform', peri
 
 
 class MissingPointsBehaviour(enum.Enum):
-    IGNORE = "ignore"
+    IGNORE = None
     ERROR = "error"
     WARN = "warn"
 
diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index 69c46c4a32..d7aba3c92f 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -185,7 +185,7 @@ def test_point_evaluator_scalar(mesh_and_points):
     assert np.allclose(f_at_points, [0.2, 0.4, 0.6])
 
     # Test standard scalar function with missing points
-    eval_missing = PointEvaluator(mesh, np.append(points, [[1.5, 1.5]], axis=0), missing_points_behaviour="ignore")
+    eval_missing = PointEvaluator(mesh, np.append(points, [[1.5, 1.5]], axis=0), missing_points_behaviour=None)
     f_at_points_missing = eval_missing.evaluate(f)
     assert np.isnan(f_at_points_missing[-1])
 

From 09b733d6a433e73c45c1a00557ac4489b8ff4b1b Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Wed, 27 Aug 2025 17:39:11 +0100
Subject: [PATCH 23/35] clarify docstring

fix comment
---
 firedrake/function.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 2f476b84c8..78d185ea72 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -721,7 +721,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             "warn": warn if a point is not found in the mesh, but continue.
             None: ignore points not found in the mesh.
         redundant : bool
-            If True, only the points on rank 0 are evaluated, and the result is broadcast to all ranks.
+            If True, only the points given to the constructor on rank 0 are evaluated, and the result is broadcast to all ranks.
             If False, each rank evaluates the points it has been given. False is useful if you are inputting
             external data that is already distributed across ranks.
             Default is True.
@@ -755,8 +755,11 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         from firedrake import assemble, interpolate
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got f{type(function).__name__}")
-        if function.function_space().mesh() is not self.mesh:
+        function_mesh = function.function_space().mesh()
+        if function_mesh is not self.mesh:
             raise ValueError("Function mesh must be the same Mesh object as the PointEvaluator mesh.")
+        # if function_mesh.coordinates.dat.dat_version != self.mesh.coordinates.dat.dat_version:
+        #     pass
         subfunctions = function.subfunctions
         if len(subfunctions) > 1:
             return tuple(self.evaluate(subfunction) for subfunction in subfunctions)
@@ -775,7 +778,7 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         f_at_points_io = Function(P0DG_io).assign(np.nan)
         f_at_points_io.interpolate(f_at_points)
         result = f_at_points_io.dat.data_ro
-        # If redundant, only rank 0 did the work. Broadcast ordered results to all ranks.
+        # If redundant, all points are now on rank 0, so we broadcast the result
         if self.redundant and self.mesh.comm.size > 1:
             if self.mesh.comm.rank != 0:
                 result = np.empty((len(self.points),) + shape, dtype=utils.ScalarType)

From 6358bf473003251e009a14209537c6db5fff01c5 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 28 Aug 2025 14:02:58 +0100
Subject: [PATCH 24/35] add to docs

fixes
---
 docs/source/point-evaluation.rst | 35 ++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/docs/source/point-evaluation.rst b/docs/source/point-evaluation.rst
index 731a0081f5..6ba398ccf0 100644
--- a/docs/source/point-evaluation.rst
+++ b/docs/source/point-evaluation.rst
@@ -214,6 +214,7 @@ be a new vertex: this is true for both ``redundant = True`` and
 and switch from ``redundant = True`` to ``redundant = False`` we will get point
 duplication.
 
+.. _missing_points:
 
 Points outside the domain
 ~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -242,6 +243,40 @@ warning or switched off entirely:
    :end-before: [test_vom_manual_points_outside_domain 6]
 
 
+``PointEvaluator`` convenience object
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The :py:class:`~.PointEvaluator` class performs point evaluation using vertex-only meshes,
+as described above. First, create a :py:class:`~.PointEvaluator` object by passing the 
+parent mesh and the points to evaluate at:
+
+.. code-block:: python3
+
+   point_evaluator = PointEvaluator(mesh, points)
+
+Then to evaluate a :py:class:`~.Function` defined on the parent mesh at the given points,
+we use :meth:`~.PointEvaluator.evaluate`:
+
+.. code-block:: python3
+
+   f_at_points = point_evaluator.evaluate(f)
+
+Under the hood, this creates the appropriate P0DG function space on the vertex-only mesh
+and performs the interpolation. The points are then reordered to match the input ordering,
+as described in :ref:`the section on the input ordering property <input_ordering>`.
+
+If ``redundant=True`` (the default) was used when creating the :py:class:`~.PointEvaluator`,
+then only the points given to the constructor on rank 0 will be evaluated. The result is then
+broadcast to all ranks. Use this option if the same points are given on all ranks.
+If ``redundant=False`` was used when creating the :py:class:`~.PointEvaluator`, then
+each rank will evaluate the points it was given. Use this option if different points are given
+on different ranks, for example when using external point data.
+
+The parameters ``missing_points_behaviour`` and ``tolerance`` (discussed :ref:`here <missing_points>` 
+and :ref:`here <tolerance>` respectively) can be set when creating the :py:class:`~.PointEvaluator` 
+and will be passed to the :func:`~.VertexOnlyMesh` it creates internally.
+
+
 Expressions with point evaluations
 ----------------------------------
 

From f1a212fc7ebf53dc5487a6a92900ccc0307791a3 Mon Sep 17 00:00:00 2001
From: Leo <leocollins511@gmail.com>
Date: Thu, 28 Aug 2025 14:13:15 +0100
Subject: [PATCH 25/35] tidy up docstrings

---
 firedrake/function.py | 47 +++++++++++++++++++++++++------------------
 1 file changed, 27 insertions(+), 20 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 78d185ea72..c5914de419 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -705,61 +705,68 @@ class PointEvaluator:
     r"""Convenience class for evaluating a :class:`Function` at a set of points."""
 
     def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: float | None = None,
-                 missing_points_behaviour: str | None = "error", redundant: bool = True) -> None:
+                 missing_points_behaviour: str = "error", redundant: bool = True) -> None:
         r"""
         Parameters
         ----------
-        mesh : :class:`MeshGeometry`
-            The :class:`Mesh` on which the :class:`Function` is defined.
-        points : np.ndarray | list
-            Array of points to evaluate the function at.
+        mesh : MeshGeometry
+            The mesh on which to embed the points.
+        points : numpy.ndarray | list
+            Array or list of points to evaluate at.
         tolerance : float | None
-            Tolerance to use when checking if a point is in a cell. Default is the `tolerance` of the :class:`MeshGeometry`.
-        missing_points_behaviour : str | None
+            Tolerance to use when checking if a point is in a cell.
+            If ``None`` (the default), the ``tolerance`` of the ``mesh`` is used.
+        missing_points_behaviour : str
             Behaviour when a point is not found in the mesh. Options are:
-            "error": raise a :class:`VertexOnlyMeshMissingPointsError` if a point is not found in the mesh.
+            "error": raise a :class:`~.VertexOnlyMeshMissingPointsError` if a point is not found in the mesh.
             "warn": warn if a point is not found in the mesh, but continue.
-            None: ignore points not found in the mesh.
+            "ignore": ignore points not found in the mesh.
         redundant : bool
             If True, only the points given to the constructor on rank 0 are evaluated, and the result is broadcast to all ranks.
             If False, each rank evaluates the points it has been given. False is useful if you are inputting
-            external data that is already distributed across ranks.
-            Default is True.
+            external data that is already distributed across ranks. Default is True.
         """
         self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
         self.redundant = redundant
+        self.missing_points_behaviour = missing_points_behaviour
+        self.tolerance = tolerance
         self.vom = VertexOnlyMesh(
             mesh, points, missing_points_behaviour=missing_points_behaviour,
             redundant=redundant, tolerance=tolerance
         )
 
     def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
-        r"""Evaluate the :class:`Function`.
+        r"""Evaluate the given :class:`Function`.
         Points that were not found in the mesh will be evaluated to np.nan.
 
         Parameters
         ----------
         function :
-            The `Function` to evaluate.
+            The :class:`Function` to evaluate.
 
         Returns
         -------
-        np.ndarray | Tuple[np.ndarray, ...]
+        numpy.ndarray | Tuple[numpy.ndarray, ...]
             A Numpy array of values at the points. If the function is scalar-valued, the Numpy array
-            has shape (len(points),). If the function is vector-valued with shape (n,), the Numpy array has shape
-            (len(points), n). If the function is tensor-valued with shape (n, m), the Numpy array has shape
-            (len(points), n, m). If the function is a mixed function, a tuple of Numpy arrays is returned,
+            has shape ``(len(points),)``. If the function is vector-valued with shape ``(n,)``, the Numpy array has shape
+            ``(len(points), n)``. If the function is tensor-valued with shape ``(n, m)``, the Numpy array has shape
+            ``(len(points), n, m)``. If the function is a mixed function, a tuple of Numpy arrays is returned,
             one for each subfunction.
         """
         from firedrake import assemble, interpolate
         if not isinstance(function, Function):
-            raise TypeError(f"Expected a Function, got f{type(function).__name__}")
+            raise TypeError(f"Expected a Function, got {type(function).__name__}")
         function_mesh = function.function_space().mesh()
         if function_mesh is not self.mesh:
             raise ValueError("Function mesh must be the same Mesh object as the PointEvaluator mesh.")
-        # if function_mesh.coordinates.dat.dat_version != self.mesh.coordinates.dat.dat_version:
-        #     pass
+        if function_mesh.coordinates.dat.dat_version != self.mesh._saved_coordinate_dat_version:
+            # The mesh coordinates have changed since the PointEvaluator was created
+            self.mesh = function_mesh
+            self.vom = VertexOnlyMesh(
+                self.mesh, self.points, missing_points_behaviour=self.missing_points_behaviour,
+                redundant=self.redundant, tolerance=self.tolerance
+            )
         subfunctions = function.subfunctions
         if len(subfunctions) > 1:
             return tuple(self.evaluate(subfunction) for subfunction in subfunctions)

From 46e1185d1f2a48c16dc9a983a5233162afcce94b Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 28 Aug 2025 15:47:27 +0100
Subject: [PATCH 26/35] fix missing points behaviour

---
 tests/firedrake/regression/test_point_eval_api.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index d7aba3c92f..69c46c4a32 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -185,7 +185,7 @@ def test_point_evaluator_scalar(mesh_and_points):
     assert np.allclose(f_at_points, [0.2, 0.4, 0.6])
 
     # Test standard scalar function with missing points
-    eval_missing = PointEvaluator(mesh, np.append(points, [[1.5, 1.5]], axis=0), missing_points_behaviour=None)
+    eval_missing = PointEvaluator(mesh, np.append(points, [[1.5, 1.5]], axis=0), missing_points_behaviour="ignore")
     f_at_points_missing = eval_missing.evaluate(f)
     assert np.isnan(f_at_points_missing[-1])
 

From 367a769b8d8b6e46aae1128758a529f132b26431 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Mon, 1 Sep 2025 16:30:51 +0100
Subject: [PATCH 27/35] add moving mesh test

---
 .../regression/test_point_eval_api.py         | 20 +++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index 69c46c4a32..1f1451f587 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -3,6 +3,7 @@
 import pytest
 
 from firedrake import *
+from firedrake.mesh import VertexOnlyMeshMissingPointsError
 
 cwd = abspath(dirname(__file__))
 
@@ -241,3 +242,22 @@ def test_point_evaluator_nonredundant(mesh_and_points):
         assert np.allclose(f_at_points, [0.8, 1.0])
     else:
         assert np.allclose(f_at_points, [1.6])
+
+
+def test_point_evaluator_moving_mesh(mesh_and_points):
+    mesh, evaluator = mesh_and_points
+    V = FunctionSpace(mesh, "CG", 1)
+    f = Function(V)
+    x, y = SpatialCoordinate(mesh)
+    f.interpolate(x + y)
+
+    mesh.coordinates.dat.data[:, 0] += 1.0
+
+    with pytest.raises(VertexOnlyMeshMissingPointsError):
+        # The VOM is reimmersed, but the points
+        # are now outside of the mesh.
+        f_at_points = evaluator.evaluate(f)
+
+    mesh.coordinates.dat.data[:, 0] -= 1.0
+    f_at_points = evaluator.evaluate(f)
+    assert np.allclose(f_at_points, [0.2, 0.4, 0.6])

From 6ef7f7982603ed1c215fdf461129766a50b28961 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Mon, 1 Sep 2025 16:42:06 +0100
Subject: [PATCH 28/35] moving mesh docs

---
 docs/source/point-evaluation.rst | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/docs/source/point-evaluation.rst b/docs/source/point-evaluation.rst
index 6ba398ccf0..0ceca58504 100644
--- a/docs/source/point-evaluation.rst
+++ b/docs/source/point-evaluation.rst
@@ -254,7 +254,8 @@ parent mesh and the points to evaluate at:
 
    point_evaluator = PointEvaluator(mesh, points)
 
-Then to evaluate a :py:class:`~.Function` defined on the parent mesh at the given points,
+Internally, this creates a vertex-only mesh at the given points, immersed in the given mesh.
+To evaluate a :py:class:`~.Function` defined on the parent mesh at the given points,
 we use :meth:`~.PointEvaluator.evaluate`:
 
 .. code-block:: python3
@@ -276,6 +277,9 @@ The parameters ``missing_points_behaviour`` and ``tolerance`` (discussed :ref:`h
 and :ref:`here <tolerance>` respectively) can be set when creating the :py:class:`~.PointEvaluator` 
 and will be passed to the :func:`~.VertexOnlyMesh` it creates internally.
 
+If the :ref:`mesh coordinates are changed <changing_coordinate_fs>`, then the vertex-only mesh
+will be re-immersed in the new mesh the next time :meth:`~.PointEvaluator.evaluate` is called.
+
 
 Expressions with point evaluations
 ----------------------------------

From 95359aa907f87eae0656427da35f0efcac206259 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Mon, 1 Sep 2025 18:25:39 +0100
Subject: [PATCH 29/35] fix

---
 firedrake/mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index d9220879bc..7358702adc 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -3319,7 +3319,7 @@ def ExtrudedMesh(mesh, layers, layer_height=None, extrusion_type='uniform', peri
 
 
 class MissingPointsBehaviour(enum.Enum):
-    IGNORE = None
+    IGNORE = "ignore"
     ERROR = "error"
     WARN = "warn"
 

From be2d89845658c572853cc023fcd00f1559312f55 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 2 Sep 2025 14:08:50 +0100
Subject: [PATCH 30/35] tidy

tidy
---
 firedrake/function.py | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index c5914de419..44e52632a9 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -757,16 +757,20 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         from firedrake import assemble, interpolate
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got {type(function).__name__}")
+
         function_mesh = function.function_space().mesh()
         if function_mesh is not self.mesh:
             raise ValueError("Function mesh must be the same Mesh object as the PointEvaluator mesh.")
-        if function_mesh.coordinates.dat.dat_version != self.mesh._saved_coordinate_dat_version:
-            # The mesh coordinates have changed since the PointEvaluator was created
+        if coord_changed := function_mesh.coordinates.dat.dat_version != self.mesh._saved_coordinate_dat_version:
             self.mesh = function_mesh
+        if tol_changed := self.mesh.tolerance != self.tolerance:
+            self.tolerance = self.mesh.tolerance
+        if coord_changed or tol_changed:
             self.vom = VertexOnlyMesh(
                 self.mesh, self.points, missing_points_behaviour=self.missing_points_behaviour,
                 redundant=self.redundant, tolerance=self.tolerance
             )
+
         subfunctions = function.subfunctions
         if len(subfunctions) > 1:
             return tuple(self.evaluate(subfunction) for subfunction in subfunctions)
@@ -779,12 +783,13 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         else:
             fs = partial(TensorFunctionSpace, shape=shape)
         P0DG = fs(self.vom, "DG", 0)
-        f_at_points = assemble(interpolate(function, P0DG))
-
         P0DG_io = fs(self.vom.input_ordering, "DG", 0)
+
+        f_at_points = assemble(interpolate(function, P0DG))
         f_at_points_io = Function(P0DG_io).assign(np.nan)
         f_at_points_io.interpolate(f_at_points)
         result = f_at_points_io.dat.data_ro
+
         # If redundant, all points are now on rank 0, so we broadcast the result
         if self.redundant and self.mesh.comm.size > 1:
             if self.mesh.comm.rank != 0:

From b16232a3ed029f62e992e8213d1dfbc321ecbabd Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 2 Sep 2025 14:14:29 +0100
Subject: [PATCH 31/35] add changing tolerance test

lint
---
 .../regression/test_point_eval_api.py         | 21 +++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index 1f1451f587..2de6b80030 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -261,3 +261,24 @@ def test_point_evaluator_moving_mesh(mesh_and_points):
     mesh.coordinates.dat.data[:, 0] -= 1.0
     f_at_points = evaluator.evaluate(f)
     assert np.allclose(f_at_points, [0.2, 0.4, 0.6])
+
+
+def test_point_evaluator_tolerance():
+    mesh = UnitSquareMesh(1, 1)
+    old_tol = mesh.tolerance
+    f = Function(VectorFunctionSpace(mesh, "CG", 1)).interpolate(SpatialCoordinate(mesh))
+    ev = PointEvaluator(mesh, [[0.2, 0.4]])
+    assert np.allclose([0.2, 0.4], ev.evaluate(f))
+    # tolerance of mesh is not changed
+    assert mesh.tolerance == old_tol
+    # Outside mesh, but within tolerance
+    ev = PointEvaluator(mesh, [[-0.1, 0.4]], tolerance=0.2)
+    assert np.allclose([-0.1, 0.4], ev.evaluate(f))
+    # tolerance of mesh is changed
+    assert mesh.tolerance == 0.2
+    # works if mesh tolerance is changed
+    mesh.tolerance = 1e-11
+    with pytest.raises(VertexOnlyMeshMissingPointsError):
+        ev.evaluate(f)
+    ev = PointEvaluator(mesh, [[-1e-12, 0.4]])
+    assert np.allclose([-1e-12, 0.4], ev.evaluate(f))

From bbe8834b702ce4a943d669e3944a512e6e2ca0ff Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Wed, 3 Sep 2025 18:01:55 +0100
Subject: [PATCH 32/35] add warning

---
 docs/source/point-evaluation.rst | 2 +-
 firedrake/function.py            | 8 ++++++++
 2 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/docs/source/point-evaluation.rst b/docs/source/point-evaluation.rst
index 0ceca58504..0d2c12ae8b 100644
--- a/docs/source/point-evaluation.rst
+++ b/docs/source/point-evaluation.rst
@@ -278,7 +278,7 @@ and :ref:`here <tolerance>` respectively) can be set when creating the :py:class
 and will be passed to the :func:`~.VertexOnlyMesh` it creates internally.
 
 If the :ref:`mesh coordinates are changed <changing_coordinate_fs>`, then the vertex-only mesh
-will be re-immersed in the new mesh the next time :meth:`~.PointEvaluator.evaluate` is called.
+will be reconstructed on the new mesh the next time :meth:`~.PointEvaluator.evaluate` is called.
 
 
 Expressions with point evaluations
diff --git a/firedrake/function.py b/firedrake/function.py
index 44e52632a9..1f72ca49c5 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -753,6 +753,13 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
             ``(len(points), n)``. If the function is tensor-valued with shape ``(n, m)``, the Numpy array has shape
             ``(len(points), n, m)``. If the function is a mixed function, a tuple of Numpy arrays is returned,
             one for each subfunction.
+
+
+        .. warning::
+
+            This method returns a numpy array and hence isn't taped for use with firedrake-adjoint.
+            If you want to use point evaluation with the adjoint, create a :class:`~.VertexOnlyMesh`
+            as described in the manual.
         """
         from firedrake import assemble, interpolate
         if not isinstance(function, Function):
@@ -762,6 +769,7 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         if function_mesh is not self.mesh:
             raise ValueError("Function mesh must be the same Mesh object as the PointEvaluator mesh.")
         if coord_changed := function_mesh.coordinates.dat.dat_version != self.mesh._saved_coordinate_dat_version:
+            # TODO: This is here until https://github.com/firedrakeproject/firedrake/issues/4540 is solved
             self.mesh = function_mesh
         if tol_changed := self.mesh.tolerance != self.tolerance:
             self.tolerance = self.mesh.tolerance

From 1af9c8f3594d80b93f45427fc8548d9ba2033bd8 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 4 Sep 2025 14:18:08 +0100
Subject: [PATCH 33/35] add explanation

---
 docs/source/point-evaluation.rst | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/docs/source/point-evaluation.rst b/docs/source/point-evaluation.rst
index 0d2c12ae8b..03e612cf67 100644
--- a/docs/source/point-evaluation.rst
+++ b/docs/source/point-evaluation.rst
@@ -242,12 +242,16 @@ warning or switched off entirely:
    :start-after: [test_vom_manual_points_outside_domain 5]
    :end-before: [test_vom_manual_points_outside_domain 6]
 
+.. _point_evaluator:
 
 ``PointEvaluator`` convenience object
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The :py:class:`~.PointEvaluator` class performs point evaluation using vertex-only meshes,
-as described above. First, create a :py:class:`~.PointEvaluator` object by passing the 
+as described above. This is a convenience object for users who want to evaluate
+functions at points without writing the vertex-only mesh boilerplate code each time.
+
+First, create a :py:class:`~.PointEvaluator` object by passing the 
 parent mesh and the points to evaluate at:
 
 .. code-block:: python3
@@ -264,7 +268,9 @@ we use :meth:`~.PointEvaluator.evaluate`:
 
 Under the hood, this creates the appropriate P0DG function space on the vertex-only mesh
 and performs the interpolation. The points are then reordered to match the input ordering,
-as described in :ref:`the section on the input ordering property <input_ordering>`.
+as described in :ref:`the section on the input ordering property <input_ordering>`. The result
+is a Numpy array containing the values of ``f`` at the given points, in the order the points were
+provided to the :py:class:`~.PointEvaluator` constructor.
 
 If ``redundant=True`` (the default) was used when creating the :py:class:`~.PointEvaluator`,
 then only the points given to the constructor on rank 0 will be evaluated. The result is then
@@ -277,7 +283,8 @@ The parameters ``missing_points_behaviour`` and ``tolerance`` (discussed :ref:`h
 and :ref:`here <tolerance>` respectively) can be set when creating the :py:class:`~.PointEvaluator` 
 and will be passed to the :func:`~.VertexOnlyMesh` it creates internally.
 
-If the :ref:`mesh coordinates are changed <changing_coordinate_fs>`, then the vertex-only mesh
+If the :ref:`coordinates <changing_coordinate_fs>` or the :ref:`tolerance <tolerance>` of the parent mesh
+are changed after creating the :py:class:`~.PointEvaluator`, then the vertex-only mesh
 will be reconstructed on the new mesh the next time :meth:`~.PointEvaluator.evaluate` is called.
 
 

From 3d3ac5a17c619ace9ae64e9a362c2ea7a6429d2a Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 4 Sep 2025 15:21:30 +0100
Subject: [PATCH 34/35] sanitise input

---
 firedrake/function.py | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 1f72ca49c5..e0fe2d1b5c 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -726,13 +726,21 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
             If False, each rank evaluates the points it has been given. False is useful if you are inputting
             external data that is already distributed across ranks. Default is True.
         """
-        self.mesh = mesh
         self.points = np.asarray(points, dtype=utils.ScalarType)
+        if not self.points.shape:
+            self.points = self.points.reshape(-1)
+        gdim = mesh.geometric_dimension()
+        if self.points.shape[-1] != gdim and (len(self.points.shape) != 1 or gdim != 1):
+            raise ValueError(f"Point dimension ({self.points.shape[-1]}) does not match geometric dimension ({gdim}).")
+        self.points = self.points.reshape(-1, gdim)
+
+        self.mesh = mesh
+        
         self.redundant = redundant
         self.missing_points_behaviour = missing_points_behaviour
         self.tolerance = tolerance
         self.vom = VertexOnlyMesh(
-            mesh, points, missing_points_behaviour=missing_points_behaviour,
+            mesh, self.points, missing_points_behaviour=missing_points_behaviour,
             redundant=redundant, tolerance=tolerance
         )
 
@@ -764,6 +772,8 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         from firedrake import assemble, interpolate
         if not isinstance(function, Function):
             raise TypeError(f"Expected a Function, got {type(function).__name__}")
+        if function.function_space().ufl_element().family() == "Real":
+            return function.dat.data_ro
 
         function_mesh = function.function_space().mesh()
         if function_mesh is not self.mesh:

From 9e3ea8bc4f5ea4b16323e481c7ef03be6e498919 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Thu, 4 Sep 2025 15:22:26 +0100
Subject: [PATCH 35/35] add tests

---
 firedrake/function.py                         |  2 +-
 .../regression/test_point_eval_api.py         | 59 +++++++++++++++++++
 2 files changed, 60 insertions(+), 1 deletion(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index e0fe2d1b5c..a508dddc8b 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -735,7 +735,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
         self.points = self.points.reshape(-1, gdim)
 
         self.mesh = mesh
-        
+
         self.redundant = redundant
         self.missing_points_behaviour = missing_points_behaviour
         self.tolerance = tolerance
diff --git a/tests/firedrake/regression/test_point_eval_api.py b/tests/firedrake/regression/test_point_eval_api.py
index 2de6b80030..22d177cb04 100644
--- a/tests/firedrake/regression/test_point_eval_api.py
+++ b/tests/firedrake/regression/test_point_eval_api.py
@@ -282,3 +282,62 @@ def test_point_evaluator_tolerance():
         ev.evaluate(f)
     ev = PointEvaluator(mesh, [[-1e-12, 0.4]])
     assert np.allclose([-1e-12, 0.4], ev.evaluate(f))
+
+
+def test_point_evaluator_inputs_1d():
+    mesh = UnitIntervalMesh(1)
+    f = mesh.coordinates
+
+    # one point
+    for input in [0.2, (0.2,), [0.2], np.array([0.2])]:
+        e = PointEvaluator(mesh, input)
+        assert np.allclose(0.2, e.evaluate(f))
+
+    # multiple points as tuples/list
+    for input in [
+        (0.2, 0.3), ((0.2,), (0.3,)), ([0.2], [0.3]),
+        (np.array(0.2), np.array(0.3)), (np.array([0.2]), np.array([0.3]))
+    ]:
+        e2 = PointEvaluator(mesh, input)
+        assert np.allclose([[0.2, 0.3]], e2.evaluate(f))
+        e3 = PointEvaluator(mesh, list(input))
+        assert np.allclose([[0.2, 0.3]], e3.evaluate(f))
+
+    # multiple points as numpy array
+    for input in [np.array([0.2, 0.3]), np.array([[0.2], [0.3]])]:
+        e = PointEvaluator(mesh, input)
+        assert np.allclose([[0.2, 0.3]], e.evaluate(f))
+
+    # test incorrect inputs
+    for input in [[[0.2, 0.3]], ([0.2, 0.3], [0.4, 0.5]), np.array([[0.2, 0.3]])]:
+        with pytest.raises(ValueError):
+            PointEvaluator(mesh, input)
+
+
+def test_point_evaluator_inputs_2d():
+    mesh = UnitSquareMesh(1, 1)
+    f = mesh.coordinates
+
+    # one point
+    for input in [(0.2, 0.4), [0.2, 0.4], [[0.2, 0.4]], np.array([0.2, 0.4])]:
+        e = PointEvaluator(mesh, input)
+        assert np.allclose([0.2, 0.4], e.evaluate(f))
+
+    # multiple points as tuple
+    for input in [
+        ((0.2, 0.4), (0.3, 0.5)), ([0.2, 0.4], [0.3, 0.5]),
+        (np.array([0.2, 0.4]), np.array([0.3, 0.5]))
+    ]:
+        e1 = PointEvaluator(mesh, input)
+        assert np.allclose([[0.2, 0.4], [0.3, 0.5]], e1.evaluate(f))
+        e2 = PointEvaluator(mesh, list(input))
+        assert np.allclose([[0.2, 0.4], [0.3, 0.5]], e2.evaluate(f))
+
+    # multiple points as numpy array
+    e = PointEvaluator(mesh, np.array([[0.2, 0.4], [0.3, 0.5]]))
+    assert np.allclose([[0.2, 0.4], [0.3, 0.5]], e.evaluate(f))
+
+    # test incorrect inputs
+    for input in [0.2, [0.2]]:
+        with pytest.raises(ValueError):
+            PointEvaluator(mesh, input)