fix(meshing): on_boundary kwarg for in-cell test — accept on-face queries by default

src/underworld3/discretisation/discretisation_mesh.py

@@ -3233,7 +3233,8 @@ def _mark_faces_inside_and_out(self):
 
         return
 
-    def _test_if_points_in_cells_internal(self, points, cells):
+    def _test_if_points_in_cells_internal(self, points, cells,
+                                          on_boundary=True, tol=0.0):
         """
         Determine if the given points lie in the suggested cells.
         Uses a mesh skeletonization array to determine whether the point is
@@ -3243,10 +3244,48 @@ def _test_if_points_in_cells_internal(self, points, cells):
 
         Parameters
         ----------
-        points : array-like
-            Coordinate array in any physical unit system (will be auto-converted)
-        cells : array-like
-            Cell indices to test
+        points : numpy.ndarray
+            Coordinate array, assumed already in model units (this internal
+            helper does not perform unit conversion — use the public
+            `test_if_points_in_cells` for unit-aware input).
+        cells : numpy.ndarray
+            1-D cell indices to test, one per point.
+        on_boundary : bool, default True
+            If True (the default), a point exactly on a cell face counts as
+            inside that cell — the natural semantics for FE evaluation,
+            where the basis at a shared face/vertex is consistent across
+            the adjacent cells. A query point lying on a face shared by N
+            cells passes the test for any of those N cells.
+
+            If False, a point exactly on a face is reported as NOT inside —
+            strict-inside semantics. Use this when uniqueness matters (a
+            strict-ownership scheme where a shared-face point being claimed
+            by all adjacent cells would be a bug).
+
+            The implementation compares the squared distance from the query
+            to a mirrored inner/outer control-point pair placed ±1e-3 along
+            the face normal; a point exactly on the face has zero distance
+            difference. With on_boundary=True the test accepts diff >= -1e-12
+            (well below the 1e-3 control-point offset, well above 64-bit
+            float roundoff); with on_boundary=False the test requires diff > 0.
+        tol : float, default 0.0
+            Face-RELATIVE tolerance — overrides ``on_boundary``'s absolute
+            -1e-12 floor when nonzero. The test becomes
+            ``diff > -tol * |O - I|^2`` (i.e. relative to the
+            control-point separation squared). With ``tol == 0`` (the
+            default) ``on_boundary`` controls the test as documented above.
+            With ``tol > 0`` the test is relaxed to admit points within
+            roughly ``tol`` of the face (relative to the face-normal
+            separation), while still rejecting points that lie inside a
+            *different* cell (whose diff is strongly negative). The
+            parallel evaluation locator
+            (``Mesh._robust_owning_cells``) uses ``tol=1e-2`` to admit
+            on-face / partition-seam node queries that PR #207's
+            ``on_boundary=True`` path (absolute 1e-12) is too tight to
+            accept — for query coordinates that are slightly off the
+            face (e.g. RBF-shifted node points), the geometric scale of
+            the test needs to match the *mesh* spacing, not float
+            roundoff. See memory project_pr207_loose_boundary_clash.
         """
         # Internal version - points assumed to already be in model units
         self._mark_faces_inside_and_out()
@@ -3257,18 +3296,40 @@ def _test_if_points_in_cells_internal(self, points, cells):
         cStart, cEnd = self.dm.getHeightStratum(0)
         num_cell_faces = self.dm.getConeSize(cStart)
 
-        inside = numpy.ones_like(cells, dtype=bool)
         insiders = numpy.ndarray(shape=(cells.shape[0], num_cell_faces), dtype=bool)
 
