From 5f081d8145837d3a198fbd4592dba87e21056f9d Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Mon, 15 Jun 2026 12:15:39 +0100
Subject: [PATCH 01/15] cache bounding box coords

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

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index b6aa2e3f15..a78e8296bb 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2496,6 +2496,8 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         Bezier curves and are completely contained in the convex hull of the mesh nodes.
         Hence the bounding box will contain the entire element.
         """
+        if self._bounding_box_coords:
+            return self._bounding_box_coords
         from firedrake import function, functionspace
         from firedrake.parloops import par_loop, READ, MIN, MAX
 
@@ -2530,7 +2532,8 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         cell_node_list = mesh.coordinates.function_space().cell_node_list
         if not mesh.extruded:
             all_coords = coords.dat.data_ro_with_halos[cell_node_list]
-            return np.min(all_coords, axis=1), np.max(all_coords, axis=1)
+            self._bounding_box_coords = np.min(all_coords, axis=1), np.max(all_coords, axis=1)
+            return self._bounding_box_coords
 
         # Extruded case: calculate the bounding boxes for all cells by running a kernel
         V = functionspace.VectorFunctionSpace(mesh, "DG", 0, dim=self.geometric_dimension)
@@ -2558,8 +2561,8 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         column_list = V.cell_node_list.reshape(-1)
         coords_min = mesh._order_data_by_cell_index(column_list, coords_min.dat.data_ro_with_halos)
         coords_max = mesh._order_data_by_cell_index(column_list, coords_max.dat.data_ro_with_halos)
-
-        return coords_min, coords_max
+        self._bounding_box_coords = coords_min, coords_max
+        return self._bounding_box_coords
 
     @property
     @PETSc.Log.EventDecorator()

From 74841c6cbdb080f9d8bc4326fcf79618aa8b9abc Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Mon, 15 Jun 2026 14:16:55 +0100
Subject: [PATCH 02/15] add helper function

---
 firedrake/mesh.py | 30 +++++++++++++++++++-----------
 1 file changed, 19 insertions(+), 11 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index a78e8296bb..d446bd9b27 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2479,7 +2479,17 @@ def clear_rtree(self):
         the coordinate field)."""
         self._rtree = None
 
-    @cached_property
+    def _clear_caches(self):
+        self._bounding_box_coords = None
+        self.clear_rtree()
+
+    def _check_coordinate_dat_version(self):
+        current = self.coordinates.dat.dat_version
+        if current != self._saved_coordinate_dat_version:
+            self._clear_caches()
+            self._saved_coordinate_dat_version = current
+
+    @property
     @PETSc.Log.EventDecorator()
     def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         """Calculates bounding boxes for the mesh rtree.
@@ -2496,8 +2506,10 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         Bezier curves and are completely contained in the convex hull of the mesh nodes.
         Hence the bounding box will contain the entire element.
         """
-        if self._bounding_box_coords:
+        self._check_coordinate_dat_version()
+        if self._bounding_box_coords is not None:
             return self._bounding_box_coords
+
         from firedrake import function, functionspace
         from firedrake.parloops import par_loop, READ, MIN, MAX
 
@@ -2532,7 +2544,7 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         cell_node_list = mesh.coordinates.function_space().cell_node_list
         if not mesh.extruded:
             all_coords = coords.dat.data_ro_with_halos[cell_node_list]
-            self._bounding_box_coords = np.min(all_coords, axis=1), np.max(all_coords, axis=1)
+            self._bounding_box_coords = (np.min(all_coords, axis=1), np.max(all_coords, axis=1))
             return self._bounding_box_coords
 
         # Extruded case: calculate the bounding boxes for all cells by running a kernel
@@ -2561,7 +2573,7 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         column_list = V.cell_node_list.reshape(-1)
         coords_min = mesh._order_data_by_cell_index(column_list, coords_min.dat.data_ro_with_halos)
         coords_max = mesh._order_data_by_cell_index(column_list, coords_max.dat.data_ro_with_halos)
