From 09103bf82ee2b88bb7194355a1577adb2763b323 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 09:19:37 +1000 Subject: [PATCH] 2-D fault zones reach the domain boundary: the ribbon outcrops (#549) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The 3-D outcrop machinery brought one dimension down. The ribbon assembly is clipped to the domain box in OCC under the same specify-long contract (polylines may run past the domain), skin edges left exactly in one wall line become the outcrop band, and the carve opens onto that wall. The 2-D case is simpler than the 3-D one that ships: the clipped skin is still a closed loop, so the band splits it into band + open interior chain, and the cap degenerates from a wall-annulus meshing problem to two straight wall segments — the outcrop fill is one simply-connected polygon, the cavity ring with its wall span replaced by the interior chain (_outcrop_ring_splice). The splice segments and the band are relabelled with what the deleted wall span carried, full closures included, gated by an exact count. Refusals keep their teeth: an interior ribbon whose cavity reaches the wall still refuses, a zone meeting two walls refuses (box-edge outcrop bands are not built, matching 3-D), and a ribbon meeting one wall in two bands refuses. Wall vertices become deletable only in the open wall's line, and never at a domain corner. Tests: the band census with its interior-twin negative control, the P2 Poisson quadratic-exact oracle on the outcropped mesh, refusal coverage, and np=2 parallel forms (band present, no top-wall edge left unlabelled, refusals collective; also run at np=4). Underworld development team with AI support from Claude Code --- src/underworld3/utilities/place_surface.py | 302 ++++++++++++++++-- .../ptest_0855_place_thin_volume_parallel.py | 69 ++++ tests/test_0855_place_thin_volume.py | 118 +++++++ 3 files changed, 463 insertions(+), 26 deletions(-) diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index 201a1707..de377cfc 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -2801,21 +2801,24 @@ def _assembly_skin(points, cells): return points[node_ids], skin_local, node_ids -def _split_skin_band(skin_xyz, skin_tris, box_lo, box_hi): +def _split_skin_band(skin_xyz, skin_facets, box_lo, box_hi): """Split a skin into its interior part and the wall BAND — the strip of the assembly's clipped face lying exactly in a wall plane (the zone's outcrop). One wall only; more refuses (box-edge bands are not built). Returns ``(interior_idx, band_idx, wall_code)`` with wall_code None when nothing touches a wall. + + Dimension-free: facets are triangles against wall planes in 3-D, edges + against wall lines in 2-D. """ - codes = np.full(len(skin_tris), -1, dtype=np.int64) - for axis in range(3): + codes = np.full(len(skin_facets), -1, dtype=np.int64) + for axis in range(skin_xyz.shape[1]): for side, value in ((0, box_lo[axis]), (1, box_hi[axis])): - on = (skin_xyz[skin_tris][:, :, axis] == value).all(axis=1) + on = (skin_xyz[skin_facets][:, :, axis] == value).all(axis=1) codes[on] = 2 * axis + side walls = set(int(c) for c in codes[codes >= 0]) if not walls: - return np.arange(len(skin_tris)), np.empty(0, dtype=np.int64), None + return np.arange(len(skin_facets)), np.empty(0, dtype=np.int64), None if len(walls) > 1: raise NotImplementedError( "the zone meets more than one domain wall; box-edge outcrop " @@ -3038,7 +3041,7 @@ def cap_tag_of(k): gmsh.finalize() -def _occ_assembly_2d(polylines, width, size, assembly="fuse"): +def _occ_assembly_2d(polylines, width, size, assembly="fuse", box=None): """Thicken each polyline into a ribbon, resolve overlaps, mesh. The 2-D thin volume: a ribbon is the mitred outline of one polyline, and @@ -3048,6 +3051,13 @@ def _occ_assembly_2d(polylines, width, size, assembly="fuse"): the same region; they differ in the internal seams the mesher must honour. Returns ``(points, triangles, cad_area)``, the area being that of the resolved faces — the union — under either choice. + + ``box = (lo, hi)`` applies the specify-long contract one dimension down + from :func:`_occ_assembly_3d`: the resolved faces are INTERSECTED with + the domain rectangle, so polylines may protrude — a ribbon reaching the + top surface leaves its clipped edge exactly in the wall line (snapped to + the line value after meshing, defensively). ``cad_area`` is then the + clipped area, so the meshed-vs-CAD gate holds unchanged. """ import gmsh @@ -3115,6 +3125,13 @@ def outline(P): occ.fuse([(2, surfs[0])], [(2, t) for t in surfs[1:]]) else: occ.fragment([(2, surfs[0])], [(2, t) for t in surfs[1:]]) + if box is not None: + lo, hi = (np.asarray(b, dtype=float) for b in box) + occ.synchronize() + faces = [t for _d, t in gmsh.model.getEntities(2)] + tool = occ.addRectangle(lo[0], lo[1], 0.0, + hi[0] - lo[0], hi[1] - lo[1]) + occ.intersect([(2, t) for t in faces], [(2, tool)]) occ.synchronize() faces = gmsh.model.getEntities(2) @@ -3136,6 +3153,15 @@ def outline(P): dtype=np.int64).reshape(-1, 3)) if not tris: raise RuntimeError("the ribbon assembly meshed to no triangles") + if box is not None: + # OCC's clipped edges sit within rounding of the wall line; the + # band logic and the wall relabel need EXACT line values, so + # snap defensively (the 3-D path does the same on its planes). + lo, hi = (np.asarray(b, dtype=float) for b in box) + for axis in range(2): + for value in (lo[axis], hi[axis]): + near = np.abs(xy[:, axis] - value) < 1e-9 + xy[near, axis] = value return xy, np.vstack(tris), float(cad_area) finally: gmsh.finalize() @@ -3183,6 +3209,95 @@ def _skin_loops(skin_edges): return loops +def _outcrop_chain_2d(loops, band_pairs): + """Split the outcropping skin loop into its wall band and interior chain. + + The clipped skin is still a set of closed loops — clipping puts part of + a loop ON the wall line, it does not open the loop. The fill, though, + cannot take the outcropping loop as a hole: its boundary role is played + by the INTERIOR CHAIN, the open path left when the band's edges are + removed. ``band_pairs`` is that band as a set of frozen vertex pairs. + + Exactly one loop may carry band edges, and they must form ONE contiguous + arc of it — a ribbon meeting the wall in two separate bands is the 2-D + face of the box-edge case the 3-D band split refuses. Returns + ``(chain, holes)``: the chain as an open vertex path from one band + endpoint to the other, in the loop's own orientation, and the loops + without band edges unchanged. + """ + chain, holes = None, [] + for loop in loops: + n = len(loop) + on_band = [frozenset((loop[i], loop[(i + 1) % n])) in band_pairs + for i in range(n)] + if not any(on_band): + holes.append(loop) + continue + if chain is not None: + raise NotImplementedError( + "two skin loops meet the domain wall; only one ribbon may " + "outcrop.") + starts = sum(1 for i in range(n) + if on_band[i] and not on_band[i - 1]) + if starts != 1: + raise NotImplementedError( + "the ribbon meets the wall in more than one band; a " + "multiply-outcropping zone is not built.") + i0 = next(i for i in range(n) if on_band[i] and not on_band[i - 1]) + k = sum(on_band) + chain = [loop[(i0 + k + j) % n] for j in range(n - k + 1)] + return chain, holes + + +def _outcrop_ring_splice(ring, X, open_wall, chain_ids, chain_t): + """Replace the cavity ring's wall span with the ribbon's interior chain. + + The raw ring of an outcropping carve runs ALONG the open wall through + vertices about to be deleted. The fill's boundary instead descends + around the ribbon: the ring's single contiguous run of wall edges is + removed and the chain is spliced between the run's surviving end + vertices, oriented to meet them. The two splice segments are the 2-D + cap — what remains of the 3-D wall annulus one dimension down. + + ``chain_ids`` are the chain's rows in the fill's combined numbering, + ``chain_t`` the along-wall coordinates of its two ends. Returns + ``(spliced_ring, removed_wall_pairs)``, the removed pairs as old vertex + rows for the wall-label discovery. A second wall run refuses — the + carve spilled onto the wall away from the outcrop. + """ + axis, value = open_wall + t = 1 - axis + n = len(ring) + onw = [X[v][axis] == value for v in ring] + wall_edge = [onw[i] and onw[(i + 1) % n] for i in range(n)] + if not any(wall_edge): + raise RuntimeError( + "the outcrop band does not meet the cavity's wall span; raise " + "`clearance` so the carve reaches the wall.") + starts = sum(1 for i in range(n) if wall_edge[i] and not wall_edge[i - 1]) + if starts != 1: + raise RuntimeError( + "the carve reached the domain wall in two separate spans; " + "reduce `clearance` or move the zone off the wall.") + i0 = next(i for i in range(n) if wall_edge[i] and not wall_edge[i - 1]) + k = sum(wall_edge) + corner_l, corner_r = ring[i0], ring[(i0 + k) % n] + lo_t, hi_t = sorted((float(X[corner_l][t]), float(X[corner_r][t]))) + if not (lo_t < min(chain_t) and max(chain_t) < hi_t): + raise RuntimeError( + "the cavity's wall span does not cover the outcrop band; raise " + "`clearance`.") + removed = [(ring[(i0 + j) % n], ring[(i0 + j + 1) % n]) + for j in range(k)] + # The surviving arc, corner_r around to corner_l, then the chain with + # its corner_l-side end first — the loop's orientation is preserved. + tail = [ring[(i0 + k + j) % n] for j in range(n - k + 1)] + tl = float(X[corner_l][t]) + seq = (chain_ids if abs(tl - chain_t[0]) <= abs(tl - chain_t[1]) + else chain_ids[::-1]) + return tail + list(seq), removed + + def _ring_growing(cells, drop, held_mask): """The cavity ring, growing the drop set at pinch vertices until simple. @@ -3839,17 +3954,32 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, Serial AND parallel through the same gather-first mechanism as the 3-D volume: the assembly is meshed once (rank 0) and broadcast, the region gathers to one rank, the carve and the holes fill run there, every rank - rebuilds collectively. Ribbons are interior by construction, so the 2-D - line path's wall-end restriction does not arise. + rebuilds collectively. + + The specify-long contract holds as it does in 3-D: polylines may run + past the domain, the assembly is clipped to the box, and a clipped edge + left in a wall line is the ribbon's OUTCROP BAND — boundary edges + carrying both the skin label and the wall's labels. The 2-D cap is two + straight wall segments (the 3-D wall annulus one dimension down), so the + outcrop fill is one simply-connected polygon: the cavity ring with its + wall span replaced by the skin's interior chain. """ comm = uw.mpi.comm + # The (axis-aligned) domain box, collectively — the clip target and the + # wall lines the band is identified against. + Xb = _coords(dm) + lo_hi = np.array([Xb.min(axis=0) if len(Xb) else np.full(2, np.inf), + -(Xb.max(axis=0)) if len(Xb) else np.full(2, np.inf)]) + comm.Allreduce(MPI.IN_PLACE, lo_hi, op=MPI.MIN) + box_lo, box_hi = lo_hi[0], -lo_hi[1] + failure = None payload = None if comm.rank == 0: try: - asm_pts, asm_tris, cad_area = _occ_assembly_2d(polylines, width, - size, assembly) + asm_pts, asm_tris, cad_area = _occ_assembly_2d( + polylines, width, size, assembly, box=(box_lo, box_hi)) P = asm_pts[asm_tris] twice = ((P[:, 1, 0] - P[:, 0, 0]) * (P[:, 2, 1] - P[:, 0, 1]) - (P[:, 1, 1] - P[:, 0, 1]) * (P[:, 2, 0] - P[:, 0, 0])) @@ -3873,6 +4003,23 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, for a, b in skin_local] loops_asm = _skin_loops(skin_edges) + # The outcrop band: skin edges lying exactly in ONE wall line. The + # outcropping loop splits into band + interior chain; loops away from + # the wall stay holes of the fill. + _int_idx, band_idx, wall_code = _split_skin_band( + asm_pts, np.asarray(skin_edges, dtype=np.int64), box_lo, box_hi) + open_wall = None + chain_asm = None + hole_loops = loops_asm + band_pairs = set() + if wall_code is not None: + w_axis, w_side = divmod(int(wall_code), 2) + open_wall = (w_axis, + float(box_hi[w_axis] if w_side else box_lo[w_axis])) + band_pairs = {frozenset((int(a), int(b))) + for a, b in np.asarray(skin_edges)[band_idx]} + chain_asm, hole_loops = _outcrop_chain_2d(loops_asm, band_pairs) + vS, vE = dm.getDepthStratum(0) pStart, pEnd = dm.getChart() X = _coords(dm)[: vE - vS] @@ -3916,7 +4063,18 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, try: beside_held = np.zeros(len(X), dtype=bool) beside_held[cells[held_c].ravel()] = True - protected = on_wall | held_v | beside_held + deletable_wall = np.zeros(len(X), dtype=bool) + if open_wall is not None: + # An outcrop deletes wall vertices — in the open wall's + # line only. A vertex in a second wall line is a domain + # corner and stays protected: deleting it for one wall + # would breach the other. + deletable_wall = on_wall & (X[:, open_wall[0]] + == open_wall[1]) + other = 1 - open_wall[0] + for v2 in (box_lo[other], box_hi[other]): + deletable_wall &= ~(X[:, other] == v2) + protected = (on_wall & ~deletable_wall) | held_v | beside_held victim = (d_skin < reach_v) & ~protected drop = victim[cells].any(axis=1) @@ -3936,11 +4094,28 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, raise ValueError("the thin volume meets no cell of this mesh") ring, drop = _ring_growing(cells, drop, held_c) - if on_wall[np.asarray(ring)].any(): - raise RuntimeError( - "the ribbon's cavity reached the domain wall; the " - "volume must be interior, with clearance to spare") - if victim[np.asarray(ring)].any(): + removed_wall = [] + if open_wall is None: + if on_wall[np.asarray(ring)].any(): + raise RuntimeError( + "the ribbon's cavity reached the domain wall; the " + "volume must be interior, with clearance to spare") + else: + chain_ids = [len(X) + int(v) for v in chain_asm] + t_ax = 1 - open_wall[0] + chain_t = (float(asm_pts[chain_asm[0]][t_ax]), + float(asm_pts[chain_asm[-1]][t_ax])) + ring, removed_wall = _outcrop_ring_splice( + ring, X, open_wall, chain_ids, chain_t) + off_open = [v for v in ring if v < len(X) + and on_wall[v] + and X[v][open_wall[0]] != open_wall[1]] + if off_open: + raise RuntimeError( + "the ribbon's cavity reached a second domain " + "wall; only a one-wall outcrop is built.") + old_ring = np.asarray([v for v in ring if v < len(X)]) + if victim[old_ring].any(): raise RuntimeError( "a deleted vertex is on the cavity boundary") @@ -3948,14 +4123,35 @@ def _place_thin_volume_2d(dm, polylines, width, label, label_value, if (~drop).any(): referenced[cells[~drop].ravel()] = True orphan = ~referenced & ~victim - if orphan[on_wall].any(): + if (orphan & on_wall & ~deletable_wall).any(): raise RuntimeError( "the cavity would strand a domain-wall vertex; the " "volume must be interior, with clearance to spare") victim |= orphan + # Labels the deleted wall span carried, read before the + # rebuild forgets them — the splice segments and the band are + # relabelled with exactly these after sewing. + wall_pairs = [] + if removed_wall: + edge_pts = [int(dm_work.getFullJoin( + [int(a) + vS, int(b) + vS])[0]) + for a, b in removed_wall] + for i in range(dm_work.getNumLabels()): + name = dm_work.getLabelName(i) + if name in reconnect._TOPOLOGY_LABELS: + continue + lab = dm_work.getLabel(name) + values = lab.getValueIS() + if values is None: + continue + for val in values.getIndices(): + if all(lab.getValue(p) == int(val) + for p in edge_pts): + wall_pairs.append((name, int(val))) + Xall = np.vstack([X, asm_pts]) - holes = [[len(X) + int(v) for v in loop] for loop in loops_asm] + holes = [[len(X) + int(v) for v in loop] for loop in hole_loops] gap_tris, extra = _gmsh_fill_2d(Xall, ring, None, holes=holes) placed = np.vstack([asm_pts, extra]) if len(extra) else asm_pts @@ -3974,8 +4170,18 @@ def mixed(v): raise RuntimeError( "place_thin_volume internal: the gathered region " "touches a shared point; the gather mask under-reached.") + outcrop = None + if open_wall is not None: + # The splice ends: the surviving corner vertices and the + # chain ends they meet — the chain is the ring's tail, in + # the orientation the splice chose. + outcrop = (wall_pairs, + int(removed_wall[0][0]), + int(removed_wall[-1][1]), + int(ring[-len(chain_asm)]) - len(X), + int(ring[-1]) - len(X)) surgery = (np.flatnonzero(victim), np.flatnonzero(drop), made, - placed) + placed, outcrop) except Exception as exc: failure = f"{type(exc).__name__}: {exc}" failures = comm.allgather(failure) @@ -3985,14 +4191,15 @@ def mixed(v): f"place_thin_volume failed on the surgery rank: {real[0]}") if comm.rank == target: - victims_arr, drop_arr, made, placed = surgery + victims_arr, drop_arr, made, placed, outcrop = surgery else: victims_arr = np.empty(0, dtype=np.int64) drop_arr = np.empty(0, dtype=np.int64) made = [] placed = np.empty((0, 2), dtype=float) + outcrop = None - new_dm, _point_map, placed_points = _rebuild_sewn( + new_dm, point_map, placed_points = _rebuild_sewn( dm_work, drop_arr, victims_arr, made, placed) skin_label = label + "_skin" @@ -4028,6 +4235,47 @@ def mixed(v): if real: raise RuntimeError(real[0]) + # An outcropping ribbon's new wall coverage — the two splice segments + # and the band itself — is relabelled EXPLICITLY with what the deleted + # wall span carried, full closures included. + n_wall_local = 0 + if open_wall is not None: + pairs = comm.bcast(outcrop[0] if comm.rank == target else None, + root=target) + for name, val in pairs: + if not new_dm.hasLabel(name): + new_dm.createLabel(name) + if comm.rank == target: + _wp, corner_l, corner_r, a_first, a_last = outcrop + wall_ids = [[int(point_map[corner_l + vS - pStart]), + int(placed_points[a_first])], + [int(placed_points[a_last]), + int(point_map[corner_r + vS - pStart])]] + wall_ids += [[int(placed_points[a]), int(placed_points[b])] + for a, b in (tuple(p) for p in band_pairs)] + for ids in wall_ids: + joined = new_dm.getFullJoin(ids) + if len(joined) != 1: + failure = ("an outcrop wall edge is not an edge of the " + "sewn mesh; the splice or band was not " + "sewn.") + break + for name, val in pairs: + lab = new_dm.getLabel(name) + for q in new_dm.getTransitiveClosure(int(joined[0]))[0]: + lab.setValue(int(q), int(val)) + n_wall_local += 1 + failures = comm.allgather(failure) + real = [f for f in failures if f] + if real: + raise RuntimeError(real[0]) + n_wall = np.array([n_wall_local], dtype=np.int64) + comm.Allreduce(MPI.IN_PLACE, n_wall, op=MPI.SUM) + if int(n_wall[0]) != 2 + len(band_pairs): + raise RuntimeError( + f"{int(n_wall[0])} outcrop wall edges relabelled for " + f"{2 + len(band_pairs)} given.") + counts = np.array([n_zone_local, n_skin_local, len(victims_arr), len(placed)], dtype=np.int64) comm.Allreduce(MPI.IN_PLACE, counts, op=MPI.SUM) @@ -4116,9 +4364,11 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, patches : sequence of array_like In 3-D: one or more PLANAR polygons, ``(N, 3)`` corners each. In 2-D: one or more polylines, ``(N, 2)`` points each, thickened into - ribbons. Interior to the domain with clearance to spare; patches may - cross — that is the point — but must not touch a surface already - embedded. + ribbons. Patches may cross — that is the point — and may run PAST + the domain: the assembly is clipped to the box, and a clipped face + left in a wall plane (an edge in a wall line, in 2-D) becomes the + zone's OUTCROP BAND, carrying both the skin label and the wall's + labels. Patches must not touch a surface already embedded. width : float The layer thickness, a real mesh parameter; ``width < h`` is supported and measured. @@ -4154,8 +4404,8 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, Raises ------ NotImplementedError - In 2-D in parallel (the ribbon shares :func:`place_along_lines`' - serial scope; the 3-D form is the parallel one). + When the zone meets more than one domain wall (box-edge outcrop + bands are not built), or meets one wall in more than one band. RuntimeError, ValueError Carve/fill refusals, always collective. """ diff --git a/tests/parallel/ptest_0855_place_thin_volume_parallel.py b/tests/parallel/ptest_0855_place_thin_volume_parallel.py index c06f946e..bac9f6ec 100644 --- a/tests/parallel/ptest_0855_place_thin_volume_parallel.py +++ b/tests/parallel/ptest_0855_place_thin_volume_parallel.py @@ -82,3 +82,72 @@ def test_refusals_are_collective(): assert all(m is not None for m in messages), ( f"some rank did NOT raise: {[m is None for m in messages]}") assert len(set(messages)) == 1, "ranks raised different errors" + + +def test_a_2d_outcropping_ribbon_embeds_in_parallel(): + """The 2-D outcrop through the same gather-first mechanism, at np>=2. + + The box, the band split and the chain decomposition are collective; + the splice and relabel run on the surgery rank. Asserted: the info + dict identical on every rank, the labels of the advertised sizes, a + band present on the top wall, and no top-wall edge left without Top. + """ + comm = uw.mpi.comm + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.07, + regular=False, qdegree=2) + line = np.array([[0.35, 0.40], [0.60, 1.10]]) # past the top wall + new, info = place_thin_volume(mesh.dm, [line], width=0.04, + label="Zone", label_value=5) + assert info["n_zone_cells"] > 0 + assert _owned_label_count(new, "Zone", 5) == info["n_zone_cells"] + assert _owned_label_count(new, "Zone_skin", 5) == info["n_skin_faces"] + gathered = comm.allgather(info) + assert all(g == gathered[0] for g in gathered) + + pStart, pEnd = new.getChart() + leaves = np.zeros(pEnd - pStart, dtype=bool) + try: + _n, ilocal, _ir = new.getPointSF().getGraph() + if ilocal is not None and len(ilocal): + leaves[np.asarray(ilocal, dtype=np.int64) - pStart] = True + except (ValueError, TypeError): + pass + fS, fE = new.getHeightStratum(1) + vS, vE = new.getDepthStratum(0) + Xn = np.asarray(new.getCoordinatesLocal().array).reshape(-1, 2)[: vE - vS] + top = new.getLabel("Top") + skin = new.getLabel("Zone_skin") + n_band = n_bare = 0 + for f in range(fS, fE): + if leaves[f - pStart] or len(new.getSupport(f)) != 1: + continue + verts = [int(q) - vS for q in new.getTransitiveClosure(f)[0] + if vS <= int(q) < vE] + if all(Xn[v][1] == 1.0 for v in verts): + if top.getValue(f) < 0: + n_bare += 1 + if skin.getValue(f) == 5: + n_band += 1 + n_band = int(comm.allreduce(n_band, op=MPI.SUM)) + n_bare = int(comm.allreduce(n_bare, op=MPI.SUM)) + assert n_band > 0, "the ribbon left no band on the surface" + assert n_bare == 0, "the relabel left top-wall edges without Top" + + +def test_2d_refusals_are_collective(): + """A ribbon stopping short of the wall refuses IDENTICALLY everywhere.""" + comm = uw.mpi.comm + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.07, + regular=False, qdegree=2) + short = np.array([[0.35, 0.40], [0.60, 0.97]]) + message = None + try: + place_thin_volume(mesh.dm, [short], width=0.04, label="Bad") + except (RuntimeError, ValueError) as exc: + message = str(exc) + messages = comm.allgather(message) + assert all(m is not None for m in messages), ( + f"some rank did NOT raise: {[m is None for m in messages]}") + assert len(set(messages)) == 1, "ranks raised different errors" diff --git a/tests/test_0855_place_thin_volume.py b/tests/test_0855_place_thin_volume.py index 2e937c1b..595f9ac0 100644 --- a/tests/test_0855_place_thin_volume.py +++ b/tests/test_0855_place_thin_volume.py @@ -338,3 +338,121 @@ def test_an_outcropping_zone_leaves_a_band_on_the_surface(): err = np.abs(np.asarray(t.data[:, 0]) - (X[:, 0]**2 + X[:, 1]**2 + X[:, 2]**2)) assert float(err.max()) < 1e-8 + + +# ------------------------------------------------------------ 2-D outcrop + +def _top_edge_census_2d(new, skin_value): + """Boundary edges in the top wall line: (total, band, missing Top). + + The band is an edge carrying BOTH the zone's skin label and the wall's; + an edge missing the Top label is a hole the relabel left in the wall. + """ + fS, fE = new.getHeightStratum(1) + vS, vE = new.getDepthStratum(0) + Xn = np.asarray(new.getCoordinatesLocal().array).reshape(-1, 2)[: vE - vS] + top = new.getLabel("Top") + skin = new.getLabel("Zone_skin") + n_top = n_band = n_bare = 0 + for f in range(fS, fE): + if len(new.getSupport(f)) != 1: + continue + verts = [int(q) - vS for q in new.getTransitiveClosure(f)[0] + if vS <= int(q) < vE] + if all(Xn[v][1] == 1.0 for v in verts): + n_top += 1 + if top.getValue(f) < 0: + n_bare += 1 + if skin.getValue(f) == skin_value: + n_band += 1 + return n_top, n_band, n_bare + + +def test_an_outcropping_ribbon_leaves_a_band_on_the_surface(): + """The 2-D zone outcrop: a ribbon specified past the top wall embeds. + + Specify-long: the polyline protrudes, the assembly is clipped in OCC, + and the clipped edge is the BAND — boundary edges carrying both the + skin label and the wall's. The 2-D cap is the two splice segments, so + every top-line boundary edge must still be labelled Top (relabelled, + not stripped) and the domain area is conserved. + + The negative control comes first: the identical census on an INTERIOR + twin counts no band, so the band counted here is produced by the + outcrop, not by a probe that fires on anything. + """ + import sympy + from underworld3.utilities.line_cut import cell_areas + + base = _box2(0.05) + bounds = base._boundaries_with("Zone") + zv = bounds["Zone"].value + + interior, _ = place_thin_volume( + base.dm, [np.array([[0.35, 0.40], [0.60, 0.80]])], width=0.03, + label="Zone", label_value=zv) + n_top0, n_band0, _bare0 = _top_edge_census_2d(interior, zv) + assert n_top0 > 0 and n_band0 == 0, ( + "the census counted a band on an interior ribbon; it cannot " + "validate the outcrop") + + before = float(cell_areas(base.dm).sum()) + line = np.array([[0.35, 0.40], [0.60, 1.10]]) # past the top wall + new, info = place_thin_volume(base.dm, [line], width=0.03, + label="Zone", label_value=zv) + assert info["n_zone_cells"] > 0 + n_top, n_band, n_bare = _top_edge_census_2d(new, zv) + assert n_band > 0, "the ribbon left no band on the surface" + assert n_bare == 0, "the relabel left top-wall edges without Top" + assert float(cell_areas(new).sum()) == pytest.approx(before, rel=1e-12) + + mesh = uw.discretisation.Mesh( + new, simplex=True, qdegree=3, boundaries=bounds, + coordinate_system_type=base.CoordinateSystem.coordinate_type) + x, y = mesh.X + exact = x**2 + y**2 + t = uw.discretisation.MeshVariable("T_band2", mesh, 1, degree=2) + poisson = uw.systems.Poisson(mesh, u_Field=t) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.f = -4.0 + for wall in ("Bottom", "Top", "Left", "Right"): + poisson.add_dirichlet_bc(sympy.Matrix([exact]), wall) + poisson.tolerance = 1e-11 + poisson.solve() + X = np.asarray(t.coords) + err = np.abs(np.asarray(t.data[:, 0]) - (X[:, 0]**2 + X[:, 1]**2)) + assert float(err.max()) < 1e-8, ( + f"the outcropped mesh assembles a wrong operator: max |u - exact| " + f"= {float(err.max()):.3e}") + + +def test_a_ribbon_stopping_short_of_the_wall_still_refuses(): + """No band, no outcrop: the interior contract keeps its refusal. + + A ribbon ending just inside the wall has no clipped edge, so the carve + must still refuse when its cavity reaches the wall — the outcrop path + must not have widened what interior ribbons may do. + """ + mesh = _box2(0.05) + short = np.array([[0.35, 0.40], [0.60, 0.98]]) + with pytest.raises(RuntimeError, match="domain wall"): + place_thin_volume(mesh.dm, [short], width=0.03, label="Zone") + + +def test_an_outcrop_through_two_walls_is_refused(): + mesh = _box2(0.05) + out_top = np.array([[0.30, 0.50], [0.30, 1.10]]) + out_right = np.array([[0.70, 0.50], [1.10, 0.50]]) + with pytest.raises(NotImplementedError, match="more than one"): + place_thin_volume(mesh.dm, [out_top, out_right], width=0.03, + label="Zone") + + +def test_a_ribbon_meeting_the_wall_in_two_bands_is_refused(): + """An arch out the top twice: one skin loop, two bands — refused.""" + mesh = _box2(0.05) + arch = np.array([[0.30, 1.10], [0.35, 0.50], [0.65, 0.50], + [0.70, 1.10]]) + with pytest.raises(NotImplementedError, match="more than one band"): + place_thin_volume(mesh.dm, [arch], width=0.03, label="Zone")