-        for f in range(num_cell_faces):
-            control_points_o = self.faces_outer_control_points[f, cells]
-            control_points_i = self.faces_inner_control_points[f, cells]
-            inside = (
-                ((control_points_o - points) ** 2).sum(axis=1)
-                - ((control_points_i - points) ** 2).sum(axis=1)
-            ) > 0
-
-            insiders[:, f] = inside[:]
+        if tol > 0.0:
+            # Face-relative tolerance branch (parallel evaluation locator).
+            # Takes precedence over on_boundary: a non-zero tol expresses
+            # a geometric tolerance on the face-normal separation²; the
+            # absolute on_boundary floor (~1e-12) is far below it and the
+            # strict on_boundary=False (>0) is what tol is widening.
+            for f in range(num_cell_faces):
+                control_points_o = self.faces_outer_control_points[f, cells]
+                control_points_i = self.faces_inner_control_points[f, cells]
+                value = (
+                    ((control_points_o - points) ** 2).sum(axis=1)
+                    - ((control_points_i - points) ** 2).sum(axis=1)
+                )
+                sep2 = ((control_points_o - control_points_i) ** 2).sum(axis=1)
+                insiders[:, f] = value > -tol * sep2
+        elif on_boundary:
+            _face_tol = -1e-12
+            for f in range(num_cell_faces):
+                control_points_o = self.faces_outer_control_points[f, cells]
+                control_points_i = self.faces_inner_control_points[f, cells]
+                insiders[:, f] = (
+                    ((control_points_o - points) ** 2).sum(axis=1)
+                    - ((control_points_i - points) ** 2).sum(axis=1)
+                ) >= _face_tol
+        else:
+            for f in range(num_cell_faces):
+                control_points_o = self.faces_outer_control_points[f, cells]
+                control_points_i = self.faces_inner_control_points[f, cells]
+                insiders[:, f] = (
+                    ((control_points_o - points) ** 2).sum(axis=1)
+                    - ((control_points_i - points) ** 2).sum(axis=1)
+                ) > 0
 
         return numpy.all(insiders, axis=1)
 
@@ -3469,7 +3530,12 @@ def get_closest_cells(self, coords: numpy.ndarray) -> numpy.ndarray:
             # CRITICAL: Must return 1D array, not 2D, for Cython buffer compatibility
             return numpy.array([], dtype=numpy.int64)
 
-    def _get_closest_local_cells_internal(self, coords: numpy.ndarray) -> numpy.ndarray:
+    def _get_closest_local_cells_internal(
+        self,
+        coords: numpy.ndarray,
+        on_boundary: bool = True,
+        tol: float = 0.0,
+    ) -> numpy.ndarray:
         """
         This method uses a kd-tree algorithm to find the closest
         cells to the provided coords. For a regular mesh, this should
@@ -3482,7 +3548,24 @@ def _get_closest_local_cells_internal(self, coords: numpy.ndarray) -> numpy.ndar
         coords:
             An array of the coordinates for which we wish to determine the
             closest cells. This should be a 2-dimensional array of
-            shape (n_coords,dim) in any physical unit system (will be auto-converted).
+            shape (n_coords,dim), assumed already in model units (this
+            internal helper does not perform unit conversion — use the
+            public `get_closest_local_cells` for unit-aware input).
+        on_boundary : bool, default True
+            Forwarded to `_test_if_points_in_cells_internal`. If True (the
+            default), queries exactly on a cell face count as inside that
+            cell — the natural semantics for FE-evaluation hints (every mesh
+            vertex sits on the faces of every cell containing it). If False,
+            strict-inside semantics; boundary queries come back as -1.
+        tol : float, default 0.0
+            Face-relative tolerance, forwarded to
+            `_test_if_points_in_cells_internal`. When > 0, the in-cell
+            test becomes ``diff > -tol * |O-I|²`` and takes precedence
+            over ``on_boundary``. Used by the parallel evaluation
+            locator (`Mesh._robust_owning_cells`) to admit on-face /
+            partition-seam node queries at a mesh-spacing-relative
+            tolerance — wider than ``on_boundary=True``'s absolute
+            -1e-12 floor.
 
         Returns:
         --------
@@ -3514,7 +3597,8 @@ def _get_closest_local_cells_internal(self, coords: numpy.ndarray) -> numpy.ndar
         cells = self._indexMap[closest_points]
         cStart, cEnd = self.dm.getHeightStratum(0)
 
-        inside = self._test_if_points_in_cells_internal(coords, cells)
+        inside = self._test_if_points_in_cells_internal(
+            coords, cells, on_boundary=on_boundary, tol=tol)
         cells[~inside] = -1
         lost_points = np.where(inside == False)[0]
 
@@ -3533,7 +3617,8 @@ def _get_closest_local_cells_internal(self, coords: numpy.ndarray) -> numpy.ndar
         for i in range(0, num_testable_neighbours):
 
             inside = self._test_if_points_in_cells_internal(
-                coords[lost_points], closest_centroids[:, i]
+                coords[lost_points], closest_centroids[:, i],
+                on_boundary=on_boundary, tol=tol,
             )
             cells[lost_points[inside]] = closest_centroids[inside, i]
 
@@ -3542,7 +3627,7 @@ def _get_closest_local_cells_internal(self, coords: numpy.ndarray) -> numpy.ndar
 
         return cells
 
-    def test_if_points_in_cells(self, points, cells):
+    def test_if_points_in_cells(self, points, cells, on_boundary=True, tol=0.0):
         """
         Determine if the given points lie in the suggested cells.
         Uses a mesh skeletonization array to determine whether the point is