-        self._bounding_box_coords = coords_min, coords_max
+        self._bounding_box_coords = (coords_min, coords_max)
         return self._bounding_box_coords
 
     @property
@@ -2582,12 +2594,9 @@ def rtree(self):
         can be found.
 
         """
-        if self.coordinates.dat.dat_version != self._saved_coordinate_dat_version:
-            if "bounding_box_coords" in self.__dict__:
-                del self.bounding_box_coords
-        else:
-            if self._rtree:
-                return self._rtree
+        self._check_coordinate_dat_version()
+        if self._rtree:
+            return self._rtree
         # Change min and max to refer to an n-hypercube, where n is the
         # geometric dimension of the mesh, centred on the midpoint of the
         # bounding box. Its side length is the L1 diameter of the bounding box.
@@ -2612,7 +2621,6 @@ def rtree(self):
 
         with PETSc.Log.Event("rtree_build"):
             self._rtree = rtree.build_from_aabb(coords_min, coords_max)
-        self._saved_coordinate_dat_version = self.coordinates.dat.dat_version
         return self._rtree
 
     @PETSc.Log.EventDecorator()

From 78e732c9361670b5210b74b90458bf2ee24092bd Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Mon, 15 Jun 2026 14:16:58 +0100
Subject: [PATCH 03/15] add test

---
 tests/firedrake/regression/test_point_eval_fs.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/tests/firedrake/regression/test_point_eval_fs.py b/tests/firedrake/regression/test_point_eval_fs.py
index 2a4a2c718a..d10051e0ac 100644
--- a/tests/firedrake/regression/test_point_eval_fs.py
+++ b/tests/firedrake/regression/test_point_eval_fs.py
@@ -243,3 +243,11 @@ def test_changing_coordinates_invalidates_rtree():
     saved_rtree = mesh.rtree
     mesh.coordinates.assign(mesh.coordinates * 2)
     assert mesh.rtree != saved_rtree
+
+
+def test_changing_coordinates_invalidates_bounding_box():
+    mesh = UnitSquareMesh(2, 2)
+
+    saved_bounding_box_coords = mesh.bounding_box_coords
+    mesh.coordinates.assign(mesh.coordinates * 2)
+    assert not np.allclose(mesh.bounding_box_coords, saved_bounding_box_coords)

From 12e44addfb69d28fc8c6ad8640bab6ebb6ff998c Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 16 Jun 2026 12:48:01 +0100
Subject: [PATCH 04/15] tidy imports

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

diff --git a/firedrake/function.py b/firedrake/function.py
index 9d8a219fb7..efc3578e5b 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -20,11 +20,10 @@
 from pyop2.exceptions import DataTypeError, DataValueError
 
 from finat.ufl import MixedElement
-from firedrake.utils import ScalarType, IntType, as_ctypes
+from firedrake.utils import ScalarType, IntType, as_ctypes, cached_property_until, _new_uid, complex_mode, ScalarType_c
 
 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
 from firedrake.mesh import MeshGeometry, VertexOnlyMesh
@@ -75,7 +74,7 @@ def __init__(self, function_space, val=None, name=None, dtype=ScalarType):
         # User comm
         self.comm = function_space.comm
         self._function_space = function_space
-        self.uid = utils._new_uid(self.comm)
+        self.uid = _new_uid(self.comm)
         self._name = name or 'function_%d' % self.uid
         self._label = "a function"
 
@@ -564,7 +563,7 @@ def evaluate(self, coord, mapping, component, index_values):
         # Called by UFL when evaluating expressions at coordinates
         if component or index_values:
             raise NotImplementedError("Unsupported arguments when attempting to evaluate Function.")
-        coord = np.asarray(coord, dtype=utils.ScalarType)
+        coord = np.asarray(coord, dtype=ScalarType)
         evaluator = PointEvaluator(self.function_space().mesh(), coord)
         result = evaluator.evaluate(self)
         if len(coord.shape) == 1:
@@ -602,8 +601,8 @@ def _at(self, arg, *args, **kwargs):
 
         if args:
             arg = (arg,) + args
-        arg = np.asarray(arg, dtype=utils.ScalarType)
-        if utils.complex_mode:
+        arg = np.asarray(arg, dtype=ScalarType)
+        if complex_mode:
             if not np.allclose(arg.imag, 0):
                 raise ValueError("Provided points have non-zero imaginary part")
             arg = arg.real.copy()
@@ -734,7 +733,7 @@ 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.points = np.asarray(points, dtype=utils.ScalarType)
+        self.points = np.asarray(points, dtype=ScalarType)
         if not self.points.shape:
             self.points = self.points.reshape(-1)
         gdim = mesh.geometric_dimension
@@ -823,7 +822,7 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         # 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)
+                result = np.empty((len(self.points),) + shape, dtype=ScalarType)
             self.mesh.comm.Bcast(result)
         return result
 
@@ -851,7 +850,7 @@ def make_c_evaluate(function, c_name="evaluate", ldargs=None, tolerance=None):
     arg = function.dat(op2.READ, function.cell_node_map())
     args.append(arg)
 
-    p_ScalarType_c = f"{utils.ScalarType_c}*"
+    p_ScalarType_c = f"{ScalarType_c}*"
     src.append(generate_single_cell_wrapper(mesh.cell_set, args,
                                             forward_args=[p_ScalarType_c,
                                                           p_ScalarType_c],

From 29ccc9919c390779b5bfd439d3d845075af6fff6 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 16 Jun 2026 12:48:35 +0100
Subject: [PATCH 05/15] add decorator, update mesh.py

---
 firedrake/mesh.py  | 27 +++------------------------
 firedrake/utils.py | 19 +++++++++++++++++++
 2 files changed, 22 insertions(+), 24 deletions(-)

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index d446bd9b27..b3e1899f23 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -37,7 +37,7 @@
 import firedrake.extrusion_utils as eutils
 import firedrake.cython.rtree as rtree
 import firedrake.utils as utils
-from firedrake.utils import as_cstr, IntType, RealType
+from firedrake.utils import as_cstr, IntType, RealType, cached_property_until
 from firedrake.logging import logger
 from firedrake.parameters import parameters
 from firedrake.petsc import PETSc, DEFAULT_PARTITIONER
@@ -2380,10 +2380,6 @@ def __init__(self, coordinates):
         # submesh
         self.submesh_parent = None
 
-        self._bounding_box_coords = None
-        self._rtree = None
-        self._saved_coordinate_dat_version = coordinates.dat.dat_version
-
         # Cache mesh object on the coordinateless coordinates function
         coordinates._as_mesh_geometry = weakref.ref(self)
 
@@ -2479,17 +2475,7 @@ def clear_rtree(self):
         the coordinate field)."""
         self._rtree = None
 
