From bcbce7c042e047763e68b01cce8f12be0f7701d9 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 13 Aug 2026 14:43:38 +1000 Subject: [PATCH] Place a tip vertex from the face that holds it, not from the mesh at large MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit pull_vertex_onto moved the nearest vertex in the whole mesh onto each target. That vertex need not belong to the cell containing the target: it can be the apex of a neighbouring cell, and moving it leaves the containing cell with no corner on the line. cut_along_lines then refuses the tip with "entered but not left" — for a tip that does have a vertex on it. The candidates are now the vertices of the smallest face whose closure contains the target: the three corners of the cell holding it, the two ends of the edge it lies on, or the vertex it coincides with. The face is read off the barycentric coordinates of one containing cell, which is what makes it right. Taking the union of the corners of every containing cell is not enough — a point on an edge lies in two cells and the union offers each cell's opposite apex, which covers one cell and not the other. That was the first attempt at this fix and it did not clear the measured case. Measured inside a ribbon zone of width 0.02: the target sat on an edge whose ends were at exactly 1.000e-02, being the ribbon's own edge vertices, and an apex won by 0.4% at 9.963e-03. Regular vertex rows at half-width make such near-ties the norm rather than bad luck, which is why this showed up in a zone and not on a plain mesh. The reduction now leads with an outside-flag so a rank holding no part of the containing cell cannot win with its own nearer vertex; all ranks flagging outside means the target is off the mesh, and the fallback is the previous nearest-anywhere behaviour. Closes #542. The regression test flattens the mesh first, because on a well-shaped mesh an edge midpoint's nearest vertex is one of that edge's own ends and the defect never appears; a negative control asserts the mesh really does present the case. Against the old rule the invariant fails at 38 of 71 midpoints tested. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/line_cut.py | 92 +++++++++++++++++++++++++-- tests/test_0844_line_cut.py | 71 +++++++++++++++++++++ 2 files changed, 156 insertions(+), 7 deletions(-) diff --git a/src/underworld3/utilities/line_cut.py b/src/underworld3/utilities/line_cut.py index ece69a6e5..8e27293b8 100644 --- a/src/underworld3/utilities/line_cut.py +++ b/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 diff --git a/tests/test_0844_line_cut.py b/tests/test_0844_line_cut.py index 5392c26c8..9640e296a 100644 --- a/tests/test_0844_line_cut.py +++ b/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")