@@ -3556,6 +3641,19 @@ def test_if_points_in_cells(self, points, cells):
             Coordinate array in any physical unit system (will be auto-converted)
         cells : array-like
             Cell indices to test
+        on_boundary : bool, default True
+            If True (the default), points exactly on a cell face count as
+            inside the cell (natural for FE evaluation, where the basis at
+            a shared face/vertex is consistent across adjacent cells). If
+            False, points on the closure of a cell are reported as NOT in
+            it (strict-inside semantics — useful when uniqueness matters).
+        tol : float, default 0.0
+            Face-relative tolerance forwarded to
+            `_test_if_points_in_cells_internal`. When ``> 0`` takes
+            precedence over ``on_boundary``: the test admits points
+            within ``tol`` of the face relative to the control-point
+            separation² — used by the parallel evaluation locator
+            for on-face / near-face queries at the mesh-spacing scale.
 
         Returns
         -------
@@ -3574,10 +3672,22 @@ def test_if_points_in_cells(self, points, cells):
         else:
             model_points = model_quantity
 
+        # Coerce cells to a 1-D numpy array — accept list/tuple input as the
+        # docstring promises ("array-like") even though the internal helper
+        # calls cells.reshape(-1) directly.
+        cells = numpy.asarray(cells).reshape(-1)
+
         # Call internal implementation
-        return self._test_if_points_in_cells_internal(model_points, cells)
+        return self._test_if_points_in_cells_internal(
+            model_points, cells, on_boundary=on_boundary, tol=tol,
+        )
 
-    def get_closest_local_cells(self, coords: numpy.ndarray) -> numpy.ndarray:
+    def get_closest_local_cells(
+        self,
+        coords: numpy.ndarray,
+        on_boundary: bool = True,
+        tol: float = 0.0,
+    ) -> numpy.ndarray:
         """
         This method uses a kd-tree algorithm to find the closest
         cells to the provided coords. For a regular mesh, this should
@@ -3591,6 +3701,18 @@ def get_closest_local_cells(self, coords: numpy.ndarray) -> numpy.ndarray:
             An array of the coordinates for which we wish to determine the
             closest cells. This should be a 2-dimensional array of
             shape (n_coords,dim) in any physical unit system (will be auto-converted).
+        on_boundary : bool, default True
+            If True (the default), queries exactly on a cell face are
+            treated as inside that cell (natural for FE-evaluation hints —
+            mesh vertices sit on cell faces by definition). If False,
+            strict-inside semantics; boundary queries return -1.
+        tol : float, default 0.0
+            Face-relative tolerance forwarded to
+            `_test_if_points_in_cells_internal`. When ``> 0`` takes
+            precedence over ``on_boundary``: the test admits points
+            within ``tol`` of the face relative to the control-point
+            separation² — used by the parallel evaluation locator
+            for on-face / near-face queries at the mesh-spacing scale.
 
         Returns:
         --------
@@ -3612,7 +3734,9 @@ def get_closest_local_cells(self, coords: numpy.ndarray) -> numpy.ndarray:
             model_coords = model_quantity
 
         # Call internal implementation
-        return self._get_closest_local_cells_internal(model_coords)
+        return self._get_closest_local_cells_internal(
+            model_coords, on_boundary=on_boundary, tol=tol,
+        )
 
     def _get_mesh_sizes(self, verbose=False):
         """

tests/test_0820_in_cell_test_loose_semantics.py