-    def _clear_caches(self):
-        self._bounding_box_coords = None
-        self.clear_rtree()
-
-    def _check_coordinate_dat_version(self):
-        current = self.coordinates.dat.dat_version
-        if current != self._saved_coordinate_dat_version:
-            self._clear_caches()
-            self._saved_coordinate_dat_version = current
-
-    @property
+    @cached_property_until(lambda self: self.coordinates.dat.dat_version)
     @PETSc.Log.EventDecorator()
     def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         """Calculates bounding boxes for the mesh rtree.
@@ -2506,10 +2492,6 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         Bezier curves and are completely contained in the convex hull of the mesh nodes.
         Hence the bounding box will contain the entire element.
         """
-        self._check_coordinate_dat_version()
-        if self._bounding_box_coords is not None:
-            return self._bounding_box_coords
-
         from firedrake import function, functionspace
         from firedrake.parloops import par_loop, READ, MIN, MAX
 
@@ -2576,7 +2558,7 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         self._bounding_box_coords = (coords_min, coords_max)
         return self._bounding_box_coords
 
-    @property
+    @cached_property_until(lambda self: self.coordinates.dat.dat_version)
     @PETSc.Log.EventDecorator()
     def rtree(self):
         """Builds an rtree from bounding box coordinates, expanding
@@ -2594,9 +2576,6 @@ def rtree(self):
         can be found.
 
         """
