Place a tip vertex from the face that holds it, not from the mesh at large

src/underworld3/utilities/line_cut.py

@@ -379,6 +379,53 @@ def _cell_vertices(dm):
                     dtype=np.int64).reshape(cE - cS, 3)
 
 
+def _minimal_face_vertices(X, cell_verts, point):
+    """Vertices of the smallest mesh face whose closure contains ``point``.
+
+    A point sits in the relative interior of exactly one face: a triangle, an
+    edge, or a vertex. This returns that face's vertices — the three corners
+    of the cell holding it, the two ends of the edge it lies on, or the single
+    vertex it coincides with. Empty when the point lies outside the mesh.
+
+    The face is read off the BARYCENTRIC coordinates of the containing cell: a
+    coordinate is zero exactly when the point lies on the opposite edge, so
+    keeping the corners with a nonzero coordinate drops the corners the point
+    is not "on". Doing it this way means one cell answers the question. Taking
+    the union of the corners of every containing cell does NOT: a point on an
+    edge lies in two cells, and the union then offers each cell's opposite
+    apex, which belongs to one of them and not the other. Moving such an apex
+    leaves the other cell without a corner on the point, which is the defect
+    this whole path exists to avoid (issue #542).
+
+    The coordinates are scale-free, so one relative tolerance serves cells of
+    any size — which embedded meshes need, their cell sizes differing by an
+    order of magnitude either side of a zone skin.
+    """
+    A = X[cell_verts[:, 0], :2]
+    B = X[cell_verts[:, 1], :2]
+    C = X[cell_verts[:, 2], :2]
+    v0, v1, v2 = C - A, B - A, np.asarray(point, dtype=float)[:2] - A
+    den = v0[:, 0] * v1[:, 1] - v1[:, 0] * v0[:, 1]
+    # A degenerate cell has no interior to be inside of; it cannot hold the
+    # point and must not divide.
+    alive = np.abs(den) > 0.0
+    den = np.where(alive, den, 1.0)
+    lc = (v2[:, 0] * v1[:, 1] - v1[:, 0] * v2[:, 1]) / den
+    lb = (v0[:, 0] * v2[:, 1] - v2[:, 0] * v0[:, 1]) / den
+    la = 1.0 - lb - lc
+
+    tol = 1e-9
+    inside = np.flatnonzero(alive & (la >= -tol) & (lb >= -tol) & (lc >= -tol))
+    if not len(inside):
+        return np.empty(0, dtype=np.int64)
+    # One containing cell settles it; the first keeps the answer independent
+    # of how the cells happen to be ordered, since every containing cell names
+    # the same face.
+    c = int(inside[0])
+    on = np.array([la[c] > tol, lb[c] > tol, lc[c] > tol])
+    return cell_verts[c][on]
+
+
 def _cell_quality(X, cell_verts):
     """Scale-free triangle quality: 1 equilateral, 0 degenerate, <0 inverted.
 
@@ -839,6 +886,14 @@ def pull_vertex_onto(dm, targets):
     behaviour as ``snap_frac``, and anything optimising over fault geometry has
     to live with it.
 
+    The vertex is chosen from the corners of the cells CONTAINING the target,
+    not from the mesh at large. The globally nearest vertex can belong to a
+    neighbouring cell, and moving that one leaves the containing cell with no
+    corner on the line, so :func:`cut_along_lines` refuses the very tip this
+    call was made to place (issue #542). Restricting the candidates ties the
+    choice to the containment relation the cut actually tests. A target outside
+    the mesh falls back to the nearest vertex anywhere.
+
     The choice is made from the coordinates alone and reduced globally, so every
     rank moves the same vertex: a rank-local nearest-vertex search picks a
     different one on each rank, which is how the mesh stops being
@@ -849,15 +904,38 @@ def pull_vertex_onto(dm, targets):
     X = _coords(dm)
     arr = X.copy()
     scale = _global_extent(dm)
+    cell_verts = _cell_vertices(dm)
 
     for t in np.atleast_2d(np.asarray(targets, dtype=float))[:, :2]:
-        d = np.linalg.norm(X[:, :2] - t, axis=1)
-        # Reduced as (distance, x, y) so the tie-break is part of the same
-        # reduction: tuples compare lexicographically, so MIN gives the closest
-        # vertex and, among equals, the one lowest in coordinate order.
-        local = ((float(d.min()), *X[int(d.argmin()), :2]) if d.size
-                 else (np.inf, np.inf, np.inf))
-        _dist, tx, ty = uw.mpi.comm.allreduce(local, op=MPI.MIN)
+        # Candidates are the vertices of the smallest face CONTAINING the
+        # target, not every vertex in the mesh. The nearest vertex overall can
+        # belong to a neighbouring cell, and moving it leaves the containing
+        # cell without a corner on the line — precisely the "entered but not
+        # left" state cut_along_lines refuses, reported for a tip that does
+        # have a vertex on it (issue #542). Measured inside a ribbon zone of
+        # width 0.02: the target sat on an edge whose two ends were at exactly
+        # 1.000e-02, being the ribbon's own edge vertices, and the apex of one
+        # of the two cells sharing that edge won by 0.4% at 9.963e-03.
+        face = _minimal_face_vertices(X, cell_verts, t)
+        held = len(face) > 0
+        candidates = face if held else np.arange(len(X))
+        d = np.linalg.norm(X[candidates, :2] - t, axis=1)
+        # Reduced as (outside, distance, x, y) so both decisions ride one
+        # reduction: tuples compare lexicographically, so MIN prefers ANY rank
+        # that holds a containing cell over every rank that does not, then the
+        # closest corner, then the one lowest in coordinate order. The
+        # outside-flag has to lead: a rank with no part of the containing cell
+        # proposes its own nearest vertex, which can easily be closer than the
+        # right answer, and a reduction on distance alone would take it.
+        #
+        # All ranks flagging 1 means the target lies outside the mesh, and the
+        # fallback is the nearest vertex anywhere — the behaviour before the
+        # containment rule, kept because a target off the mesh still has to
+        # land somewhere.
+        local = ((0 if held else 1, float(d.min()),
+                  *X[candidates[int(d.argmin())], :2])
+                 if d.size else (1, np.inf, np.inf, np.inf))
+        _outside, _dist, tx, ty = uw.mpi.comm.allreduce(local, op=MPI.MIN)
 
         # Move it by POSITION, not by index: the chosen vertex may be a ghost
         # here and an owned point there, and both copies have to end up in the

tests/test_0844_line_cut.py

@@ -879,3 +879,74 @@ def test_repair_fixes_the_cells_the_cut_leaves_thin():
     assert cell_areas(fixed.dm).sum() == pytest.approx(
         cell_areas(plain.dm).sum(), rel=1e-13)
     assert (cell_areas(fixed.dm) > 0).all()
+
+
+def test_the_pulled_vertex_covers_every_cell_holding_the_target():
+    """A target ON AN EDGE must take one of that edge's ends, not an apex.
+
+    `pull_vertex_onto` used to move the nearest vertex in the mesh at large.
+    That vertex can be the apex of one of the two cells sharing the edge the
+    target lies on: it wins on distance, it belongs to one of the two cells,
+    and the OTHER cell is left with no corner on the target. A line ending
+    there is then "entered but not left" and `cut_along_lines` refuses a tip
+    that does have a vertex on it — issue #542, measured inside a ribbon zone
+    where the edge ends sit at exactly the zone half-width and an apex won by
+    0.4%.
+
+    The invariant asserted is the one the cut actually needs: EVERY cell whose
+    closure contains the target ends up with a corner exactly on it. Edge
+    midpoints are used because that is the configuration that fails.
+
+    The mesh is FLATTENED in y first, and that is not cosmetic: on a
+    well-shaped mesh an edge midpoint's nearest vertex is one of that edge's
+    own ends, so the defect never appears and the test would pass without
+    testing anything. Flat cells put an apex within reach of the opposite
+    edge, which is the shape a ribbon zone has along its length. The negative
+    control below asserts the mesh really does present the case.
+    """
+    from underworld3.utilities.line_cut import (_cell_vertices, _coords,
+                                                _minimal_face_vertices,
+                                                _set_coordinates)
+
+    base = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1 / 12,
+        regular=False, qdegree=2)
+    dm = base.dm.clone()
+    squashed = _coords(dm).copy()
+    squashed[:, 1] *= 0.12
+    _set_coordinates(dm, np.arange(len(squashed)), squashed)
+
+    X = _coords(dm)
+    cells = _cell_vertices(dm)
+
+    interior = [e for e in range(*dm.getDepthStratum(1))
+                if len(dm.getSupport(e)) == 2]
+    vS = dm.getDepthStratum(0)[0]
+    midpoints = []
+    for e in interior[::7]:
+        a, b = (int(v) - vS for v in dm.getCone(e))
+        midpoints.append(0.5 * (X[a, :2] + X[b, :2]))
+    assert len(midpoints) > 5
+
+    exposed = 0
+    for t in midpoints:
+        face = _minimal_face_vertices(X, cells, t)
+        assert len(face) == 2, "an edge midpoint lies on an edge, two ends"
+        nearest = int(np.linalg.norm(X[:, :2] - t, axis=1).argmin())
+        if nearest not in face:
+            exposed += 1
+
+        out = pull_vertex_onto(dm, t)
+        Xo = _coords(out)
+        holders = [c for c in range(len(cells))
+                   if len(_minimal_face_vertices(X, cells[c: c + 1], t))]
+        assert holders, "the midpoint is inside the mesh"
+        for c in holders:
+            on = np.linalg.norm(Xo[cells[c], :2] - t, axis=1).min()
+            assert on < 1e-12, (
+                f"cell {c} holds the target and has no corner on it "
+                f"(closest {on:.3e}); a line ending here would be refused")
+
+    assert exposed > 0, (
+        "no midpoint had its nearest vertex off the edge, so this mesh does "
+        "not exercise the defect and the test proves nothing")