@@ -0,0 +1,121 @@
+"""Regression test for `Mesh._test_if_points_in_cells_internal` on_boundary modes.
+
+Locks in the contract for the `on_boundary` kwarg added to
+`_test_if_points_in_cells_internal` (and forwarded through
+`_get_closest_local_cells_internal`, `get_closest_local_cells`, and the
+public `test_if_points_in_cells`):
+
+- on_boundary=True (default): a point exactly on a cell face counts as
+  inside the cell — the natural semantics for FE evaluation, where the
+  basis at a shared face/vertex is consistent across the adjacent cells.
+- on_boundary=False: strict-inside semantics — a point on the face is
+  reported as NOT inside. Useful when uniqueness matters.
+"""
+
+import numpy as np
+import pytest
+
+import underworld3 as uw
+
+
+pytestmark = pytest.mark.level_1
+
+
+def test_default_accepts_vertices_simplex_2d():
+    """Default (on_boundary=True): every 2D simplex vertex resolves to a containing cell."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25
+    )
+    verts = np.asarray(mesh.X.coords)
+    cells = mesh._get_closest_local_cells_internal(verts)
+    assert (cells == -1).sum() == 0, (
+        f"default loose mode rejected {(cells == -1).sum()}/{len(verts)} vertices"
+    )
+
+
+def test_default_accepts_vertices_simplex_3d():
+    """Default (on_boundary=True): every 3D simplex vertex resolves to a containing cell."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), cellSize=0.4
+    )
+    verts = np.asarray(mesh.X.coords)
+    cells = mesh._get_closest_local_cells_internal(verts)
+    assert (cells == -1).sum() == 0, (
+        f"default loose mode rejected {(cells == -1).sum()}/{len(verts)} vertices"
+    )
+
+
+def test_default_accepts_vertices_quad():
+    """Default (on_boundary=True): every structured-quad vertex resolves to a containing cell."""
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)
+    )
+    verts = np.asarray(mesh.X.coords)
+    cells = mesh._get_closest_local_cells_internal(verts)
+    assert (cells == -1).sum() == 0, (
+        f"default loose mode rejected {(cells == -1).sum()}/{len(verts)} vertices"
+    )
+
+
+def test_on_boundary_false_rejects_vertices_simplex_3d():
+    """on_boundary=False reproduces strict-inside semantics — most 3D simplex vertices come back -1."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), cellSize=0.4
+    )
+    verts = np.asarray(mesh.X.coords)
+    cells = mesh._get_closest_local_cells_internal(verts, on_boundary=False)
+    assert (cells == -1).sum() > 0, (
+        "expected strict mode to reject at least some boundary-vertex queries"
+    )
+
+
+def test_on_boundary_modes_diverge_at_face_queries():
+    """Strict and loose mode must give a distinguishable result on vertex queries."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), cellSize=0.4
+    )
+    verts = np.asarray(mesh.X.coords)
+    hint = mesh.get_closest_cells(verts)
+    inside_strict = mesh._test_if_points_in_cells_internal(verts, hint, on_boundary=False)
+    inside_loose = mesh._test_if_points_in_cells_internal(verts, hint, on_boundary=True)
+    assert (~inside_strict).sum() > (~inside_loose).sum(), (
+        f"strict-vs-loose distinction lost: strict rejected {(~inside_strict).sum()}, "
+        f"loose rejected {(~inside_loose).sum()}"
+    )
+    assert (~inside_loose).sum() == 0, (
+        f"loose mode rejected {(~inside_loose).sum()} kdtree-nearest cells of vertices"
+    )
+
+
+def test_default_returns_containing_cells():
+    """The cell id returned by the default must have the query in its closure."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), cellSize=0.4
+    )
+    verts = np.asarray(mesh.X.coords)
+    cells = mesh._get_closest_local_cells_internal(verts)
+    assert (cells == -1).sum() == 0
+
+    cStart, _ = mesh.dm.getHeightStratum(0)
+    pStart, pEnd = mesh.dm.getDepthStratum(0)
+    for v, c in zip(verts, cells):
+        closure = mesh.dm.getTransitiveClosure(int(c) + cStart)[0]
+        vp = closure[(closure >= pStart) & (closure < pEnd)]
+        vc = mesh._coords[vp - pStart]
+        assert np.linalg.norm(vc - v, axis=1).min() < 1e-10, (
+            f"vertex {v} returned cell {c} whose closure does not contain it"
+        )
+
+
+def test_get_closest_local_cells_public_forwards_kwarg():
+    """The public `get_closest_local_cells` wrapper forwards on_boundary."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25
+    )
+    verts = np.asarray(mesh.X.coords)
+    # Default (True): no vertices rejected
+    cells_loose = mesh.get_closest_local_cells(verts)
+    # Opt out (False): expect some vertices rejected
+    cells_strict = mesh.get_closest_local_cells(verts, on_boundary=False)
+    assert (cells_loose == -1).sum() == 0
+    assert (cells_strict == -1).sum() > 0