-        self._check_coordinate_dat_version()
-        if self._rtree:
-            return self._rtree
         # Change min and max to refer to an n-hypercube, where n is the
         # geometric dimension of the mesh, centred on the midpoint of the
         # bounding box. Its side length is the L1 diameter of the bounding box.
diff --git a/firedrake/utils.py b/firedrake/utils.py
index fa58f1f7a0..92a08be315 100644
--- a/firedrake/utils.py
+++ b/firedrake/utils.py
@@ -1,6 +1,8 @@
 # Some generic python utilities not really specific to our work.
 import collections.abc
 import warnings
+import functools
+from typing import Callable, Self, Hashable
 from decorator import decorator
 from pyop2.datatypes import ScalarType, as_cstr
 from pyop2.datatypes import RealType     # noqa: F401
@@ -264,3 +266,20 @@ def check_netgen_installed() -> None:
             "are installed and available to Firedrake (see "
             "https://www.firedrakeproject.org/install.html#netgen)."
         )
+
+
+def cached_property_until(value: Callable[[Self], Hashable]):
+    def decorator(func):
+        @property
+        @functools.wraps(func)
+        def wrapper(self):
+            value_attribute = f"_{func.__name__}_cache"
+            current_value = value(self)
+            cached_value = getattr(self, value_attribute, None)
+            if cached_value is None or cached_value[0] != current_value:
+                result = func(self)
+                setattr(self, value_attribute, (current_value, result))
+                return result
+            return cached_value[1]
+        return wrapper
+    return decorator
\ No newline at end of file

From 6aad43834e5d79db552b6a138e7834a495cb18ef Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 16 Jun 2026 13:12:18 +0100
Subject: [PATCH 06/15] tidy

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

diff --git a/firedrake/utils.py b/firedrake/utils.py
index 92a08be315..d59df59a0d 100644
--- a/firedrake/utils.py
+++ b/firedrake/utils.py
@@ -269,16 +269,19 @@ def check_netgen_installed() -> None:
 
 
 def cached_property_until(value: Callable[[Self], Hashable]):
+    """Decorator for a property that is cached until some value changes.
+    For examples of usage, see PointEvaluator and MeshGeometry.
+    """
     def decorator(func):
         @property
         @functools.wraps(func)
         def wrapper(self):
-            value_attribute = f"_{func.__name__}_cache"
+            cache_attribute = f"_{func.__name__}_cache"
             current_value = value(self)
-            cached_value = getattr(self, value_attribute, None)
+            cached_value = getattr(self, cache_attribute, None)
             if cached_value is None or cached_value[0] != current_value:
                 result = func(self)
-                setattr(self, value_attribute, (current_value, result))
+                setattr(self, cache_attribute, (current_value, result))
                 return result
             return cached_value[1]
         return wrapper

From bd6330067c31f043461c0da971e108352ed47474 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 16 Jun 2026 13:13:07 +0100
Subject: [PATCH 07/15] use for PointEvaluator

---
 firedrake/function.py | 34 +++++++++++++++-------------------
 1 file changed, 15 insertions(+), 19 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index efc3578e5b..3744c0bb0b 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -742,13 +742,19 @@ 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
-        self.vom = VertexOnlyMesh(
-            mesh, self.points, missing_points_behaviour=missing_points_behaviour,
-            redundant=redundant, tolerance=tolerance
+        if tolerance is not None:
+            mesh.tolerance = tolerance
+
+    @cached_property_until(
+        lambda self: (self.mesh.coordinates.dat.dat_version, self.mesh.tolerance)
+    )
+    def vom(self) -> MeshGeometry:
+        """The VOM used for point evaluation. This is cached until the mesh coordinates or tolerance change."""
+        return VertexOnlyMesh(
+            self.mesh, self.points, missing_points_behaviour=self.missing_points_behaviour,
+            redundant=self.redundant, tolerance=None
         )
 
     def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
@@ -786,19 +792,9 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         if function.function_space().ufl_element().family() == "Real":
             return function.dat.data_ro
 
-        function_mesh = function.function_space().mesh().unique()
-        if function_mesh is not self.mesh:
+        if function.function_space().mesh().unique() 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
-        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
-            )
+        vom = self.vom
 
         subfunctions = function.subfunctions
         if len(subfunctions) > 1:
