From 9434a63c160b15a3d39d21be85248a71741d196b Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 4 Aug 2026 12:34:59 +1000 Subject: [PATCH 1/3] Place an embedded surface instead of cutting for one (2-D, serial) The cut represents a surface by splitting every edge it crosses, and every restriction it carries follows from that: an edge can be split at one point, so two flanks closer than one element compete for the same edge and the cut is refused; the surface can never be finer than the local h; and a triangle the surface enters but does not leave has no split that represents it. place_along_lines does the same job with the opposite move. It asserts the surface's own points as vertices, deletes the mesh vertices in the way, and retriangulates the cavity so the placed segments survive as element edges. A fault tip terminates inside the mesh, the point spacing is a parameter, and two surfaces may run at any separation. Measured on a 1/16 box: the cut accepts two parallel surfaces one element apart and refuses them at half an element; placement carries them to a tenth of an element. How the cavity is filled. How many ends reach the domain wall decides its shape - none an annulus, one a disc, two a disc per flank - and all three are one walk between two chains, ordered by arc length around the surface's own boundary. At a zero-thickness tip the two flanks meet at a point, so the turn through 180 degrees is given a window of that parameter to itself and is interpolated across by angle about the tip: the tip comes out as a fan of one placed vertex against many cavity vertices, with no width floor. Two things measurement forced that the prototype never met. The walk needs a third move, because a cavity ring is not convex: clipping a protruding corner off as an ear, and where even that fails, swallowing the spike of surviving mesh the walk wedged on and re-clearing. Over 100 random traces on a uniform mesh and 100 on a graded one, that is 2 and 8 failures without it against none with, area exact to 2e-16 throughout. And the walk needs a quality floor rather than only an orientation test - where the cavity reaches the wall the ring runs along it, and two wall vertices plus the surface's end on that wall are three collinear points, which a wall differing in the last bit resolves confidently into a cell of area 1e-15 and a zero angle that every positivity test passes. An end reaching the wall slides the boundary VERTEX along the wall onto it rather than moving the surface, so the surface stays where it was asked for, and the slide is refused where the wall turns so the domain is never deformed. Snapping the trace instead put the chain 9 % of h off. A split wall facet inherits the labels the whole one carried, or a boundary condition steps over the hole left behind. reconnect.rebuild_without_vertices becomes rebuild_cavities and takes placed coordinates: it is now the one rebuild that changes the point chart in both directions. It does not extend the star-forest's leaf set, so placement refuses in parallel rather than returning a mesh whose forest is silently wrong. The cut remains the parallel path, and add_conforming_surface is untouched. Underworld development team with AI support from Claude Code --- .../conforming-surfaces-and-fault-zones.md | 79 ++ src/underworld3/utilities/place_surface.py | 1061 +++++++++++++++++ src/underworld3/utilities/reconnect.py | 113 +- .../parallel/ptest_0844_reconnect_parallel.py | 4 +- tests/test_0848_place_surface.py | 382 ++++++ 5 files changed, 1596 insertions(+), 43 deletions(-) create mode 100644 src/underworld3/utilities/place_surface.py create mode 100644 tests/test_0848_place_surface.py diff --git a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md index ecd77d9b..5a94fffa 100644 --- a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md +++ b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md @@ -210,6 +210,82 @@ Every one of these is raised **collectively**. A rank-local refusal aborts one rank while its peers walk on into the next collective and block there, turning a clear error into a hang. +## The other way: place the surface, do not cut for it + +Everything above describes **cutting** — splitting the edges the surface crosses. +There is a second implementation, +`underworld3.utilities.place_surface.place_along_lines`, which reaches the same +end by the opposite move: it asserts the surface's **own points** as mesh +vertices, deletes the mesh vertices in the way, and retriangulates the cavity so +that the placed segments survive as element edges. + +Both produce a chain of labelled facets with no straddling cell, and both leave +the base mesh untouched. What differs is which restrictions come with the method. +Every one of the cut's refusals above is a consequence of *splitting*, and none of +them applies to placing: + +| | cut | place | +|---|---|---| +| a surface ending **inside** the mesh | refused — a triangle entered and not left has no split that represents it | the tip is a placed vertex like any other, and the cavity closes round it | +| two surfaces **closer than one element** | refused — they cross the same edge, and an edge can be split once | free; neither is competing for an edge of the original mesh | +| surface **finer than the local `h`** | impossible — the surface's vertices *are* the mesh's crossings | `spacing` is a parameter | +| what it does to the mesh | moves vertices onto the surface (`snap_frac`) and splits edges | deletes vertices near the surface and refills the hole | +| parallel | yes, and partition-independent | **serial only** — see below | + +Measured on the same box, viscosity contrast irrelevant: the cut accepts two +parallel surfaces one element apart and refuses them at half an element, while +placement carries them down to a tenth of an element with a worst angle of 4.2 +degrees. + +The construction is the same operation in 2-D on a curve and in 3-D on a sheet — +place, delete, refill — which matters because cutting tetrahedra along a surface +is an unsolved pattern problem while filling a cavity is standard meshing +practice. Only the 2-D half exists. + +### How the cavity is filled + +Clearing the vertices in the way leaves one hole, and what has to be triangulated +inside it depends only on how many ends of the surface reach the domain boundary: +none leaves an annulus, one leaves a disc, two leave a disc per flank. All three +are the same **walk** between two chains — the cavity ring outside, the surface +traversed *out and back* inside — advancing whichever chain is further behind in a +shared parameter. + +The parameter is arc length around the surface's own boundary: down one flank, +around the tip, back up the other. Two cheaper choices were tried and both fail. +Arc length along each chain measures the cavity ring's *wiggle* rather than its +progress, and the two rings drifted four surface widths apart. Position along +strike is discontinuous at the tips, which is exactly where the difficulty lives. + +At a zero-thickness tip the two flanks meet at a point, so the parameter would go +flat over the whole turn through 180 degrees and the walk would have nothing left +to order by. The turn is therefore given a window of the parameter to itself, one +point spacing wide, interpolated by **angle about the tip** — which turns the tip +into a fan of one placed vertex against many cavity vertices. + +The walk has a third move, and without it it wedges. A cavity ring is not convex, +so a corner of it can protrude into the region; clipping such a corner off as an +**ear** consumes both of its ring edges and lets the walk carry on. Where even +that fails the ring vertex is a *spike* of surviving mesh poking into the cavity — +the only triangle that would fill the notch uses a vertex the walk has already +passed — and the answer is to swallow that vertex and re-clear, which is what the +routine does, up to eight times. Measured over 100 random traces on a uniform mesh +and 100 on a graded adapt-on-top mesh: 2 and 8 failures respectively without the +growth step, none with it, and the total area exact to 2e-16 throughout. + +### What it costs + +The walk fills the cavity by parameter, not by shape, so the cells it leaves are +worse than the cut's before repair. Flipping (`reconnect.flip_to_reduce_max_angle`) +fixes that, and both repair passes refuse to touch a labelled edge, so the surface +itself comes through with the same facets. This is the normal composition, not a +workaround. + +Placement is **serial**. It adds points, and a point added across a partition seam +needs the star-forest's leaf set extended — the chart-expansion rebuild, which does +not exist. A parallel call is refused rather than returning a mesh whose +star-forest is silently wrong. The cut remains the parallel path. + ## Limitations * **Two dimensions.** The 3-D mechanism is validated on single tets and small @@ -228,3 +304,6 @@ clear error into a hang. * {doc}`meshing` — mesh construction and `Surface` * {doc}`mesh-metric-redistribution` — the adapt metric that sets the local `h` * `underworld3.utilities.line_cut` — the cutting mechanism +* `underworld3.utilities.place_surface` — the placing mechanism +* `underworld3.utilities.reconnect` — the repair passes, and `rebuild_cavities`, + the rebuild both deletion and placement go through diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py new file mode 100644 index 00000000..d756b938 --- /dev/null +++ b/src/underworld3/utilities/place_surface.py @@ -0,0 +1,1061 @@ +"""Embed a surface in a mesh by PLACING its points, not by cutting for them. + +:mod:`underworld3.utilities.line_cut` embeds a surface by **splitting** every +edge the surface crosses. That works, but every restriction the cut carries is a +consequence of it: an edge can be split at exactly one point, so two flanks +running closer than one element apart compete for the *same* edge and the cut is +refused; the surface can never be finer than the local ``h``; and a triangle the +surface enters but does not leave — a fault tip — has no split that represents +it. + +This module does the same job with the opposite move. It **asserts the +surface's own points**, deletes the mesh vertices in the way, and retriangulates +the cavity so that the placed segments survive as element edges. Nothing is +crossed, so nothing has to be split twice, and none of those restrictions apply: + +* a tip terminates *inside* the mesh, because the tip is a placed vertex like + any other and the cavity closes around it; +* the surface's point spacing is chosen by the caller, not by wherever the mesh + happened to put its edges; +* two surfaces may run arbitrarily close together, because neither is competing + for an edge of the original mesh. + +The operation is the same on a curve in 2-D and on a sheet in 3-D — place, +delete, retriangulate — which is the other reason to prefer it: cutting +tetrahedra along a surface is an unsolved pattern problem, while filling a +cavity is standard meshing practice. Only the 2-D half exists here. + +The three shapes of the cavity +------------------------------ +Deleting the vertices in the way leaves one hole. What has to be triangulated +inside it depends only on how many ends of the surface reach the domain +boundary, and all three cases are the **same walk** between two chains: + +==================== ========================================================= +ends on the boundary the region to fill +==================== ========================================================= +none (a fault) an annulus — the cavity ring outside, the surface + traversed out and back inside. Opened at one rung so that + it, too, is a pair of chains. +one (one tip inside) a disc, its boundary the cavity ring and the out-and-back + traverse, the two meeting at the boundary end. +two (a crossing) two discs, one per flank, each bounded by half the cavity + ring and by the surface. +==================== ========================================================= + +The surface is traversed **out and back**, and the two passes name the SAME +vertices. That is what makes a placed segment one edge with a cell on each side +rather than two edges with a gap between them, and it is the whole of "zero +thickness": the finite-width ribbon this construction was prototyped on is the +same walk with the two passes held apart. + +Ordering the walk: the parameter +-------------------------------- +The two chains have to be marched in step, which needs one parameter that orders +*both* — a point of the surface and a point of the cavity ring alike. Two +cheaper choices fail, for instructive reasons. Arc length along each chain +measures the cavity ring's wiggle rather than its progress, and the two rings +drift apart. Position along strike is discontinuous at the tips, which is +exactly where the difficulty lives. + +:func:`_traverse_parameter` uses arc length around the surface's own boundary — +down one flank, **around the tip**, back up the other. It has neither problem, +because it goes round the tip for the same reason the surface does. At a +zero-thickness tip the two flanks meet at a point, so the turn through 180 +degrees is given a window of its own (``cap``, one point spacing wide) and the +parameter is interpolated across it by ANGLE about the tip. That window is what +turns the tip into a **fan** — one placed vertex, many cavity vertices — instead +of a plateau on which the walk has nothing left to order by. + +Scope +----- +Two dimensions, serial. Placement ADDS points, and the chart-expansion rebuild a +shared point would need does not exist yet, so a parallel call is refused rather +than silently returning a mesh whose star-forest is wrong. See +:func:`~underworld3.utilities.reconnect.rebuild_cavities`. + +Several surfaces are placed **one at a time**, each against the result of the +last. The already-placed segments carry an edge label, and a labelled edge is an +interface, so a later placement will not delete a vertex out of an earlier one. +""" + +import numpy as np + +import underworld3 as uw +from underworld3.utilities import reconnect +from underworld3.utilities.line_cut import (CUT_LABEL, _coords, + _distance_to_lines, _edge_vertices, + _vertex_h, cell_areas, min_angles) + + +# ------------------------------------------------------------- mesh reads (2-D) + +def _cells_anticlockwise(dm, X): + """(n_cells, 3) local vertex indices of every triangle, wound anticlockwise. + + The winding is not cosmetic: :func:`_cavity_ring` reads the boundary of the + dropped set off the cells' own directed edges, and a clockwise cell would + contribute its edges the wrong way round and break the cancellation. + """ + vS, vE = dm.getDepthStratum(0) + cS, cE = dm.getHeightStratum(0) + cells = np.array([[int(p) - vS for p in dm.getTransitiveClosure(c)[0] + if vS <= p < vE] for c in range(cS, cE)], + dtype=np.int64).reshape(cE - cS, 3) + P = X[cells] + turn = ((P[:, 1, 0] - P[:, 0, 0]) * (P[:, 2, 1] - P[:, 0, 1]) + - (P[:, 1, 1] - P[:, 0, 1]) * (P[:, 2, 0] - P[:, 0, 0])) + cells[turn < 0] = cells[turn < 0][:, [0, 2, 1]] + return cells + + +def _boundary_edges(dm): + """Every facet of the domain boundary, as ``(edge point, va, vb)``. + + A facet with one cell in its support is on the boundary of the domain. Its + vertices may never be deleted — that would change the domain — so they are + the one class of vertex the placement has to work around rather than through. + """ + vS, _vE = dm.getDepthStratum(0) + out = [] + for e in range(*dm.getDepthStratum(1)): + if len(dm.getSupport(e)) == 1: + a, b = (int(v) - vS for v in dm.getCone(e)) + out.append((int(e), a, b)) + return out + + +def _boundary_vertices(dm, n_vertices): + """Mask of the vertices lying on the domain boundary.""" + on = np.zeros(n_vertices, dtype=bool) + for _e, a, b in _boundary_edges(dm): + on[a] = on[b] = True + return on + + +def _interface_vertices(dm, n_vertices): + """Vertices carrying an edge of an already-embedded surface. + + Deleting one would put a gap in a surface placed earlier, so these are + refused as victims. Read through :func:`reconnect._interface_edges`, which + already knows which labels mean *interface* and which are PETSc's or UW3's + own bookkeeping — a distinction that has produced three separate silent + failures when guessed at. + """ + pStart, _pEnd = dm.getChart() + vS, _vE = dm.getDepthStratum(0) + locked = reconnect._interface_edges(dm) + on = np.zeros(n_vertices, dtype=bool) + for e in range(*dm.getDepthStratum(1)): + if locked[e - pStart]: + for v in dm.getCone(e): + on[int(v) - vS] = True + return on + + +# ------------------------------------------------------------------ the polyline + +def _arc_length(pts): + """Segment vectors, their lengths, and the cumulative arc length.""" + seg = pts[1:] - pts[:-1] + seglen = np.linalg.norm(seg, axis=1) + return seg, seglen, np.concatenate([[0.0], np.cumsum(seglen)]) + + +def _point_at(pts, a): + """The point of the polyline at arc length ``a``.""" + _seg, seglen, cum = _arc_length(pts) + k = int(np.clip(np.searchsorted(cum, a, side="right") - 1, 0, len(seglen) - 1)) + u = (a - cum[k]) / seglen[k] + return pts[k] + u * (pts[k + 1] - pts[k]) + + +def _nearest_foot(pts, X): + """Arc length of, and signed distance to, the nearest point of the polyline. + + The sign is taken from the nearest SEGMENT's normal. Where a polyline turns + sharply the two segments meeting at the corner disagree about the sign in + the reflex wedge; a fault trace turns gently, and the cavity keeps its ring + a clear element away from the corner. + """ + seg, seglen, cum = _arc_length(pts) + best = np.full(len(X), np.inf) + a = np.zeros(len(X)) + d = np.zeros(len(X)) + for k in range(len(seg)): + u = np.clip(((X - pts[k]) @ seg[k]) / (seg[k] @ seg[k]), 0.0, 1.0) + foot = pts[k] + u[:, None] * seg[k] + dist = np.linalg.norm(X - foot, axis=1) + closer = dist < best + normal = np.array([-seg[k][1], seg[k][0]]) / seglen[k] + best[closer] = dist[closer] + a[closer] = (cum[k] + u * seglen[k])[closer] + d[closer] = ((X - pts[k]) @ normal)[closer] + return a, d + + +def _resample(pts, spacing): + """The polyline with every segment cut into pieces no longer than ``spacing``. + + The original control points are kept, so the geometry is exact and only the + point DENSITY changes. Spacing is what sets the placed surface's resolution + — under the cut that was decided for you by wherever the mesh happened to + put its edges. + """ + out = [pts[0]] + for A, B in zip(pts[:-1], pts[1:]): + span = float(np.linalg.norm(B - A)) + # A repeated control point contributes no segment, and a segment of zero + # length has no direction: it would give the tip cap a normal of 0/0 and + # the whole walk a parameter of NaN, silently. + if span == 0.0: + continue + n = max(int(np.ceil(span / spacing)), 1) + for i in range(1, n + 1): + out.append(A + (i / n) * (B - A)) + return np.array(out) + + +def _inside_mesh(X, cells, points): + """Which of ``points`` lie inside some cell of the mesh.""" + P = X[cells] + v0, v1 = P[:, 1] - P[:, 0], P[:, 2] - P[:, 0] + det = v0[:, 0] * v1[:, 1] - v0[:, 1] * v1[:, 0] + inside = np.zeros(len(points), dtype=bool) + for i, q in enumerate(points): + w = q - P[:, 0] + s = (w[:, 0] * v1[:, 1] - w[:, 1] * v1[:, 0]) / det + t = (v0[:, 0] * w[:, 1] - v0[:, 1] * w[:, 0]) / det + inside[i] = bool(((s >= 0.0) & (t >= 0.0) & (s + t <= 1.0)).any()) + return inside + + +def _clip_to_domain(dm, X, cells, pts): + """The single run of the polyline that lies inside the mesh. + + A trace is normally specified with its ends outside the domain, so that it + crosses cleanly rather than stopping a hair short of the wall. Placement + needs the ends themselves — they become vertices — so the polyline is cut + where it meets the boundary. + + Crossings are found against the boundary FACETS rather than against an + assumed box, so this works on whatever domain the mesh describes. A trace + that leaves and re-enters is refused: it is two surfaces, and they want two + names. + """ + _seg, seglen, cum = _arc_length(pts) + ends = _boundary_edges(dm) + + breaks = [0.0, cum[-1]] + for k in range(len(seglen)): + A, d = pts[k], pts[k + 1] - pts[k] + for _e, ia, ib in ends: + p, e = X[ia], X[ib] - X[ia] + det = d[0] * (-e[1]) - d[1] * (-e[0]) + if det == 0.0: + continue # parallel: no single crossing point + rhs = p - A + t = (rhs[0] * (-e[1]) - rhs[1] * (-e[0])) / det + u = (d[0] * rhs[1] - d[1] * rhs[0]) / det + if 0.0 <= t <= 1.0 and 0.0 <= u <= 1.0: + breaks.append(cum[k] + t * seglen[k]) + + breaks = np.unique(np.round(breaks, 12)) + mids = np.array([_point_at(pts, 0.5 * (breaks[i] + breaks[i + 1])) + for i in range(len(breaks) - 1)]) + live = _inside_mesh(X, cells, mids) + if not live.any(): + raise ValueError( + "no part of the surface lies inside the mesh: the trace and the " + "mesh do not overlap.") + starts = np.flatnonzero( + np.diff(np.concatenate([[0], live.astype(int), [0]])) == 1) + if len(starts) > 1: + raise ValueError( + f"the surface leaves the domain and re-enters it, in {len(starts)} " + "separate pieces. Each piece is a surface of its own and wants its " + "own name; place them one at a time.") + + live_at = np.flatnonzero(live) + # Clamped back onto the polyline: the break positions are rounded so that a + # crossing found twice on two segments is one break, and a rounded end can + # land a hair PAST the trace — which puts the last control point strictly + # inside the kept run, repeats it, and leaves a zero-length final segment. + lo = max(float(breaks[live_at[0]]), 0.0) + hi = min(float(breaks[live_at[-1] + 1]), float(cum[-1])) + interior = [p for a, p in zip(cum, pts) if lo < a < hi] + return np.array([_point_at(pts, lo)] + interior + [_point_at(pts, hi)]) + + +def _spacing_near(dm, X, cells, pts): + """The mean cell diameter along the surface: the mesh's own resolution there.""" + from underworld3.utilities.edge_split import cell_diameters + + diameter = cell_diameters(dm) + near = _distance_to_lines(X[cells].mean(axis=1), [pts]) < diameter + return float(diameter[near].mean() if near.any() else diameter.mean()) + + +# ------------------------------------------------------------------- the cavity + +def _cavity_ring(cells, drop): + """The boundary of the dropped cell set, as a closed anticlockwise ring. + + Each dropped cell contributes its directed edges; an edge whose reverse also + belongs to a dropped cell is interior to the cavity and cancels. What is left + traverses the hole with the cavity on its left — the orientation the walk + needs, obtained without a single geometric test. + + ``None`` when the survivors do not leave one simple hole: a vertex appearing + twice on the ring, or two disconnected cavities. Either means the caller has + to widen the clearance, not that this can be patched up. + """ + directed = {} + for ci in drop: + v0, v1, v2 = cells[ci] + for a, b in ((v0, v1), (v1, v2), (v2, v0)): + directed[(int(a), int(b))] = int(ci) + ring_edges = [(a, b) for (a, b) in directed if (b, a) not in directed] + if not ring_edges: + return None + step = dict(ring_edges) + if len(step) != len(ring_edges): + return None # a vertex leaves the cavity twice + start = min(step) + ring, cur = [start], step[start] + while cur != start: + if cur not in step or len(ring) > len(step): + return None + ring.append(cur) + cur = step[cur] + return ring if len(ring) == len(step) else None + + +def _cells_meeting(X, cells, pts, candidates): + """Which candidate cells a segment of the polyline enters. + + The clearance test alone is not enough. It deletes the vertices NEAR the + surface, but a cell can be crossed by the surface while all three of its + corners sit further away than the clearance — on an obtuse cell, or where + the surface runs nearly parallel to an edge. Such a cell would survive and + straddle, which is the one thing the construction exists to prevent, so the + cells the surface passes through are added outright. + """ + P = X[cells[candidates]] + hit = np.zeros(len(candidates), dtype=bool) + + def side(a, b, c): + return ((b[..., 0] - a[..., 0]) * (c[..., 1] - a[..., 1]) + - (b[..., 1] - a[..., 1]) * (c[..., 0] - a[..., 0])) + + for A, B in zip(pts[:-1], pts[1:]): + for i in range(3): + p, q = P[:, i], P[:, (i + 1) % 3] + # Two segments cross when each separates the other's ends. Touching + # counts, which over-selects by at most a ring of cells and is the + # conservative direction. + hit |= ((side(A, B, p) * side(A, B, q) <= 0.0) + & (side(p, q, A) * side(p, q, B) <= 0.0)) + return np.asarray(candidates)[hit] + + +# ---------------------------------------------------------------- the parameter + +def _traverse_parameter(pts, X, cap, cap_near, cap_far): + """Arc length around the surface's boundary, for points off the surface. + + The surface's boundary is its two flanks and, where an end terminates inside + the mesh, a cap joining them. The parameter runs from 0 along the flank on + the NEGATIVE side, through the far cap, back along the positive side, and + through the near cap. Its two properties are the ones the walk needs and + nothing cheaper has: it is continuous across a tip, and it orders a cavity + vertex and a surface vertex on the same scale. + + ``cap`` is the arc length allotted to a tip. Any positive value orders the + walk correctly; one point spacing gives the tip the same share of the + parameter as an ordinary segment, so the fan around it comes out with about + as many triangles as a segment gets. + + Returns the parameter and the raw arc length of each foot, the latter + because a caller that knows a point's flank from the TOPOLOGY has a better + source for the side than the sign of a distance. + """ + seg, seglen, cum = _arc_length(pts) + length = cum[-1] + + a, d = _nearest_foot(pts, X) + s = np.where(d < 0.0, a, 2.0 * length + cap * cap_far - a) + + # Beyond a tip every point has the SAME foot — the tip — so arc length says + # nothing about where it sits around the turn. The angle about the tip does, + # and it agrees with the flank formula at both ends of the window. + if cap_far: + t = seg[-1] / seglen[-1] + n = np.array([-t[1], t[0]]) + r = X - pts[-1] + phi = np.arctan2(r @ n, r @ t) + s = np.where(r @ t > 0.0, length + cap * (phi + 0.5 * np.pi) / np.pi, s) + if cap_near: + t = seg[0] / seglen[0] + n = np.array([-t[1], t[0]]) + r = X - pts[0] + psi = np.arctan2(r @ n, -(r @ t)) + s = np.where(r @ t < 0.0, + 2.0 * length + cap * cap_far + + cap * (0.5 * np.pi - psi) / np.pi, s) + return s, a + + +def _traverse(index, pts, cap, cap_near, cap_far): + """The surface out along one flank and back, as ``(vertex, parameter)``. + + ``index`` gives the vertex carrying each point of the polyline. A point in + the middle of the surface appears TWICE — once per flank — as the same + vertex; a tip inside the mesh appears once, and an end on the domain + boundary opens the chain there. + """ + _seg, _seglen, cum = _arc_length(pts) + length = cum[-1] + last = len(pts) - 1 + shift = 2.0 * length + cap * cap_far + + out = [] if cap_near else [(int(index[0]), 0.0)] + out += [(int(index[k]), cum[k]) for k in range(1, last)] + out.append((int(index[last]), length + 0.5 * cap if cap_far else length)) + out += [(int(index[k]), shift - cum[k]) for k in range(last - 1, 0, -1)] + out.append((int(index[0]), shift + 0.5 * cap if cap_near else shift)) + return out + + +# --------------------------------------------------------------------- the walk + +#: A triangle the walk may close has to have some shape to it, not merely a +#: resolvable sign. The case that forces this is a cavity reaching the domain +#: wall: the ring then runs ALONG the wall, and any triangle made of two wall +#: vertices and the surface's end on that same wall is three collinear points. +#: The orientation predicate declines an exactly-collinear triple, but a wall +#: whose vertices differ in the last bit gives a resolvable sign and an area of +#: 1e-15 — which passes every positivity test and is a cell with a zero angle. +#: Refusing it makes the walk advance along the SURFACE instead, which is the +#: triangulation that was wanted. Well below anything a real cell reaches: a +#: cell of a fifth of a degree measures 9e-3 on this scale. +_MIN_QUALITY = 1.0e-6 + + +def _quality(P): + """Scale-free triangle quality: 1 equilateral, 0 degenerate, <0 inverted.""" + e = np.array([P[2] - P[1], P[0] - P[2], P[1] - P[0]]) + twice_area = ((P[1, 0] - P[0, 0]) * (P[2, 1] - P[0, 1]) + - (P[1, 1] - P[0, 1]) * (P[2, 0] - P[0, 0])) + return float(2.0 * np.sqrt(3.0) * twice_area + / max(float((e ** 2).sum()), np.finfo(float).tiny)) + + +def _usable(tri, X): + """Whether the walk may close this triangle: right way round, and not flat.""" + return (reconnect._orient2d(X[tri[0]], X[tri[1]], X[tri[2]]) > 0 + and _quality(X[list(tri)]) > _MIN_QUALITY) + + +def _encloses(tri, X, points): + """Whether any of ``points`` lies strictly inside the triangle.""" + A, B, C = X[tri[0]], X[tri[1]], X[tri[2]] + v0, v1 = B - A, C - A + det = v0[0] * v1[1] - v0[1] * v1[0] + for p in points: + if p in tri: + continue + w = X[p] - A + s = (w[0] * v1[1] - w[1] * v1[0]) / det + t = (v0[0] * w[1] - v0[1] * w[0]) / det + if s > 0.0 and t > 0.0 and s + t < 1.0: + return True + return False + + +def _zip(outer, inner, s_out, s_in, X): + """Triangulate the region between two chains, advancing one side at a time. + + Both chains run in increasing parameter and bound the same region, which + lies on the left of each. Every step closes one triangle on the chain that + is further behind, so the two are consumed in step; the other side is tried + when the preferred one would invert. The chains may SHARE an end vertex — + where the surface meets the domain boundary — and the triangle that would + close there is degenerate, so it is skipped rather than emitted. + + Choosing the side by parameter rather than by shape is not a detail. With + the choice made on shape alone one chain runs away from the other: measured + on the ribbon prototype, the inner advanced 23 points against the outer's 10 + before the correspondence was 14 surface widths out and every triangle + inverted. + + There is a third move, and without it the walk wedges more often. A cavity + ring is not convex — it follows whatever cells were cleared — so a corner of + it can protrude into the region, and neither cross-triangle at such a corner + turns the right way. Clipping that corner off the ring as an EAR consumes + both of its ring edges and lets the walk carry on. It is a last resort + rather than a preference, because it advances the ring without advancing the + surface, and doing that by choice would let the two run out of step. The ear + is refused if it would swallow another vertex of either chain, which would + overlap the triangles that vertex still has to be given. + + Only the ring may be clipped this way. An ear on the surface's own chain + would drop one of its points, which is the whole thing being placed. + + Returns ``(triangles, None)``, or ``(None, vertex)`` naming the ring vertex + the walk could not get past. That vertex is a SPIKE of surviving mesh poking + into the cavity: the cavity wraps more than half way round it, so the only + triangle that could fill the notch uses a vertex the walk has already + consumed, and no forward walk can reach it. The caller's answer is to + sacrifice the vertex — it is one the cavity should have swallowed in the + first place — rather than to make the walk able to go backwards. + """ + outer, s_out = list(outer), list(s_out) + inner, s_in = list(inner), list(s_in) + tris = [] + i = j = 0 + while i < len(outer) - 1 or j < len(inner) - 1: + m, n = len(outer) - 1, len(inner) - 1 + prefer_outer = i < m and (j >= n or s_out[i + 1] <= s_in[j + 1]) + taken = None + for side in (("o", "i") if prefer_outer else ("i", "o")): + if side == "o" and i < m: + tri = (outer[i], outer[i + 1], inner[j]) + elif side == "i" and j < n: + tri = (outer[i], inner[j + 1], inner[j]) + else: + continue + if len(set(tri)) < 3: + taken = (side, None) # the chains meet: nothing to close + break + if _usable(tri, X): + taken = (side, tri) + break + + if taken is None: + ear = (outer[i], outer[i + 1], outer[i + 2]) if i + 2 <= m else None + if ear is None or not _usable(ear, X) or _encloses( + ear, X, outer[i:] + inner[j:]): + return None, outer[min(i + 1, m)] + tris.append(ear) + del outer[i + 1], s_out[i + 1] + continue + + side, tri = taken + if tri is not None: + tris.append(tri) + if side == "o": + i += 1 + else: + j += 1 + return tris, None + + +def _opened(chain, s): + """A closed chain cut open at its lowest parameter, ready for :func:`_zip`. + + The cut is a bridge the walk crosses twice, once at each end, so the edge it + introduces is an ordinary interior edge with a cell on either side and not a + seam in the result. + """ + start = int(np.argmin(s)) + rolled = list(chain[start:]) + list(chain[:start]) + t = np.roll(np.asarray(s, dtype=float), -start) - float(np.min(s)) + return rolled + [rolled[0]], np.concatenate([t, [1.0]]) + + +def _walk_cavity(ring, s_ring, traverse, s_in, X, corners): + """Fill the cavity, in the one or two pieces the surface leaves of it. + + ``corners`` are the ring vertices where the surface meets the domain + boundary, in polyline order: none for a fault, one for a fault with one tip + inside, two for a surface crossing the domain. Returns what :func:`_zip` + returns — the triangles, or the ring vertex the walk could not get past. + """ + if not corners: + outer, s_out = _opened(ring, s_ring) + inner, t_in = _opened(traverse, s_in) + return _zip(outer, inner, s_out, t_in, X) + + start = ring.index(corners[0]) + ring = ring[start:] + ring[:start] + s_ring = np.concatenate([s_ring[start:], s_ring[:start]]) + # A corner sits ON the surface, where the flank test that gives every other + # ring vertex its parameter has no side to report. Its value is known + # exactly from the traverse instead. + s_ring[0] = s_in[0] + + if len(corners) == 1: + return _zip(ring + [ring[0]], traverse, + np.concatenate([s_ring, [s_in[-1]]]), s_in, X) + + q = ring.index(corners[1]) + j = traverse.index(corners[1]) + s_ring[q] = s_in[j] + out = [] + for chain, s_chain, sub, t_sub in ( + (ring[:q + 1], s_ring[:q + 1], traverse[:j + 1], s_in[:j + 1]), + (ring[q:] + [ring[0]], np.concatenate([s_ring[q:], [s_in[-1]]]), + traverse[j:], s_in[j:])): + piece, stuck = _zip(chain, sub, s_chain, t_sub, X) + if piece is None: + return None, stuck + out += piece + return out, None + + +# ------------------------------------------------------------- labels and edges + +def _edge_lookup(dm): + """``{(v0, v1): edge point}`` over the whole mesh, keyed low vertex first.""" + out = {} + for e in range(*dm.getDepthStratum(1)): + a, b = (int(v) for v in dm.getCone(e)) + out[(a, b) if a < b else (b, a)] = int(e) + return out + + +def _label_placed_edges(dm, chain, name, value): + """Mark the edges the surface's segments became; return how many. + + Read off the VERTEX NUMBERS the rebuild handed back rather than by geometry: + the placed points are known exactly, so a geometric test could only weaken + an identity that is already exact. A segment that is not an edge of the + result is a defect, and shows up here immediately. + """ + edge_of = _edge_lookup(dm) + if not dm.hasLabel(name): + dm.createLabel(name) + label = dm.getLabel(name) + label.setDefaultValue(0) + for a, b in zip(chain[:-1], chain[1:]): + e = edge_of.get((a, b) if a < b else (b, a)) + if e is None: + raise RuntimeError( + "a segment of the placed surface is not an edge of the result: " + "the cavity was triangulated across the surface rather than up " + "to it.") + label.setValue(e, int(value)) + return len(chain) - 1 + + +def _inherit_boundary_labels(new_dm, dm, splits): + """Give the halves of a split boundary facet the labels the whole one had. + + Placing a surface's end on the domain boundary replaces one boundary facet + by two. The rebuild carries labels by point id, and the facet that was split + has no point id in the result, so its ``Left`` / ``UW_Boundaries`` / named + boundary values would simply be dropped — leaving a hole in the wall that + every boundary condition applied there would quietly step over. + + ``splits`` is ``(source edge point, end vertex, placed vertex, end vertex)`` + in the RESULT's numbering for the three vertices. + """ + edge_of = _edge_lookup(new_dm) + for source, a, v, b in splits: + for i in range(dm.getNumLabels()): + name = dm.getLabelName(i) + if name in reconnect._TOPOLOGY_LABELS: + continue + value = dm.getLabel(name).getValue(source) + if value < 0: + continue # this label says nothing about it + if not new_dm.hasLabel(name): + new_dm.createLabel(name) + target = new_dm.getLabel(name) + for x, y in ((a, v), (v, b)): + target.setValue(edge_of[(x, y) if x < y else (y, x)], int(value)) + + +# -------------------------------------------------------------------- placement + +def _place_one(dm, pts, label, label_value, clearance, spacing, end_snap): + """Place one polyline into ``dm``. Returns the new mesh and its counts.""" + X = _coords(dm) + cells = _cells_anticlockwise(dm, X) + wall = [X[[a, b]] for _e, a, b in _boundary_edges(dm)] + + pts = _clip_to_domain(dm, X, cells, pts) + if spacing is None: + spacing = _spacing_near(dm, X, cells, pts) + pts = _resample(pts, spacing) + + # An end that reaches the domain boundary becomes a vertex ON it; an end + # inside the mesh is a tip and the cavity closes round it. Where exactly one + # end is on the boundary the polyline is turned so that it is the FIRST, + # which is what lets the traverse open there and stay a single chain. + at_wall = list(_distance_to_lines(pts[[0, -1]], wall) < 1.0e-9 * spacing) + if at_wall[1] and not at_wall[0]: + pts, at_wall = np.ascontiguousarray(pts[::-1]), at_wall[::-1] + + dm, X, corner, split_at = _settle_ends(dm, X, pts, spacing, end_snap, at_wall) + vS, _vE = dm.getDepthStratum(0) + cS, _cE = dm.getHeightStratum(0) + cells = _cells_anticlockwise(dm, X) + on_boundary = _boundary_vertices(dm, len(X)) + + # Which vertices are in the way. A domain-boundary vertex is never one: + # deleting it would change the domain. Nor is a vertex of a surface already + # embedded, which would put a gap in that surface. + protected = on_boundary | _interface_vertices(dm, len(X)) + victim = ((_distance_to_lines(X, [pts]) + < clearance * _vertex_h(X, _edge_vertices(dm))) & ~protected) + + # One local index space covering the surviving mesh vertices and the placed + # points, so that the two can appear in one triangle without a special case. + placed_rows, index = [], np.empty(len(pts), dtype=np.int64) + for k in range(len(pts)): + if corner.get(k) is not None: + index[k] = corner[k] + else: + index[k] = len(X) + len(placed_rows) + placed_rows.append(pts[k]) + placed = np.array(placed_rows).reshape(-1, 2) + Xall = np.vstack([X, placed]) + + cap = spacing + tip_near, tip_far = 0 not in corner, (len(pts) - 1) not in corner + length = _arc_length(pts)[2][-1] + total = 2.0 * length + cap * (tip_near + tip_far) + traverse = _traverse(index, pts, cap, tip_near, tip_far) + corners = [int(index[k]) for k in sorted(corner)] + + edge = np.linalg.norm(X[cells[:, 0]] - X[cells[:, 1]], axis=1) + reachable = np.flatnonzero( + _distance_to_lines(X[cells].mean(axis=1), [pts]) < edge + spacing) + crossed = _cells_meeting(X, cells, pts, reachable) + + # The cavity is CLEARED and filled in one loop, because whether it is big + # enough is not knowable in advance: the walk can only fail at a spike of + # surviving mesh poking into it, and the answer to a spike is to swallow the + # vertex at its point. Each round strictly grows the victim set, so this + # terminates; the cap is against a spike the mesh will not give up — one on + # the domain wall, or on a surface already embedded. + for _attempt in range(8): + drop = np.union1d(np.flatnonzero(victim[cells].any(axis=1)), crossed) + if not len(drop): + raise ValueError( + "the surface meets no cell of this mesh: there is nothing to " + "place it in.") + ring = _cavity_ring(cells, drop) + if ring is None: + raise RuntimeError( + "the cells cleared for the surface do not leave one simple " + "hole. Raise `clearance` so the cavity is wider than the shapes " + "pinching it.") + if victim[ring].any(): + raise RuntimeError( + "a deleted vertex is on the cavity boundary; the cavity is not " + "the union of the victims' stars.") + + # An end placed part-way along a boundary facet splits it, so it belongs + # on the cavity ring between that facet's two vertices. + for k, (_e, a, b) in split_at.items(): + ring = _insert_on_ring(ring, a, b, int(index[k])) + + s_ring, foot = _traverse_parameter(pts, Xall[ring], cap, tip_near, + tip_far) + # With both ends on the boundary the ring is cut in two and each piece + # belongs to a known flank, so the side comes from the TOPOLOGY rather + # than from the sign of a distance — which a ring vertex sitting a hair + # on the wrong side of the trace would otherwise get wrong, and with it + # its whole position in the walk. + if len(corners) == 2: + s_ring = _flank_parameter(ring, foot, length, corners) + + tris, stuck = _walk_cavity(list(ring), s_ring / total, + [v for v, _s in traverse], + np.array([s for _v, s in traverse]) / total, + Xall, corners) + if tris is not None: + break + if stuck is None or stuck >= len(X) or protected[stuck]: + raise RuntimeError( + "the cavity could not be triangulated around the surface, and " + "the vertex it wedged on may not be deleted — it is on the " + "domain boundary, on a surface already embedded, or a point of " + "this one. Move the surface off it, or raise `clearance`.") + victim[stuck] = True + else: + raise RuntimeError( + "the cavity could not be triangulated around the surface after " + "growing it 8 times. Raise `clearance`, or place the surface's " + "points further apart.") + + def mixed(v): + return int(v) + vS if v < len(X) else -(int(v) - len(X) + 1) + + new_dm, point_map, placed_points = reconnect.rebuild_cavities( + dm, np.flatnonzero(victim) + vS, drop + cS, + [tuple(mixed(v) for v in t) for t in tris], placed) + + def new_point(v): + return (int(placed_points[v - len(X)]) if v >= len(X) + else int(point_map[v + vS])) + + n_facets = _label_placed_edges(new_dm, [new_point(int(v)) for v in index], + label, label_value) + _inherit_boundary_labels( + new_dm, dm, [(e, new_point(a), new_point(int(index[k])), new_point(b)) + for k, (e, a, b) in split_at.items()]) + return new_dm, {"n_placed": len(placed), + "n_on_surface": len(pts) - len(placed), + "n_removed": int(victim.sum()), + "n_surface_facets": n_facets} + + +def _wall_facet_at(dm, X, point, spacing): + """The boundary facet a placed end lands on: ``(edge point, va, vb, u)``. + + ``u`` is the position along the facet, so a caller can tell an end landing + in the middle of a facet from one landing all but on top of a vertex. + """ + best, hit = np.inf, None + for e, a, b in _boundary_edges(dm): + d = X[b] - X[a] + u = float(np.clip(((point - X[a]) @ d) / (d @ d), 0.0, 1.0)) + away = float(np.linalg.norm(point - (X[a] + u * d))) + if away < best: + best, hit = away, (e, a, b, u) + if best > 1.0e-9 * spacing: + raise RuntimeError( + "the surface's end was clipped to the domain boundary but lies on " + "no boundary facet; the trace and the mesh geometry disagree.") + return hit + + +def _wall_is_straight(dm, X, v, spacing): + """Whether the domain boundary runs straight through vertex ``v``. + + A vertex on a straight run of wall may be slid ALONG the wall without + changing the domain at all. A vertex where the wall turns may not: moving it + would move the corner, so it is left where it is and the surface's end is + placed beside it instead. + """ + facets = [(a, b) for _e, a, b in _boundary_edges(dm) if v in (a, b)] + if len(facets) != 2: + return False + directions = [] + for a, b in facets: + other = b if a == v else a + d = X[other] - X[v] + directions.append(d / np.linalg.norm(d)) + cross = abs(directions[0][0] * directions[1][1] + - directions[0][1] * directions[1][0]) + return bool(cross < 1.0e-9) + + +def _settle_ends(dm, X, pts, spacing, end_snap, at_wall): + """Give every end that reaches the wall a mesh vertex, MOVING THE MESH. + + An end on the domain boundary has to be a vertex of the result. Placing one + a hair from an existing boundary vertex would carve a sliver against the + wall that no later pass may repair, since a boundary vertex can never be + deleted — so an end landing within ``end_snap`` of a facet's end slides that + vertex along the wall onto it instead. + + Moving the mesh rather than the surface is the same choice + :func:`~underworld3.utilities.line_cut.pull_vertex_onto` makes, for the same + reason: the surface's position is a design variable in an outer + optimisation, and it has to end up where it was asked for. The cost is mesh + displacement of at most half a boundary facet, and it is paid only where the + wall is locally straight, so the domain itself is never deformed. + + Returns the working mesh (a copy, when a vertex moved), its coordinates, the + vertex carrying each end, and the facets an end has to be inserted into. + """ + from underworld3.utilities.line_cut import _set_coordinates + + corner, split_at, moves = {}, {}, {} + for k, here in ((0, at_wall[0]), (len(pts) - 1, at_wall[1])): + if not here: + continue # a tip: the cavity closes round it + e, a, b, u = _wall_facet_at(dm, X, pts[k], spacing) + # Put the end exactly ON the facet. Clipping solves a 2x2 system, so it + # lands within rounding of the wall rather than on it, and a vertex a + # picometre outside the domain is still outside it. + pts[k] = X[a] + u * (X[b] - X[a]) + near = a if u < 0.5 else b + if np.linalg.norm(X[near] - pts[k]) < 1.0e-9 * spacing: + corner[k] = near # already exactly there + elif min(u, 1.0 - u) < end_snap and _wall_is_straight(dm, X, near, spacing): + corner[k] = near + moves[near] = pts[k] + else: + corner[k] = None + split_at[k] = (e, a, b) + + if not moves: + return dm, X, corner, split_at + work = dm.clone() + X = X.copy() + for v, target in moves.items(): + X[v] = target + _set_coordinates(work, np.array(sorted(moves), dtype=np.int64), + X[sorted(moves)]) + return work, X, corner, split_at + + +def _insert_on_ring(ring, a, b, v): + """Put ``v`` between the consecutive ring vertices ``a`` and ``b``.""" + for i, pair in enumerate(zip(ring, ring[1:] + ring[:1])): + if set(pair) == {a, b}: + return ring[:i + 1] + [v] + ring[i + 1:] + raise RuntimeError( + "the boundary facet the surface ends on is not on the cavity ring, so " + "the cell holding it was not cleared. Raise `clearance`.") + + +def _flank_parameter(ring, foot, length, corners): + """Re-derive the ring's parameter from which flank each piece of it is on. + + Only for a surface crossing the domain, where the ring genuinely divides: + the two corners cut it in two, so every vertex of a piece is on one flank + whatever the sign of its distance to the trace happens to say. Reading the + side off the topology instead is what stops a ring vertex sitting a hair on + the wrong side of the trace from being given the far flank's parameter, and + with it the wrong place in the walk entirely. + """ + first = np.zeros(len(ring), dtype=bool) + i, stop = ring.index(corners[0]), ring.index(corners[1]) + while True: + first[i] = True + if i == stop: + break + i = (i + 1) % len(ring) + return np.where(first, foot, 2.0 * length - foot) + + +def place_along_lines(dm, lines, label=CUT_LABEL, label_value=1, + clearance=0.55, spacing=None, end_snap=0.25, + verbose=False): + """Embed surfaces in a mesh by placing their points and refilling the hole. + + The surface's own points become mesh vertices, the mesh vertices in the way + are deleted, and the cavity is retriangulated so that every segment of the + surface is an element edge carrying ``label``. Nothing is split, so — unlike + :func:`~underworld3.utilities.line_cut.cut_along_lines` — a surface may end + INSIDE the mesh, may be finer than the local ``h``, and may run alongside + another surface at any separation. + + Parameters + ---------- + dm : PETSc.DMPlex + A 2-D simplex mesh. **Not modified** — the result is a new mesh, so a + surface can be moved and re-placed against the same fixed base. + lines : sequence of array_like + One or more polylines, each an ``(N, 2)`` array of points. Ends outside + the mesh are clipped to the boundary, so the usual "specify it a little + long" convention works. They are placed **one at a time**, each into the + result of the last, and must not intersect one another. + label, label_value : str, int + Name and stratum value of the label put on the surface's edges. The + default is the one :mod:`~underworld3.utilities.line_cut` uses, so a + downstream pass — ``relax(pin_bands=...)``, the reconnection passes — + does not have to know which routine embedded the surface. + clearance : float + Delete a mesh vertex within this multiple of its own local ``h`` of the + surface. This is what buys room for the placed points: too small and the + cavity is too tight to triangulate around them, too large and more of + the mesh is rebuilt than needs to be. + spacing : float or None + Distance between placed points along the surface. ``None`` takes the + mean cell diameter near the surface, putting the surface's resolution + where the mesh's already is. This is the knob the cut does not have: a + placed surface may be finer than the mesh it is placed in. + end_snap : float + An end reaching the domain boundary within this fraction of a boundary + facet's length from one of its vertices slides that vertex ALONG the + wall onto the end, rather than placing a second vertex a hair from it. A + boundary vertex can never be deleted, so that sliver is one no later + pass could repair. The surface still ends exactly where it was asked to + — the mesh moves, not the surface — and the slide is refused where the + wall turns, so the domain is never deformed. + verbose : bool + Report the counts and the worst cell of the result. + + Returns + ------- + placed : PETSc.DMPlex + A new mesh in which every segment of every surface is an edge carrying + ``label``. + info : dict + ``n_placed`` vertices this created on the surfaces, ``n_on_surface`` + existing vertices it reused, ``n_removed`` vertices it deleted, + ``n_surface_facets`` edges labelled, and ``min_area`` and ``min_angle`` + of the result. + + One surface is one chain, so its facets number one fewer than the + vertices along it — the same identity the cut reports, with the placed + count standing where the split count stood: + + ``n_surface_facets == n_placed + n_on_surface - len(lines)``. + + Raises + ------ + ValueError + If a surface does not overlap the mesh, or leaves the domain and + re-enters it. + RuntimeError + If the cleared cells do not leave one simple hole, or the cavity cannot + be triangulated around the surface without inverting a cell. Both mean + ``clearance`` is too small for the surface asked for. + NotImplementedError + In 3-D, or in parallel. + + Examples + -------- + A fault that stops inside the mesh — the case the cut refuses outright: + + >>> tip = numpy.array([[0.2, 0.5], [0.6, 0.55]]) + >>> placed, info = place_along_lines(mesh.dm, [tip], label="Fault") + >>> info["n_surface_facets"] == info["n_placed"] + info["n_on_surface"] - 1 + True + + See Also + -------- + underworld3.utilities.line_cut.cut_along_lines : the same job by splitting + the edges the surface crosses. + underworld3.utilities.reconnect.remove_vertices : the repair pass that + cleans up what the walk leaves, and which will not touch a labelled edge. + """ + if dm.getDimension() != 2: + raise NotImplementedError( + f"place_along_lines is 2-D; this mesh is {dm.getDimension()}-D. The " + "cavity of a placed sheet is a polyhedron, and filling one can need " + "Steiner points that the 2-D walk never has to invent.") + if uw.mpi.size > 1: + raise NotImplementedError( + "place_along_lines is serial. Placement ADDS points, so a surface " + "crossing a partition seam needs the star-forest's leaf set " + "extended — the chart-expansion rebuild, which does not exist yet. " + "Refusing beats returning a mesh whose star-forest is wrong.") + + out = dm + totals = {"n_placed": 0, "n_on_surface": 0, "n_removed": 0, + "n_surface_facets": 0} + for pts in lines: + out, one = _place_one(out, np.asarray(pts, dtype=float)[:, :2], + label, label_value, clearance, spacing, end_snap) + for key in totals: + totals[key] += one[key] + + areas = cell_areas(out) + over = sum(1 for f in range(*out.getHeightStratum(1)) + if len(out.getSupport(f)) > 2) + if over: + raise RuntimeError( + f"{over} facet(s) of the result have more than two cells: the " + "retriangulated cavity is not conforming.") + if (areas <= 0.0).any(): + raise RuntimeError( + f"{int((areas <= 0.0).sum())} cell(s) of the result are inverted.") + + info = dict(totals, min_area=float(areas.min()), + min_angle=float(min_angles(out).min())) + if verbose: + uw.pprint(f"[place {label!r}] placed {info['n_placed']} vertices, " + f"reused {info['n_on_surface']}, removed {info['n_removed']}; " + f"{info['n_surface_facets']} surface facets, min angle " + f"{info['min_angle']:.2f} deg") + return out, info diff --git a/src/underworld3/utilities/reconnect.py b/src/underworld3/utilities/reconnect.py index d809a47c..6d5d6a4d 100644 --- a/src/underworld3/utilities/reconnect.py +++ b/src/underworld3/utilities/reconnect.py @@ -152,7 +152,7 @@ Deletion freezes the same seam, for a stronger reason: it compacts the point chart, so unlike a flip it cannot hand the star-forest across verbatim. Every point after a deleted one shifts, and each leaf's *remote* index is a number only -its owner holds. :func:`rebuild_without_vertices` renumbers locally and +its owner holds. :func:`rebuild_cavities` renumbers locally and broadcasts the new numbering root-to-leaf **once** to close that gap. One exchange of bookkeeping over the existing partition — no cell changes rank and nothing is redistributed, which is the property an external remesher costs us. @@ -446,16 +446,16 @@ def _cell_vertices_and_seam(dm, X, shared): # ------------------------------------------------------------------- the rebuild -def _write_coordinates(new, dm, vertex_range, source): - """Give ``new`` a local vertex coordinate section holding ``source``'s rows. +def _write_coordinates(new, cdim, vertex_range, values): + """Give ``new`` a local vertex coordinate section holding ``values``. - ``vertex_range`` is the new mesh's vertex stratum and ``source`` indexes the - source mesh's coordinate rows, one per new vertex — the identity for a - rebuild that preserves the numbering, and the survivor list for one that - compacts the chart. Written through a section rather than ``setCoordinates`` - because this is purely local data and the latter wants a global vector. + ``vertex_range`` is the new mesh's vertex stratum and ``values`` is one row + per vertex of it, in point order — the source's own rows for a rebuild that + preserves the numbering, the survivors' rows for one that compacts the + chart, and the survivors followed by the newly placed points for one that + grows it. Written through a section rather than ``setCoordinates`` because + this is purely local data and the latter wants a global vector. """ - cdim = dm.getCoordinateDim() vS, vE = vertex_range new.setCoordinateDim(cdim) section = new.getCoordinateSection() @@ -468,8 +468,7 @@ def _write_coordinates(new, dm, vertex_range, source): section.setUp() coords = PETSc.Vec().createSeq(section.getStorageSize(), comm=PETSc.COMM_SELF) - X = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, cdim) - coords.array[:] = X[source].reshape(-1) + coords.array[:] = np.asarray(values, dtype=float).reshape(-1) new.setCoordinatesLocal(coords) @@ -576,7 +575,7 @@ def rebuild_with_cones(dm, new_cells, new_edges): # Coordinates verbatim: the vertex points are unchanged, so this is the same # section over the same chart holding the same values. - _write_coordinates(new, dm, (vS, vE), np.arange(vE - vS)) + _write_coordinates(new, dm.getCoordinateDim(), (vS, vE), _coords(dm)[:vE - vS]) _copy_labels(new, dm) # The star-forest transfers verbatim: every rank preserves its numbering, so @@ -586,8 +585,13 @@ def rebuild_with_cones(dm, new_cells, new_edges): return new -def rebuild_without_vertices(dm, victims, drop_cells, new_cells): - """Build a fresh plex with vertices deleted and their links retriangulated. +def rebuild_cavities(dm, victims, drop_cells, new_cells, placed_coords=()): + """Build a fresh plex with cells replaced by a retriangulation of their cavity. + + The one rebuild that changes the point chart, in both directions: vertices + may be **deleted** (a repair pass dissolving a bad point) and vertices may be + **placed** (an embedded surface asserting its own points). A caller doing + only one of the two passes an empty list for the other. Parameters ---------- @@ -596,32 +600,39 @@ def rebuild_without_vertices(dm, victims, drop_cells, new_cells): victims : sequence of int Vertex points to delete. drop_cells : sequence of int - Cell points to delete — the union of the victims' stars. + Cell points to delete — every cell of the cavity being retriangulated. new_cells : sequence of tuple - Replacement cells as anticlockwise vertex triples, in the **source** - point numbering. + Replacement cells as vertex triples. A **non-negative** entry is a point + number in the source mesh; a **negative** entry ``-(k + 1)`` is row ``k`` + of ``placed_coords``, a vertex this call creates. Winding is fixed here, + so the caller need not order the triples. + placed_coords : array_like + ``(n, cdim)`` coordinates of the vertices being added, or empty. Returns ------- new : PETSc.DMPlex - The rebuilt mesh, on a compacted chart. + The rebuilt mesh, on a re-numbered chart. point_map : numpy.ndarray Chart-indexed source point -> new point, ``-1`` for a deleted point. + placed_points : numpy.ndarray + New point number of each row of ``placed_coords``. Notes ----- This is :func:`rebuild_with_cones` with its one restriction lifted. A flip adds and removes no points, so that function can preserve the numbering and - hand the star-forest across verbatim. A deletion cannot: the chart shrinks, - and every point after a deleted one shifts. So the numbering is rebuilt, and - with it the edges — which are not given, but **derived from the new cells**, - since an edge of the retriangulated cavity may be an old edge that survived - or a chord the ear-clip invented and there is no way to tell them apart - except by looking. - - Ordering is by source point number for everything that survives and by - vertex tuple for everything new, so the result is a function of the input - topology and not of the order the caller happened to accumulate it in. + hand the star-forest across verbatim. A change of chart cannot: every point + after a deleted one shifts, and a placed one has no source number at all. So + the numbering is rebuilt, and with it the edges — which are not given, but + **derived from the new cells**, since an edge of the retriangulated cavity + may be an old edge that survived or a chord the retriangulation invented and + there is no way to tell them apart except by looking. + + Ordering is by source point number for everything that survives, then by + caller order for the placed vertices, and by vertex tuple for the new cells, + so the result is a function of the input topology and not of the order the + caller happened to accumulate it in. The parallel cost is one exchange. Renumbering the star-forest's *local* indices is local knowledge, but the *remote* index of each leaf is the @@ -630,6 +641,15 @@ def rebuild_without_vertices(dm, victims, drop_cells, new_cells): is bookkeeping over the existing partition; no cell moves rank, and nothing is redistributed. + .. warning:: + + A **placed** vertex has no remote number anywhere, so the one-exchange + renumbering does not cover it: placing points across a partition seam + needs the leaf set itself extended, which is the chart-EXPANSION rebuild + and is not this function. Deleting is safe in parallel because + :func:`remove_vertices` never offers a shared point; placing is currently + refused before it gets here. + The caller is responsible for never deleting a point that the star-forest touches (see :func:`remove_vertices`), which is what lets the leaf set carry across unchanged rather than having to be recomputed. @@ -638,6 +658,7 @@ def rebuild_without_vertices(dm, victims, drop_cells, new_cells): cS, cE = dm.getHeightStratum(0) vS, vE = dm.getDepthStratum(0) eS, eE = dm.getDepthStratum(1) + cdim = dm.getCoordinateDim() dead_v = np.zeros(vE - vS, dtype=bool) dead_v[np.asarray(victims, dtype=np.int64) - vS] = True @@ -646,6 +667,7 @@ def rebuild_without_vertices(dm, victims, drop_cells, new_cells): surv_v = np.flatnonzero(~dead_v) + vS surv_c = np.flatnonzero(~dead_c) + cS + placed = np.asarray(placed_coords, dtype=float).reshape(-1, cdim) # Cells, in the source numbering, as vertex triples: survivors keep their # relative order, the replacements follow sorted by vertex tuple. Every @@ -654,9 +676,13 @@ def rebuild_without_vertices(dm, victims, drop_cells, new_cells): # negative volume without raising. X = _coords(dm) + def at(v): + """Coordinates of a mixed index: a source point, or a placed row.""" + return placed[-int(v) - 1] if v < 0 else X[int(v) - vS] + def anticlockwise(tri): a, b, c = (int(v) for v in tri) - if _orient2d(X[a - vS], X[b - vS], X[c - vS]) < 0: + if _orient2d(at(a), at(b), at(c)) < 0: return (a, c, b) return (a, b, c) @@ -687,25 +713,30 @@ def anticlockwise(tri): - {pair_of[e] for e in surv_e}) - # Strata keep the source's relative order; only their sizes change. - sizes = {"c": len(cells), "v": len(surv_v), "e": len(edges)} - offset, at = {}, pStart + # Strata keep the source's relative order; only their sizes change. The + # placed vertices sit at the END of the vertex stratum, after the survivors, + # so a source vertex's new number depends only on how many vertices before + # it died and not on how many were placed. + sizes = {"c": len(cells), "v": len(surv_v) + len(placed), "e": len(edges)} + offset, next_point = {}, pStart for _start, key in sorted(((cS, "c"), (vS, "v"), (eS, "e"))): - offset[key] = at - at += sizes[key] + offset[key] = next_point + next_point += sizes[key] point_map = np.full(pEnd - pStart, -1, dtype=np.int64) point_map[surv_c - pStart] = offset["c"] + np.arange(len(surv_c)) point_map[surv_v - pStart] = offset["v"] + np.arange(len(surv_v)) point_map[np.asarray(surv_e, dtype=np.int64) - pStart] = ( offset["e"] + np.arange(len(surv_e))) + placed_points = offset["v"] + len(surv_v) + np.arange(len(placed)) def v_new(v): - return int(point_map[v - pStart]) + return (int(placed_points[-int(v) - 1]) if v < 0 + else int(point_map[int(v) - pStart])) new = PETSc.DMPlex().create(comm=dm.comm) new.setDimension(dm.getDimension()) - new.setChart(pStart, at) + new.setChart(pStart, next_point) for i in range(len(cells)): new.setConeSize(offset["c"] + i, 3) for i in range(len(edges)): @@ -731,13 +762,13 @@ def v_new(v): new.symmetrize() new.stratify() - _write_coordinates(new, dm, (offset["v"], offset["v"] + len(surv_v)), - surv_v - vS) + _write_coordinates(new, cdim, (offset["v"], offset["v"] + sizes["v"]), + np.vstack([X[surv_v - vS], placed.reshape(-1, cdim)])) _copy_labels(new, dm, point_map) if uw.mpi.size > 1: - _rebuild_point_sf(new, dm, point_map, at - pStart) - return new, point_map + _rebuild_point_sf(new, dm, point_map, next_point - pStart) + return new, point_map, placed_points def _rebuild_point_sf(new, dm, point_map, nroots): @@ -1172,7 +1203,7 @@ def remove_vertices(dm, candidates, max_passes=3, gate="both"): n = uw.mpi.comm.allreduce(len(victims), op=MPI.SUM) if n == 0: break - dm, point_map = rebuild_without_vertices(dm, victims, drop, made) + dm, point_map, _placed = rebuild_cavities(dm, victims, drop, made) cand = point_map[cand - pStart] cand = cand[cand >= 0] total += n diff --git a/tests/parallel/ptest_0844_reconnect_parallel.py b/tests/parallel/ptest_0844_reconnect_parallel.py index d545f6ad..616be947 100644 --- a/tests/parallel/ptest_0844_reconnect_parallel.py +++ b/tests/parallel/ptest_0844_reconnect_parallel.py @@ -182,7 +182,7 @@ def _sf_coordinate_drift(dm): Deletion compacts the point chart, so the star-forest cannot be reused verbatim the way a flip's can: every point after a deleted one shifts, and a leaf's *remote* index is a number only its owner holds. - ``rebuild_without_vertices`` renumbers locally and broadcasts the new + ``rebuild_cavities`` renumbers locally and broadcasts the new numbering once to close that gap. This is the check a mis-renumbering cannot pass and nothing else catches. @@ -233,7 +233,7 @@ def test_removal_renumbers_the_star_forest(): def test_removal_leaves_the_seam_alone(): """No shared point may be deleted, which is what keeps the leaf set intact. - ``rebuild_without_vertices`` renumbers the star-forest but does not rebuild + ``rebuild_cavities`` renumbers the star-forest but does not rebuild it, so a deleted shared point would leave a leaf pointing at nothing. The pass freezes any cavity touching the seam; this is that rule as a postcondition, checked by coordinates because the numbering has moved. diff --git a/tests/test_0848_place_surface.py b/tests/test_0848_place_surface.py new file mode 100644 index 00000000..996dd88c --- /dev/null +++ b/tests/test_0848_place_surface.py @@ -0,0 +1,382 @@ +"""Surfaces embedded by PLACING their points +(:mod:`underworld3.utilities.place_surface`). + +The cut (``tests/test_0844_line_cut.py``) makes a surface a chain of element +edges by splitting every edge it crosses. Placement makes the same chain by +asserting the surface's own points, deleting the mesh vertices in the way and +retriangulating the cavity. The two must deliver the same guarantees, so this +file asserts the same properties the cut's suite does — chain identity, no +straddling cell, exact placement, the label, area — and then the cases the cut +**refuses**, which are the reason the second implementation exists: + +- a surface **ending inside** the mesh. The cut has no split that represents a + triangle the surface enters and does not leave, and says so; +- two surfaces **closer together than one element**. They cross the same edges, + an edge can be split once, and the cut says so; +- a surface **finer than the local h**, which the cut cannot express at all + because its vertices are the crossings of the mesh's own edges. + +Three of these tests assert the refusal as well as the success. That pairing is +the point: a test that only showed placement working would not show that it is +buying anything. + +The walk that fills the cavity can fail, and the failure is loud. It is measured +rather than argued: 100 random traces on a uniform mesh and 100 on a graded +adapt-on-top mesh all place, with the total area exact to 2e-16 in every case. +""" +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.utilities import reconnect +from underworld3.utilities.line_cut import (cell_areas, cut_along_lines, + min_angles) +from underworld3.utilities.place_surface import place_along_lines + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + +# Both ends outside the mesh: the usual "specify it a little long" convention. +CROSSING = np.array([[-0.2, 0.317], [1.2, 0.683]]) +# One end outside, one inside: a fault reaching the wall and stopping. +ONE_TIP = np.array([[-0.2, 0.40], [0.55, 0.52]]) +# Both ends inside: a fault. The cut refuses this outright. +FAULT = np.array([[0.25, 0.45], [0.70, 0.56]]) + + +def _box(cell_size=1 / 16, **kwargs): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=cell_size, regular=False, qdegree=2, **kwargs) + + +def _coords(dm): + return np.asarray(dm.getCoordinatesLocal().array).reshape(-1, 2) + + +def _cell_vertex_indices(dm): + vS, vE = dm.getDepthStratum(0) + cS, cE = dm.getHeightStratum(0) + return np.array([[int(p) - vS for p in dm.getTransitiveClosure(c)[0] + if vS <= p < vE] for c in range(cS, cE)]) + + +def _surface_vertices(dm, name, value): + """Local indices of the vertices carrying the labelled facets.""" + vS, _vE = dm.getDepthStratum(0) + out = set() + for e in dm.getLabel(name).getStratumIS(value).getIndices(): + out.update(int(v) - vS for v in dm.getCone(e)) + return np.array(sorted(out)) + + +def _distance_to_trace(points, trace): + """Distance from each point to the trace, measured to its SEGMENTS. + + Not to the infinite line through its ends: a fault stops, and the mesh + beyond its tip is legitimately on both sides of that line. + """ + best = np.full(len(points), np.inf) + for A, B in zip(trace[:-1], trace[1:]): + d = B - A + u = np.clip(((points - A) @ d) / (d @ d), 0.0, 1.0) + best = np.minimum(best, np.linalg.norm(points - (A + u[:, None] * d), + axis=1)) + return best + + +def _conforming(dm): + """No facet with more than two cells, no inverted cell, Euler number 1.""" + nv = dm.getDepthStratum(0)[1] - dm.getDepthStratum(0)[0] + ne = dm.getDepthStratum(1)[1] - dm.getDepthStratum(1)[0] + nc = dm.getHeightStratum(0)[1] - dm.getHeightStratum(0)[0] + over = sum(1 for f in range(*dm.getHeightStratum(1)) + if len(dm.getSupport(f)) > 2) + return over == 0 and (cell_areas(dm) > 0.0).all() and nv - ne + nc == 1 + + +@pytest.mark.parametrize("trace", [CROSSING, ONE_TIP, FAULT], + ids=["crossing", "one-tip", "fault"]) +def test_the_three_cavity_shapes_all_place(trace): + """No end, one end or two ends on the wall: three topologies, one walk. + + How many ends reach the domain boundary decides whether the cavity is an + annulus, a disc, or two discs. All three go through the same walk, so all + three are asserted together — a change that broke one of them would + otherwise be caught only by whichever case happened to be tested. + """ + base = _box() + before = cell_areas(base.dm).sum() + dm, info = place_along_lines(base.dm, [trace], label="Fault", label_value=7) + + assert _conforming(dm) + assert cell_areas(dm).sum() == pytest.approx(before, rel=1e-13) + assert info["n_surface_facets"] == (info["n_placed"] + + info["n_on_surface"] - 1) + assert dm.getLabel("Fault").getStratumSize(7) == info["n_surface_facets"] + + +@pytest.mark.parametrize("trace", [CROSSING, ONE_TIP, FAULT], + ids=["crossing", "one-tip", "fault"]) +def test_the_surface_is_where_it_was_asked_for(trace): + """The placed points are vertices of the result, not near ones. + + This is the property placement has and snapping does not: the cut moves the + MESH onto the surface and so is exact too, but it can only choose from the + crossings the mesh offers. Here the surface's own points are the vertices, + so "on the surface" is not a tolerance — the only slack is where an end had + to be projected onto a boundary facet, which is rounding. + """ + dm, _info = place_along_lines(_box().dm, [trace], label="Fault", + label_value=7) + X = _coords(dm)[: dm.getDepthStratum(0)[1] - dm.getDepthStratum(0)[0]] + on = _surface_vertices(dm, "Fault", 7) + assert _distance_to_trace(X[on], trace).max() < 1e-12 + + +@pytest.mark.parametrize("trace", [CROSSING, ONE_TIP, FAULT], + ids=["crossing", "one-tip", "fault"]) +def test_no_cell_straddles_the_surface(trace): + """The property a cell-wise material needs in order to be CORRECT. + + Asserted at the source: every segment between consecutive placed points is + an element edge, so the surface is a union of element edges and no cell can + have part of it inside. Checked by walking the chain rather than by + signed distances, which cannot tell "beyond the tip" from "straddling". + """ + dm, info = place_along_lines(_box().dm, [trace], label="Fault", + label_value=7) + X = _coords(dm) + on = _surface_vertices(dm, "Fault", 7) + order = on[np.argsort((X[on] - trace[0]) @ (trace[-1] - trace[0]))] + + vS, _vE = dm.getDepthStratum(0) + edges = {frozenset(int(v) - vS for v in dm.getCone(e)): int(e) + for e in range(*dm.getDepthStratum(1))} + labelled = set(dm.getLabel("Fault").getStratumIS(7).getIndices()) + assert len(order) == info["n_surface_facets"] + 1 + for a, b in zip(order[:-1], order[1:]): + e = edges.get(frozenset((int(a), int(b)))) + assert e is not None, "a segment of the surface is not a mesh edge" + assert e in labelled, "a segment of the surface is not labelled" + + +def test_a_surface_ending_inside_the_mesh_is_placed_where_the_cut_refuses(): + """The headline case. A tip is a placed vertex; it is not a split. + + The cut has no representation for a triangle the surface enters and does not + leave, so it refuses. Placement closes the cavity round the tip with a fan, + which is the whole reason the tip's turn through 180 degrees is given a + window of the walk's parameter to itself. + """ + base = _box() + with pytest.raises(ValueError, match="entered but not left"): + cut_along_lines(base.dm, [FAULT]) + + dm, info = place_along_lines(base.dm, [FAULT], label="Fault", label_value=7) + assert _conforming(dm) + assert info["n_surface_facets"] == (info["n_placed"] + + info["n_on_surface"] - 1) + + # The tip really is a vertex of the mesh, and it is the END of the chain. + X = _coords(dm) + tip = np.flatnonzero(np.linalg.norm(X - FAULT[-1], axis=1) < 1e-13) + assert len(tip) == 1, "the tip is not a vertex of the result" + on = set(_surface_vertices(dm, "Fault", 7).tolist()) + assert int(tip[0]) in on + + +def test_two_surfaces_closer_than_one_element_are_placed(): + """An edge can be split once, so the cut refuses converging flanks. + + This is the restriction that makes a tapering fault unmeshable by cutting, + and it is a property of the METHOD rather than of the mesh: refining shrinks + the separation at which it bites but never removes it. Placement is not + competing for the mesh's edges at all, so the separation is free. + """ + base = _box() + h = 1 / 16 + + def pair(gap): + return [np.array([[-0.2, 0.5 + gap * h / 2], [1.2, 0.5 + gap * h / 2]]), + np.array([[-0.2, 0.5 - gap * h / 2], [1.2, 0.5 - gap * h / 2]])] + + with pytest.raises(ValueError, match="crossed more than once"): + cut_along_lines(base.dm, pair(0.25)) + + dm, info = place_along_lines(base.dm, pair(0.25), label="Pair", + label_value=9, clearance=0.8) + assert _conforming(dm) + # Two chains, so two fewer facets than vertices along them. + assert info["n_surface_facets"] == (info["n_placed"] + + info["n_on_surface"] - 2) + + +def test_the_surface_may_be_finer_than_the_mesh(): + """``spacing`` is the knob the cut does not have. + + A cut's surface vertices ARE the crossings of the mesh's edges, so its + resolution is the mesh's and cannot be asked for separately. A placed + surface carries its own point spacing, which is what lets the geometry be + resolved without refining everything around it. + """ + base = _box() + counts = [] + for factor in (1.0, 0.5, 0.25): + _dm, info = place_along_lines(base.dm, [FAULT], label="Fault", + label_value=7, spacing=factor / 16) + counts.append(info["n_surface_facets"]) + assert counts[1] > 1.5 * counts[0] and counts[2] > 1.5 * counts[1], ( + f"halving the spacing did not roughly double the facets: {counts}") + + +def test_the_base_mesh_is_not_modified(): + """The surface's position is a design variable: the base must survive it.""" + base = _box() + cells_before = base.dm.getHeightStratum(0)[1] - base.dm.getHeightStratum(0)[0] + coords_before = _coords(base.dm).copy() + + place_along_lines(base.dm, [CROSSING], label="A", label_value=7) + place_along_lines(base.dm, [FAULT], label="B", label_value=8) + + after = base.dm.getHeightStratum(0)[1] - base.dm.getHeightStratum(0)[0] + assert after == cells_before + assert np.array_equal(_coords(base.dm), coords_before) + + +def test_splitting_a_wall_facet_keeps_the_wall_labelled(): + """A surface reaching the boundary replaces one wall facet with two. + + The rebuild carries labels by point id, and the facet that was split has no + point id in the result, so without inheritance its ``Left`` / ``Top`` / + ``UW_Boundaries`` values would simply be dropped. The hole that leaves is + invisible until a boundary condition steps over it, which is why the count + is asserted rather than the mesh merely being looked at. + """ + base = _box(1 / 12) + # The four named walls only. ``Null_Boundary`` is the sentinel a natural + # condition attaches to when it applies to no boundary, and it marks every + # VERTEX of every UW3 mesh — so its stratum shrinks by however many vertices + # the cavity swallowed, which says nothing about the wall. + walls = {name: base.dm.getLabel(name).getStratumSize( + base.boundaries[name].value) + for name in ("Left", "Right", "Top", "Bottom")} + dm, _info = place_along_lines(base.dm, [CROSSING], label="Fault", + label_value=99) + + for name, before in walls.items(): + after = dm.getLabel(name).getStratumSize(base.boundaries[name].value) + assert after >= before, f"{name} lost facets: {before} -> {after}" + # The two walls the surface crosses each gained exactly one facet. + gained = sum(dm.getLabel(n).getStratumSize(base.boundaries[n].value) - b + for n, b in walls.items()) + assert gained == 2, f"a split wall facet was not replaced by two: {gained}" + + +def test_a_second_surface_does_not_damage_the_first(): + """Placed surfaces compose, because a labelled edge is an interface. + + The second placement clears mesh vertices out of its own way, and the first + surface's vertices are exactly the ones it must not take. They are protected + by the same rule the repair passes use — a labelled edge is an interface — + so this also checks that the label placement writes is the one that rule + reads. + """ + base = _box() + one, info_one = place_along_lines(base.dm, [CROSSING], label="Fault", + label_value=7) + two, _info_two = place_along_lines(one, [np.array([[-0.2, 0.12], + [1.2, 0.12]])], + label="Moho", label_value=8) + assert _conforming(two) + assert two.getLabel("Fault").getStratumSize(7) == info_one["n_surface_facets"] + assert two.getLabel("Moho").getStratumSize(8) > 0 + + +@pytest.mark.level_2 +def test_repair_improves_the_shapes_and_leaves_the_surface_alone(): + """The walk fills the cavity by parameter, not by shape; repair does shape. + + Flipping and deleting are the two operations the fill does not have, and + both refuse to act on a labelled edge — so the surface has to come through + with the same facets, which is asserted rather than assumed. + + A GRADED mesh, and the fixture is checked before the claim is made. On a + uniform mesh this same fault comes out of the walk with no cell under 15 + degrees at all, so a uniform fixture would assert that repair improved + something that was not wrong — true, and vacuous. + """ + line = np.array([[-0.1, 0.37], [1.1, 0.63]]) + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1 / 6, + regular=False, qdegree=2, refinement=2) + surf = uw.meshing.Surface("Grade", base, line) + surf.discretize() + graded = base.adapt(surf.refinement_metric_function( + h_near=1 / 48, h_far=1 / 6, width=1 / 12), max_levels=3) + + dm, info = place_along_lines(graded.dm, [line], label="Fault", label_value=7) + before = min_angles(dm) + assert int((before < 15).sum()) > 0, ( + "the fill left no thin cell on this mesh, so there is nothing for " + "repair to fix and this test would pass without exercising it") + + flipped, n_flips = reconnect.flip_to_reduce_max_angle(dm) + after = min_angles(flipped) + + assert n_flips > 0 + assert int((after < 15).sum()) < int((before < 15).sum()) + assert after.min() >= before.min() - 1e-12, "repair lowered the worst angle" + assert flipped.getLabel("Fault").getStratumSize(7) == info["n_surface_facets"] + assert cell_areas(flipped).sum() == pytest.approx(cell_areas(dm).sum(), + rel=1e-13) + + +def test_a_dirichlet_condition_applies_on_a_placed_surface(): + """A label is only useful if a solver actually constrains those DOFs. + + Also the end-to-end check that the rebuilt plex is a mesh UW3 can use at + all: sections, coordinates, boundaries and the solve. + """ + base = _box(1 / 12, refinement=1) + boundaries = base._boundaries_with("Fault") + dm, _info = place_along_lines(base.dm, [np.array([[0.5, -0.2], [0.5, 1.2]])], + label="Fault", + label_value=boundaries["Fault"].value) + mesh = uw.discretisation.Mesh( + dm, simplex=True, qdegree=3, boundaries=boundaries, + coordinate_system_type=base.CoordinateSystem.coordinate_type, + verbose=False) + + u = uw.discretisation.MeshVariable("u_placed", mesh, 1, degree=1) + poisson = uw.systems.Poisson(mesh, u_Field=u) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = 0.0 + for b in ("Left", "Right", "Top", "Bottom"): + poisson.add_dirichlet_bc(0.0, b) + poisson.add_dirichlet_bc(1.0, "Fault") + poisson.solve() + + X, values = np.asarray(u.coords), np.asarray(u.data[:, 0]) + on = np.abs(X[:, 0] - 0.5) < 1e-11 + assert on.sum() > 0 + assert np.allclose(values[on], 1.0, atol=1e-10), "the surface BC was not applied" + interior = (~on) & (X[:, 0] > 0.1) & (X[:, 0] < 0.4) + assert 0.0 < values[interior].max() < 1.0 + + +def test_a_surface_that_leaves_the_domain_and_returns_is_refused(): + """Two pieces are two surfaces, and they want two names.""" + base = _box() + out_and_back = np.array([[0.2, 0.5], [0.5, 1.2], [0.8, 0.5]]) + with pytest.raises(ValueError, match="leaves the domain and re-enters"): + place_along_lines(base.dm, [out_and_back], label="Fault", label_value=7) + + +def test_three_dimensions_is_refused_with_the_reason(): + """A placed sheet's cavity is a polyhedron, and filling one is not this.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), + cellSize=0.5, qdegree=2) + with pytest.raises(NotImplementedError, match="Steiner"): + place_along_lines(mesh.dm, [CROSSING], label="Fault", label_value=7) From 8a79e3dbb3958816ea9ac322c508846219e8d2ab Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 5 Aug 2026 11:49:15 +1000 Subject: [PATCH 2/3] Hold an embedded surface's cells: a second placement was eating the first Placing a surface against one already embedded destroyed it, silently. The cavity protects interface VERTICES from deletion, but a cell is not a vertex: both cells supporting an interface facet could be cleared while every one of their corners was protected. The facet then has no support left, the refill has no reason to recreate that edge, and the earlier surface loses a facet out of the middle of its chain. Measured on a T junction: a trunk of 21 facets came back with 20, the junction vertex carried only the branch's label, and nothing raised. Three changes. Cells owning an INTERIOR labelled facet are held out of the cavity - interior, because the domain's own walls carry edge labels too and holding their cells would forbid clearing anything against a wall, which is what a surface crossing the domain must do. A vertex may only be deleted if every cell of its star can be cleared, so the cavity stays the union of its victims' stars. And every placement re-reads each earlier surface's facet count off the RESULT mesh and refuses if one dropped; the label being written may grow, since several polylines may share a name, but no other may change. This retracts a capability claim rather than adding one. The suite asserted that two surfaces a tenth of an element apart could be placed - an identity summed over each placement's own counts, which a partial corruption still satisfies. Re-measured with both surfaces checked intact afterwards, on a 1/16 box: placing one at a time accepts 1.5 h separation and refuses 1.0 h, while the cut accepts 1.0 h and refuses 0.5 h. For closely spaced surfaces the cut is currently the more capable of the two, and the docs and module docstring now say so. The two limits are not the same kind, which is what survives of the argument. The cut's is inherent: converging flanks cross the same edge and an edge splits at one point. This one is an implementation limit - one cavity holds one surface - and lifting it means placing both into a single cavity, which is the finite-width ribbon and is not built. Junctions now refuse loudly at every resolution tested (T, Y and X alike) instead of T corrupting quietly. Interior end-snapping is what will turn that refusal into an abutment. Underworld development team with AI support from Claude Code --- .../conforming-surfaces-and-fault-zones.md | 26 +++- src/underworld3/utilities/place_surface.py | 131 +++++++++++++++--- tests/test_0848_place_surface.py | 98 ++++++++++--- 3 files changed, 214 insertions(+), 41 deletions(-) diff --git a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md index 5a94fffa..15adf0ee 100644 --- a/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md +++ b/docs/developer/subsystems/conforming-surfaces-and-fault-zones.md @@ -227,15 +227,31 @@ them applies to placing: | | cut | place | |---|---|---| | a surface ending **inside** the mesh | refused — a triangle entered and not left has no split that represents it | the tip is a placed vertex like any other, and the cavity closes round it | -| two surfaces **closer than one element** | refused — they cross the same edge, and an edge can be split once | free; neither is competing for an edge of the original mesh | | surface **finer than the local `h`** | impossible — the surface's vertices *are* the mesh's crossings | `spacing` is a parameter | +| two surfaces **closer than one element** | refused (inherent) | refused (implementation) — see below | | what it does to the mesh | moves vertices onto the surface (`snap_frac`) and splits edges | deletes vertices near the surface and refills the hole | | parallel | yes, and partition-independent | **serial only** — see below | -Measured on the same box, viscosity contrast irrelevant: the cut accepts two -parallel surfaces one element apart and refuses them at half an element, while -placement carries them down to a tenth of an element with a worst angle of 4.2 -degrees. +```{warning} +An earlier version of this section claimed placement carried two parallel +surfaces down to a tenth of an element. That was wrong. The second surface's +cavity was consuming the first, and the identity being asserted summed each +placement's own facet counts, so it could not see the damage. Read such counts +back off the RESULT mesh. + +Measured on a 1/16 box, both surfaces checked intact afterwards: placing one at +a time accepts 1.5 `h` separation and refuses 1.0 `h`; the cut accepts 1.0 `h` +and refuses 0.5 `h`. **For closely spaced surfaces the cut is currently the more +capable of the two.** + +The two limits are not the same kind. The cut's is inherent — converging flanks +cross the same edge and an edge can be split at one point, so no implementation +removes it. Placement's is an implementation limit: a cavity is cleared and +filled for one surface, and the cells carrying an earlier surface's facets are +held back so that surface survives, which eventually leaves no room. Lifting it +means placing both surfaces into a single cavity — the finite-width ribbon, +which is not built. +``` The construction is the same operation in 2-D on a curve and in 3-D on a sheet — place, delete, refill — which matters because cutting tetrahedra along a surface diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index d756b938..cea556d1 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -16,9 +16,19 @@ * a tip terminates *inside* the mesh, because the tip is a placed vertex like any other and the cavity closes around it; * the surface's point spacing is chosen by the caller, not by wherever the mesh - happened to put its edges; -* two surfaces may run arbitrarily close together, because neither is competing - for an edge of the original mesh. + happened to put its edges. + +Not, however, two surfaces closer together than a cell. That was claimed here +once and it was wrong: the second surface's cavity was consuming the first, and +the identity being asserted summed each placement's own counts and so could not +see it. A cavity is cleared and filled for ONE surface, and the cells carrying +an earlier surface's facets are held back so that surface survives — which +eventually leaves no room. Measured on a 1/16 box: placing one at a time accepts +1.5 h separation and refuses 1.0 h, where the cut accepts 1.0 h and refuses +0.5 h. The two limits are not the same KIND of limit — the cut's is inherent +(converging flanks cross one edge, and an edge splits at one point) while this +one is an implementation limit, lifted by placing both surfaces into a single +cavity, which is the finite-width ribbon and is not built. The operation is the same on a curve in 2-D and on a sheet in 3-D — place, delete, retriangulate — which is the other reason to prefer it: cutting @@ -133,24 +143,72 @@ def _boundary_vertices(dm, n_vertices): return on -def _interface_vertices(dm, n_vertices): - """Vertices carrying an edge of an already-embedded surface. - - Deleting one would put a gap in a surface placed earlier, so these are - refused as victims. Read through :func:`reconnect._interface_edges`, which - already knows which labels mean *interface* and which are PETSc's or UW3's - own bookkeeping — a distinction that has produced three separate silent - failures when guessed at. +def _interface_vertices_and_cells(dm, n_vertices, n_cells): + """Vertices and cells that an already-embedded surface may not lose. + + Two masks, and the second is the one that is easy to forget. Protecting the + interface's VERTICES from deletion does not protect its FACETS: a cell is + not a vertex, and both cells supporting an interface facet can be cleared + while every one of their corners is protected. The facet then has no support + left, the refill has no reason to recreate that edge, and the earlier + surface loses a facet out of the middle of its chain — silently. Measured on + a T junction: a trunk of 21 facets came back with 20 and nothing raised. + + Read through :func:`reconnect._interface_edges`, which already knows which + labels mean *interface* and which are PETSc's or UW3's own bookkeeping — a + distinction that has produced three separate silent failures when guessed at. + Narrowed here to INTERIOR facets, which is what "a surface already embedded" + means: the domain's own walls carry edge labels too, and holding their cells + would forbid clearing anything against a wall — which is exactly what a + surface crossing the domain has to do. """ pStart, _pEnd = dm.getChart() vS, _vE = dm.getDepthStratum(0) + cS, _cE = dm.getHeightStratum(0) locked = reconnect._interface_edges(dm) - on = np.zeros(n_vertices, dtype=bool) - for e in range(*dm.getDepthStratum(1)): - if locked[e - pStart]: - for v in dm.getCone(e): - on[int(v) - vS] = True - return on + vertices = np.zeros(n_vertices, dtype=bool) + cells = np.zeros(n_cells, dtype=bool) + for e in _interior_interface_facets(dm, locked, pStart): + for v in dm.getCone(e): + vertices[int(v) - vS] = True + for c in dm.getSupport(e): + cells[int(c) - cS] = True + return vertices, cells + + +def _interior_interface_facets(dm, locked, pStart): + """The labelled facets that lie INSIDE the mesh — an embedded surface.""" + return [e for e in range(*dm.getDepthStratum(1)) + if locked[e - pStart] and len(dm.getSupport(e)) == 2] + + +def _interface_facet_counts(dm): + """How many INTERIOR facets each label holds, so a breach can be detected. + + Interior only, for the reason :func:`_interface_vertices_and_cells` gives: + a wall label's count legitimately changes when a surface's end splits a wall + facet, and comparing that would refuse the ordinary case. + """ + pStart, _pEnd = dm.getChart() + locked = reconnect._interface_edges(dm) + interior = set(_interior_interface_facets(dm, locked, pStart)) + counts = {} + for i in range(dm.getNumLabels()): + name = dm.getLabelName(i) + if name in reconnect._TOPOLOGY_LABELS: + continue + label = dm.getLabel(name) + values = label.getValueIS() + if values is None: + continue + for val in values.getIndices(): + if label.getStratumSize(int(val)) == 0: + continue + held = [p for p in label.getStratumIS(int(val)).getIndices() + if int(p) in interior] + if held: + counts[(name, int(val))] = len(held) + return counts # ------------------------------------------------------------------ the polyline @@ -697,7 +755,15 @@ def _place_one(dm, pts, label, label_value, clearance, spacing, end_snap): # Which vertices are in the way. A domain-boundary vertex is never one: # deleting it would change the domain. Nor is a vertex of a surface already # embedded, which would put a gap in that surface. - protected = on_boundary | _interface_vertices(dm, len(X)) + held_v, held_c = _interface_vertices_and_cells(dm, len(X), len(cells)) + held_counts = _interface_facet_counts(dm) + # A vertex may only be deleted if every cell of its star may be cleared. + # Held cells stay, so a vertex beside a surface already embedded would end + # up deleted while still on the cavity's boundary — a cavity that is not the + # union of its victims' stars, which the walk cannot fill. + beside_held = np.zeros(len(X), dtype=bool) + beside_held[cells[held_c].ravel()] = True + protected = on_boundary | held_v | beside_held victim = ((_distance_to_lines(X, [pts]) < clearance * _vertex_h(X, _edge_vertices(dm))) & ~protected) @@ -733,6 +799,14 @@ def _place_one(dm, pts, label, label_value, clearance, spacing, end_snap): # the domain wall, or on a surface already embedded. for _attempt in range(8): drop = np.union1d(np.flatnonzero(victim[cells].any(axis=1)), crossed) + # A cell owning a facet of a surface already embedded is never cleared. + # Clearing BOTH cells of such a facet destroys it — the facet's support + # is gone, so the refill has no reason to recreate that edge — and the + # earlier surface loses a facet out of the middle of its chain without + # anything raising. Holding the cells instead stops the cavity at the + # earlier surface, which is also what makes the ligament of an offset + # junction survive. + drop = drop[~held_c[drop]] if not len(drop): raise ValueError( "the surface meets no cell of this mesh: there is nothing to " @@ -798,6 +872,27 @@ def new_point(v): _inherit_boundary_labels( new_dm, dm, [(e, new_point(a), new_point(int(index[k])), new_point(b)) for k, (e, a, b) in split_at.items()]) + + # Every surface already embedded must come through with the facets it had. + # Held cells are what makes that true; this is the check that it IS true, + # and it is not redundant — losing a facet of an earlier surface produces a + # mesh that passes conformity, area and orientation, and is wrong only in + # the one place nothing else looks. + after = _interface_facet_counts(new_dm) + for key, before in held_counts.items(): + now = after.get(key, 0) + # The label being written may GROW — several polylines may share a name, + # which is how a fault with more than one segment is labelled. Any other + # label must come through with exactly the facets it had. + if now < before or (now != before and key != (label, int(label_value))): + raise RuntimeError( + f"placing {label!r} would leave the surface {key[0]!r} with " + f"{now} facets instead of {before}: this surface's cavity " + "reached one already embedded. Surfaces must be separated by " + "at least a cell when placed one at a time; two that run closer " + "than that have to be placed together, into one cavity, which " + "is not yet implemented.") + return new_dm, {"n_placed": len(placed), "n_on_surface": len(pts) - len(placed), "n_removed": int(victim.sum()), diff --git a/tests/test_0848_place_surface.py b/tests/test_0848_place_surface.py index 996dd88c..4333455b 100644 --- a/tests/test_0848_place_surface.py +++ b/tests/test_0848_place_surface.py @@ -185,30 +185,59 @@ def test_a_surface_ending_inside_the_mesh_is_placed_where_the_cut_refuses(): assert int(tip[0]) in on -def test_two_surfaces_closer_than_one_element_are_placed(): - """An edge can be split once, so the cut refuses converging flanks. +def _parallel_pair(gap, h=1 / 16): + return [np.array([[-0.2, 0.5 + gap * h / 2], [1.2, 0.5 + gap * h / 2]]), + np.array([[-0.2, 0.5 - gap * h / 2], [1.2, 0.5 - gap * h / 2]])] - This is the restriction that makes a tapering fault unmeshable by cutting, - and it is a property of the METHOD rather than of the mesh: refining shrinks - the separation at which it bites but never removes it. Placement is not - competing for the mesh's edges at all, so the separation is free. - """ - base = _box() - h = 1 / 16 - def pair(gap): - return [np.array([[-0.2, 0.5 + gap * h / 2], [1.2, 0.5 + gap * h / 2]]), - np.array([[-0.2, 0.5 - gap * h / 2], [1.2, 0.5 - gap * h / 2]])] +def test_two_surfaces_placed_one_at_a_time_survive_each_other(): + """Well separated, both come through with every facet they were given. + Read the counts back off the RESULT, not off the returned info. The info is + accumulated as each surface is placed, so it reports what was labelled at + the time and cannot see the second cavity eating the first — which is + exactly how a capability that never worked was reported as working. Any + identity that a partial corruption still satisfies is not a test. + """ + base = _box() + first, info_one = place_along_lines(base.dm, [_parallel_pair(3.0)[0]], + label="A", label_value=9) + both, info_two = place_along_lines(first, [_parallel_pair(3.0)[1]], + label="B", label_value=10) + + assert _conforming(both) + assert both.getLabel("A").getStratumSize(9) == info_one["n_surface_facets"], ( + "the second surface ate facets from the first") + assert both.getLabel("B").getStratumSize(10) == info_two["n_surface_facets"] + + +def test_surfaces_closer_than_a_cell_are_refused_rather_than_merged(): + """One cavity holds one surface, so placing them one at a time has a limit. + + The two methods reach their limits for different reasons, and the difference + matters for where each can be taken. The CUT's limit is inherent: converging + flanks cross the same edge and an edge can be split at one point, so no + amount of implementation removes it. Placement's is an implementation limit + — a cavity is cleared and filled for ONE surface, and the cells carrying an + earlier surface's facets are held back so that surface survives, which + eventually leaves no room. Lifting it means placing both into one cavity, + which is the finite-width ribbon and is not built. + + Measured on the 1/16 box: the cut accepts 1.0 h and refuses 0.5 h; placing + one at a time accepts 1.5 h and refuses 1.0 h. Placement is currently the + MORE restrictive of the two here, and this test pins that rather than the + opposite claim, which an earlier version of this file asserted on an + identity that could not see the damage. + """ + base = _box() with pytest.raises(ValueError, match="crossed more than once"): - cut_along_lines(base.dm, pair(0.25)) + cut_along_lines(base.dm, _parallel_pair(0.5)) - dm, info = place_along_lines(base.dm, pair(0.25), label="Pair", + first, _ = place_along_lines(base.dm, [_parallel_pair(0.5)[0]], label="A", label_value=9, clearance=0.8) - assert _conforming(dm) - # Two chains, so two fewer facets than vertices along them. - assert info["n_surface_facets"] == (info["n_placed"] - + info["n_on_surface"] - 2) + with pytest.raises(RuntimeError): + place_along_lines(first, [_parallel_pair(0.5)[1]], label="B", + label_value=10, clearance=0.8) def test_the_surface_may_be_finer_than_the_mesh(): @@ -292,6 +321,39 @@ def test_a_second_surface_does_not_damage_the_first(): assert two.getLabel("Moho").getStratumSize(8) > 0 +def test_a_surface_placed_against_an_existing_one_cannot_break_it(): + """A second surface reaching an existing one must not eat its facets. + + Protecting interface VERTICES from deletion — which placement does — is not + enough. A cell is not a vertex: both cells supporting an interface facet can + be cleared while every one of their corners is protected, and then the facet + has no support left, the refill has no reason to recreate that edge, and the + first surface loses a facet in the middle of its chain. + + Measured before the guard existed, on the T fixture below: the trunk went 21 + facets to 20 and the junction vertex came out carrying only the branch's + label. Nothing raised. That is the failure mode this whole construction is + meant to prevent — a mesh that looks plausible and is quietly wrong — so it + is checked on every placement, not just here. + """ + base = _box(1 / 20) + trunk = np.array([[-0.2, 0.34], [1.2, 0.66]]) + branch = np.array([[0.5, 0.5], [0.62, 1.2]]) + + one, info = place_along_lines(base.dm, [trunk], label="F0", label_value=20) + assert one.getLabel("F0").getStratumSize(20) == info["n_surface_facets"] + + # Either the branch places and the trunk is untouched, or it is refused. + # What must never happen is placing while the trunk quietly loses facets. + try: + two, _ = place_along_lines(one, [branch], label="F1", label_value=21) + except RuntimeError as exc: + assert "already embedded" in str(exc) or "would lose" in str(exc) + return + assert two.getLabel("F0").getStratumSize(20) == info["n_surface_facets"], ( + "the trunk lost facets to the branch's cavity") + + @pytest.mark.level_2 def test_repair_improves_the_shapes_and_leaves_the_surface_alone(): """The walk fills the cavity by parameter, not by shape; repair does shape. From 54fc61fc3bc6c32c48f78430d593fbe37f49257d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 9 Aug 2026 08:45:33 +1000 Subject: [PATCH 3/3] Reconcile reconnect.py with the fault-split-node stream; renumber the placement test The two fault streams evolved utilities/reconnect.py independently after the 6d494db1 branch point: the contact stream added _install_point_sf (a rebuilt star-forest must also be installed on the coordinate DM, or a parallel HDF5 save writes every shared vertex as owned on every rank and the serial reload dies in coordinatesLoad) while this stream renamed rebuild_without_vertices to rebuild_cavities and taught it to add placed vertices as well as delete. This commit folds their three hunks in verbatim - the helper, the two call sites in _rebuild_point_sf, and the docstring sentence - so ONE version of the file exists and their branch's later merge diff on it is exactly the rename adaptation and nothing else. rebuild_with_cones' verbatim star-forest handoff is routed through the same helper: it creates its coordinate DM before the SF is installed, so it carries the same trap, and leaving one of three sites on the old pattern is how the next divergence starts. The placement test moves from test_0848 to test_0853: 0848 was claimed twice across the branches (their test_0848_fault_split_3d.py is pushed), and the number registry in the planning-hub liaison ledger now records 0848 as the contact stream's and 0853 as this stream's. Cross-branch takings, per the ledger rules: the _install_point_sf hunks are from c680e010 / c8693579 / 1d487319 on feature/fault-split-node. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/reconnect.py | 26 ++++++++++++++++--- ..._surface.py => test_0853_place_surface.py} | 0 2 files changed, 22 insertions(+), 4 deletions(-) rename tests/{test_0848_place_surface.py => test_0853_place_surface.py} (100%) diff --git a/src/underworld3/utilities/reconnect.py b/src/underworld3/utilities/reconnect.py index 6d5d6a4d..9b7e14e0 100644 --- a/src/underworld3/utilities/reconnect.py +++ b/src/underworld3/utilities/reconnect.py @@ -581,7 +581,7 @@ def rebuild_with_cones(dm, new_cells, new_edges): # The star-forest transfers verbatim: every rank preserves its numbering, so # the remote point numbers it carries are still the right ones. if uw.mpi.size > 1: - new.setPointSF(dm.getPointSF()) + _install_point_sf(new, dm.getPointSF()) return new @@ -771,6 +771,22 @@ def v_new(v): return new, point_map, placed_points +def _install_point_sf(new, new_sf): + """Install a rebuilt star-forest on the plex AND its coordinate DM. + + The coordinate DM is created when the rebuilt chart's coordinates + are written — BEFORE any point SF exists — so it snapshots an empty + one. The solve never notices (field sections are created later, + from the plex SF), but a parallel HDF5 save then writes every + shared vertex as owned on every rank, and the serial reload of such + a checkpoint dies in ``coordinatesLoad`` (measured: a split mesh + written at np = 4 carried owned+shared = 7121 vertex rows where the + true owned count is 6334). + """ + new.setPointSF(new_sf) + new.getCoordinateDM().setPointSF(new_sf) + + def _rebuild_point_sf(new, dm, point_map, nroots): """Carry the point star-forest onto a renumbered chart, in one exchange. @@ -780,7 +796,9 @@ def _rebuild_point_sf(new, dm, point_map, nroots): delivers exactly that, one value per leaf. The leaf set itself is unchanged: :func:`remove_vertices` never deletes a - shared point, so every leaf still exists and only its number has moved. + shared point, and the split-fault rebuilds never duplicate one (a + seam-touching fault is redistributed onto one rank before splitting), so + every leaf still exists and only its number has moved. """ pStart, pEnd = dm.getChart() sf = dm.getPointSF() @@ -802,7 +820,7 @@ def _rebuild_point_sf(new, dm, point_map, nroots): if ilocal is None or not len(ilocal): new_sf.setGraph(nroots, np.zeros(0, dtype=PETSc.IntType), np.zeros(0, dtype=PETSc.IntType)) - new.setPointSF(new_sf) + _install_point_sf(new, new_sf) return leaves = np.asarray(ilocal, dtype=np.int64) @@ -817,7 +835,7 @@ def _rebuild_point_sf(new, dm, point_map, nroots): remote[:, 0] = np.asarray(iremote).reshape(-1, 2)[:, 0] remote[:, 1] = remote_index new_sf.setGraph(nroots, local.astype(PETSc.IntType), remote.reshape(-1)) - new.setPointSF(new_sf) + _install_point_sf(new, new_sf) # ---------------------------------------------------------------- the flip pass diff --git a/tests/test_0848_place_surface.py b/tests/test_0853_place_surface.py similarity index 100% rename from tests/test_0848_place_surface.py rename to tests/test_0853_place_surface.py