@@ -812,8 +808,8 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         else:
             fs = partial(TensorFunctionSpace, shape=shape)
 
-        P0DG = fs(self.vom, "DG", 0)
-        P0DG_io = fs(self.vom.input_ordering, "DG", 0)
+        P0DG = fs(vom, "DG", 0)
+        P0DG_io = fs(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)

From 534b1fd5654c17107a2202ee688cf118161b33fa Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 16 Jun 2026 13:13:41 +0100
Subject: [PATCH 08/15] lint

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

diff --git a/firedrake/utils.py b/firedrake/utils.py
index d59df59a0d..12df341511 100644
--- a/firedrake/utils.py
+++ b/firedrake/utils.py
@@ -285,4 +285,4 @@ def wrapper(self):
                 return result
             return cached_value[1]
         return wrapper
-    return decorator
\ No newline at end of file
+    return decorator

From b0d466c46028719a255be6d17f483a2c1354b190 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 16 Jun 2026 13:18:28 +0100
Subject: [PATCH 09/15] tidy

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

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index b3e1899f23..d4d8db8c59 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2526,8 +2526,7 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         cell_node_list = mesh.coordinates.function_space().cell_node_list
         if not mesh.extruded:
             all_coords = coords.dat.data_ro_with_halos[cell_node_list]
-            self._bounding_box_coords = (np.min(all_coords, axis=1), np.max(all_coords, axis=1))
-            return self._bounding_box_coords
+            return np.min(all_coords, axis=1), np.max(all_coords, axis=1)
 
         # Extruded case: calculate the bounding boxes for all cells by running a kernel
         V = functionspace.VectorFunctionSpace(mesh, "DG", 0, dim=self.geometric_dimension)
@@ -2555,8 +2554,7 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         column_list = V.cell_node_list.reshape(-1)
         coords_min = mesh._order_data_by_cell_index(column_list, coords_min.dat.data_ro_with_halos)
         coords_max = mesh._order_data_by_cell_index(column_list, coords_max.dat.data_ro_with_halos)
-        self._bounding_box_coords = (coords_min, coords_max)
-        return self._bounding_box_coords
+        return coords_min, coords_max
 
     @cached_property_until(lambda self: self.coordinates.dat.dat_version)
     @PETSc.Log.EventDecorator()

From 54fab4773c9caba66221d03a315cc40c70eeefd0 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 23 Jun 2026 11:12:56 +0100
Subject: [PATCH 10/15] fixes

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

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index d4d8db8c59..757092cf34 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2473,7 +2473,7 @@ def clear_rtree(self):
 
         Use this if you move the mesh (for example by reassigning to
         the coordinate field)."""
-        self._rtree = None
+        self._rtree_cache = None
 
     @cached_property_until(lambda self: self.coordinates.dat.dat_version)
     @PETSc.Log.EventDecorator()
@@ -2556,7 +2556,7 @@ def bounding_box_coords(self) -> Tuple[np.ndarray, np.ndarray]:
         coords_max = mesh._order_data_by_cell_index(column_list, coords_max.dat.data_ro_with_halos)
         return coords_min, coords_max
 
-    @cached_property_until(lambda self: self.coordinates.dat.dat_version)
+    @cached_property_until(lambda self: (self.coordinates.dat.dat_version, self.tolerance))
     @PETSc.Log.EventDecorator()
     def rtree(self):
         """Builds an rtree from bounding box coordinates, expanding

From aa43b3b478574c5894d08bf82b1099681fffa531 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Wed, 24 Jun 2026 11:41:44 +0100
Subject: [PATCH 11/15] review suggestions

---
 firedrake/function.py                         |  6 +++---
 firedrake/utils.py                            | 20 ++++++++++++++++---
 .../regression/test_point_eval_fs.py          | 12 +++++++++--
 3 files changed, 30 insertions(+), 8 deletions(-)

diff --git a/firedrake/function.py b/firedrake/function.py
index 3744c0bb0b..886a4db5ac 100644
--- a/firedrake/function.py
+++ b/firedrake/function.py
@@ -752,6 +752,7 @@ def __init__(self, mesh: MeshGeometry, points: np.ndarray | list, tolerance: flo
     )
     def vom(self) -> MeshGeometry:
         """The VOM used for point evaluation. This is cached until the mesh coordinates or tolerance change."""
+        # We invalidate when the parent mesh moves until https://github.com/firedrakeproject/firedrake/issues/4540 is fixed.
         return VertexOnlyMesh(
             self.mesh, self.points, missing_points_behaviour=self.missing_points_behaviour,
             redundant=self.redundant, tolerance=None
@@ -794,7 +795,6 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
 
         if function.function_space().mesh().unique() is not self.mesh:
             raise ValueError("Function mesh must be the same Mesh object as the PointEvaluator mesh.")
-        vom = self.vom
 
         subfunctions = function.subfunctions
         if len(subfunctions) > 1:
@@ -808,8 +808,8 @@ def evaluate(self, function: Function) -> np.ndarray | Tuple[np.ndarray, ...]:
         else:
             fs = partial(TensorFunctionSpace, shape=shape)
 
-        P0DG = fs(vom, "DG", 0)
-        P0DG_io = fs(vom.input_ordering, "DG", 0)
+        P0DG = fs(self.vom, "DG", 0)
+        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)
diff --git a/firedrake/utils.py b/firedrake/utils.py
index 12df341511..b3103e11c3 100644
--- a/firedrake/utils.py
+++ b/firedrake/utils.py
@@ -268,16 +268,30 @@ def check_netgen_installed() -> None:
         )
 
 
-def cached_property_until(value: Callable[[Self], Hashable]):
+def cached_property_until(key: Callable[[Self], Hashable]):
     """Decorator for a property that is cached until some value changes.
-    For examples of usage, see PointEvaluator and MeshGeometry.
+    
+    For example, the `expensive_property` below will be cached until `self.value` changes,
+    and will be recomputed with the new `self.value` when accessed again.
+
+    ```
+    class MyClass:
+        
+        def __init__(self, value):
+            self.value = value
+
+        @cached_property_until(lambda self: self.value)
+        def expensive_property(self):
+            # Some expensive computation that depends on self.value
+            ...
+    ```
     """
     def decorator(func):
         @property
         @functools.wraps(func)
         def wrapper(self):
             cache_attribute = f"_{func.__name__}_cache"
-            current_value = value(self)
+            current_value = key(self)
             cached_value = getattr(self, cache_attribute, None)
             if cached_value is None or cached_value[0] != current_value:
                 result = func(self)
diff --git a/tests/firedrake/regression/test_point_eval_fs.py b/tests/firedrake/regression/test_point_eval_fs.py
index d10051e0ac..e2ba9553ce 100644
--- a/tests/firedrake/regression/test_point_eval_fs.py
+++ b/tests/firedrake/regression/test_point_eval_fs.py
@@ -239,7 +239,7 @@ def test_point_reset_works():
 
 def test_changing_coordinates_invalidates_rtree():
     mesh = UnitSquareMesh(2, 2)
-
+    assert mesh.rtree is mesh.rtree
     saved_rtree = mesh.rtree
     mesh.coordinates.assign(mesh.coordinates * 2)
     assert mesh.rtree != saved_rtree
@@ -247,7 +247,15 @@ def test_changing_coordinates_invalidates_rtree():
 
 def test_changing_coordinates_invalidates_bounding_box():
     mesh = UnitSquareMesh(2, 2)
-
+    assert mesh.bounding_box_coords is mesh.bounding_box_coords
     saved_bounding_box_coords = mesh.bounding_box_coords
     mesh.coordinates.assign(mesh.coordinates * 2)
     assert not np.allclose(mesh.bounding_box_coords, saved_bounding_box_coords)
+
+
+def test_changing_tolerance_invalidates_rtree():
+    mesh = UnitSquareMesh(2, 2)
+    assert mesh.rtree is mesh.rtree
+    saved_rtree = mesh.rtree
+    mesh.tolerance = 1e-5
+    assert mesh.rtree != saved_rtree

From 0f9e94f965edf0cf1fb839778f9df8362156743d Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Wed, 24 Jun 2026 11:43:24 +0100
Subject: [PATCH 12/15] wording

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

diff --git a/firedrake/utils.py b/firedrake/utils.py
index b3103e11c3..be6d48dc37 100644
--- a/firedrake/utils.py
+++ b/firedrake/utils.py
@@ -272,7 +272,7 @@ def cached_property_until(key: Callable[[Self], Hashable]):
     """Decorator for a property that is cached until some value changes.
     
     For example, the `expensive_property` below will be cached until `self.value` changes,
-    and will be recomputed with the new `self.value` when accessed again.
+    and will be recomputed with the new `self.value` and cached when accessed again.
 
     ```
     class MyClass:

From 4b3f3a3c15bc8517f5fa17ee14d7a89d35cf401d Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Wed, 24 Jun 2026 11:43:47 +0100
Subject: [PATCH 13/15] lint

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

diff --git a/firedrake/utils.py b/firedrake/utils.py
index be6d48dc37..46b59099e6 100644
--- a/firedrake/utils.py
+++ b/firedrake/utils.py
@@ -270,13 +270,13 @@ def check_netgen_installed() -> None:
 
 def cached_property_until(key: Callable[[Self], Hashable]):
     """Decorator for a property that is cached until some value changes.
-    
+
     For example, the `expensive_property` below will be cached until `self.value` changes,
     and will be recomputed with the new `self.value` and cached when accessed again.
 
     ```
     class MyClass:
-        
+
         def __init__(self, value):
             self.value = value
 

From 7441d2dca5ac74e2fe73f811a2b19f34809a56a1 Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 7 Jul 2026 10:45:52 +0100
Subject: [PATCH 14/15] fix docstring

---
 firedrake/utils.py | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/firedrake/utils.py b/firedrake/utils.py
index 46b59099e6..068e557bfc 100644
--- a/firedrake/utils.py
+++ b/firedrake/utils.py
@@ -271,20 +271,21 @@ def check_netgen_installed() -> None:
 def cached_property_until(key: Callable[[Self], Hashable]):
     """Decorator for a property that is cached until some value changes.
 
-    For example, the `expensive_property` below will be cached until `self.value` changes,
-    and will be recomputed with the new `self.value` and cached when accessed again.
+    For example, the ``expensive_property`` below will be cached until
+    ``self.value`` changes, and will be recomputed with the new ``self.value``
+    and cached when accessed again.
 
-    ```
-    class MyClass:
+    .. code-block:: python
 
-        def __init__(self, value):
-            self.value = value
+        class MyClass:
 
-        @cached_property_until(lambda self: self.value)
-        def expensive_property(self):
-            # Some expensive computation that depends on self.value
-            ...
-    ```
+            def __init__(self, value):
+                self.value = value
+
+            @cached_property_until(lambda self: self.value)
+            def expensive_property(self):
+                # Some expensive computation that depends on self.value
+                ...
     """
     def decorator(func):
         @property

From d3bc21cc6c99bcebff2682b75201965f9317e0cb Mon Sep 17 00:00:00 2001
From: Leo Collins <collins.leo@proton.me>
Date: Tue, 7 Jul 2026 13:03:38 +0100
Subject: [PATCH 15/15] add comment on `clear_rtree`

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

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index 757092cf34..57dc5771c5 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2473,6 +2473,8 @@ def clear_rtree(self):
 
         Use this if you move the mesh (for example by reassigning to
         the coordinate field)."""
+        # `cached_property_until` stores the cached rtree in self._rtree_cache
+        # setting it to None will force the rtree to be rebuilt on next access.
         self._rtree_cache = None
 
     @cached_property_until(lambda self: self.coordinates.dat.dat_version)