From 9d16872ac3951cf0e062e96a9094999db8f07142 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 16 Aug 2026 21:44:00 +1000 Subject: [PATCH 01/10] The general 3-D outcrop cap: per-region collar, crease overlay, smooth walls MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The zone outcrop no longer needs an axis-aligned wall. The cap over the bowl is re-triangulated PER COPLANAR BOUNDARY REGION, each piece meshed flat in its region's own plane, so the cap stays exactly on the faceted surface and volume conservation survives. Where the band crosses a crease the two sides segment the line differently (mesh vertices against assembly nodes), so a 1-D overlay merges both node sets by arc coordinate and keeps each elementary sub-segment for the side(s) whose bowl covers it and the band does not; a one-side segment must be a whole mesh edge, because a cavity-shell face matches it. The frame rule (_outcrop_frame_3d) separates the two notions the third dimension forces apart. Where the cavity may OPEN is by SMOOTH WALL — coplanar regions grouped across low-dihedral creases — so a faceted sphere's bowl legitimately spills onto facets beside the band's own while a box's other walls stay refused. What the carve may DELETE is by shape: band-covered vertices, single-region interiors, and interior vertices of straight creases between touched regions; a vertex where three or more regions meet is the domain's shape and is protected. A protected vertex whose whole cell star drops is not stranded: it is a collar node, and the gap fill's tets reference it — that is how a curved boundary keeps its faceting through the surgery. Wall labels are restored per removed face (the 2-D per-segment rule one level up): each new wall triangle takes the labels of the removed face it lies on, so a bowl spanning several walls restores each wall's own. The Euler gate now demands CONSERVATION of the input's Euler number rather than 1 — a spherical shell is S^2 x I, Euler 2, and refused before for its topology, not for any defect. The same assumption in place_sheet and remove_embedded is marked TODO(BUG). _trace_wall_code and its refusals are deleted; the box path runs through the general machinery and the box oracle tests are unchanged. Driven on a rotated box (no wall axis-aligned), a band across the Top/Front box edge (trace on both walls), and a spherical shell outer surface (75 trace facets, 6 vertices removed): volume conserved to 1e-12 in each. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/place_surface.py | 765 +++++++++++++++++---- 1 file changed, 624 insertions(+), 141 deletions(-) diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index e285dddb..d37718d6 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -1384,6 +1384,38 @@ def _sheet_distance(X, pts, tris): return best +def _nearest_facet(X, pts, tris): + """``(distance, facet index)`` from each point of ``X`` to a + triangulated sheet — :func:`_sheet_distance` keeping the argmin, for + when the answer is WHICH facet a point lies on, not how far it is.""" + best = np.full(len(X), np.inf) + which = np.zeros(len(X), dtype=np.int64) + for k, t in enumerate(tris): + A, B, C = pts[t[0]], pts[t[1]], pts[t[2]] + ab, ac = B - A, C - A + n = np.cross(ab, ac) + nn = float(n @ n) + rel = X - A + d00, d01, d11 = float(ab @ ab), float(ab @ ac), float(ac @ ac) + d20, d21 = rel @ ab, rel @ ac + denom = d00 * d11 - d01 * d01 + v = (d11 * d20 - d01 * d21) / denom + w = (d00 * d21 - d01 * d20) / denom + inside = (v >= 0.0) & (w >= 0.0) & (v + w <= 1.0) + d_plane = np.abs(rel @ n) / np.sqrt(nn) + d_edge = np.full(len(X), np.inf) + for P, Q in ((A, B), (B, C), (C, A)): + e = Q - P + u = np.clip(((X - P) @ e) / float(e @ e), 0.0, 1.0) + d_edge = np.minimum( + d_edge, np.linalg.norm(X - (P + u[:, None] * e), axis=1)) + d = np.where(inside, d_plane, d_edge) + closer = d < best + best[closer] = d[closer] + which[closer] = k + return best, which + + def _propagate_vertex(dm, chart_values, mpi_op, np_combine): """Reconcile a chart-length per-point array over the point star-forest. @@ -1612,7 +1644,7 @@ def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, shell, cap_faces, drop = _closed_shell_3d( dm, X, cells, drop, victim, held_cells, shared_chart, "sheet", - open_wall=open_wall) + open_vertex=on_open if open_wall is not None else None) # The straddle rule — or the shell growth — can drop the whole star of # a vertex that is NOT a victim; such a vertex is on no shell face and @@ -1750,7 +1782,7 @@ def _outcrop_chain(pts, tris, on_plane): def _closed_shell_3d(dm, X, cells, drop, victim, held_cells, shared_chart, - noun, open_wall=None): + noun, open_vertex=None): """The cavity shell over ``drop``, GROWN at pinch edges until manifold. Shared by every 3-D carve. The union of victim stars around any object @@ -1761,16 +1793,18 @@ def _closed_shell_3d(dm, X, cells, drop, victim, held_cells, shared_chart, original mesh did not). Growing the drop at every non-manifold edge merges the wedges; dropping more cells only enlarges the fill. - ``open_wall = (axis, value)`` lets the cavity OPEN onto one flat wall — - the outcrop bowl. A dropped cell's wall face lying in that plane becomes - a CAP face (returned separately; the caller pre-meshes the cap); any - other wall contact still refuses. The manifold check runs on shell - ∪ cap, which together must close. + ``open_vertex`` (a vertex mask) lets the cavity OPEN onto the boundary + near an outcrop — the bowl. A dropped cell's wall face whose vertices + all carry the mask becomes a CAP face (returned separately; the caller + pre-meshes the cap); any other wall contact still refuses. The mask is + the caller's frame rule: the wall plane for a box-framed sheet, the + band-touched coplanar regions for a general zone. The manifold check + runs on shell ∪ cap, which together must close. Refusals stay per-object worded via ``noun``; a wall or seam contact and a growth that would need a held cell are refusals, not growth. Returns ``(shell, cap_faces, drop)`` with ``drop`` possibly grown; ``cap_faces`` - is empty when ``open_wall`` is None. + is empty when ``open_vertex`` is None. """ from collections import Counter @@ -1805,8 +1839,8 @@ def _closed_shell_3d(dm, X, cells, drop, victim, held_cells, shared_chart, "after the gather; the region marking under-reached. " "A defect, not a configuration error.") verts = face_verts[f] - if open_wall is not None and all( - X[v][open_wall[0]] == open_wall[1] for v in verts): + if open_vertex is not None and all( + open_vertex[v] for v in verts): cap_faces.append((f, verts)) continue raise RuntimeError( @@ -2628,6 +2662,9 @@ def place_sheet(dm, points, triangles, label=CUT_LABEL, label_value=1, f"the placement changed the domain volume: " f"{volume_before[0]:.12f} -> {volume_after[0]:.12f}") + # TODO(BUG): Euler 1 assumes a ball-topology domain; a spherical + # shell is Euler 2 and refuses here. place_thin_volume gates on + # CONSERVATION of the input's Euler number instead — do the same. owned = np.asarray(_owned_stratum_counts(new), dtype=np.int64) comm.Allreduce(MPI.IN_PLACE, owned, op=MPI.SUM) nv_g, ne_g, nf_g, nc_g = (int(x) for x in owned) @@ -2800,21 +2837,15 @@ def _snap_to_boundary_2d(xy, loops, tol=1e-9): return xy -def _occ_domain_3d(occ, verts, tris): - """The domain as ONE OCC solid built from its boundary triangles. +def _coplanar_regions(verts, tris): + """The boundary's COPLANAR REGIONS: adjacent coplanar facets merged. - Adjacent coplanar triangles merge into single planar faces first — a - box wall becomes one rectangle, so a box builds the same six-faced tool - the analytic box gave, and the boolean does not imprint the wall's mesh - spacing onto the clipped faces — and each merged face's rim compresses - to its corners. The chain shared by two merged faces lies in both - planes, hence on a straight line, so both faces compress it to the same - endpoints; a vertex where three or more faces meet is always kept, or - two faces whose planes cross the same line would disagree about it. - Faces share OCC points and lines, so the shells they close into are - topologically sewn. Returns ``(volume_tag, planes)`` with ``planes`` a - list of ``(anchor, unit_normal)`` per merged face, for the snap; the - caller synchronizes. + Two triangles sharing an edge join when their planes agree to rounding, + so a box wall is one region while a faceted sphere keeps every triangle + its own — the region is the 3-D analogue of the compressed straight + segment of a 2-D boundary loop. Region ids are member triangle indices, + arbitrary but deterministic. Returns ``(region, planes)``: the region id + per triangle and ``{region: (anchor, unit_normal)}``. """ tris = np.asarray(tris, dtype=np.int64) T = verts[tris] @@ -2850,6 +2881,35 @@ def find(a): if coplanar: parent[find(t)] = find(s) + region = np.array([find(t) for t in range(len(tris))], dtype=np.int64) + planes = {int(r): (verts[int(tris[r][0])], unit_n[r]) + for r in np.unique(region)} + return region, planes + + +def _occ_domain_3d(occ, verts, tris): + """The domain as ONE OCC solid built from its boundary triangles. + + Adjacent coplanar triangles merge into single planar faces first — a + box wall becomes one rectangle, so a box builds the same six-faced tool + the analytic box gave, and the boolean does not imprint the wall's mesh + spacing onto the clipped faces — and each merged face's rim compresses + to its corners. The chain shared by two merged faces lies in both + planes, hence on a straight line, so both faces compress it to the same + endpoints; a vertex where three or more faces meet is always kept, or + two faces whose planes cross the same line would disagree about it. + Faces share OCC points and lines, so the shells they close into are + topologically sewn. Returns ``(volume_tag, planes)`` with ``planes`` a + list of ``(anchor, unit_normal)`` per merged face, for the snap; the + caller synchronizes. + """ + tris = np.asarray(tris, dtype=np.int64) + T = verts[tris] + region, planes_of = _coplanar_regions(verts, tris) + + def find(t): + return int(region[t]) + regions = {} for t in range(len(tris)): regions.setdefault(find(t), []).append(t) @@ -2900,8 +2960,7 @@ def compress_rim(loop): rim[e] += 1 rim_edges = [tuple(e) for e, k in rim.items() if k == 1] loops = _skin_loops(rim_edges, what="a domain face's rim") - m = unit_n[members[0]] - anchor = verts[int(tris[members[0]][0])] + anchor, m = planes_of[r] planes.append((anchor, m)) # In-plane coordinates for the outer-vs-hole ranking only. u = T[members[0], 1] - T[members[0], 0] @@ -3238,32 +3297,414 @@ def _outcrop_frame_2d(comp_loops, asm_pts, band_pairs, X, on_wall): return deletable, near_touched -def _trace_wall_code(skin_xyz, skin_facets, trace_idx, box_lo, box_hi): - """The single axis-aligned wall the trace lies in, as the legacy - ``2*axis + side`` code — the frame the carve and the sew still run - against. One wall only; more refuses (box-edge bands are not built), - and a trace on a non-axis-aligned boundary refuses until the general - carve/sew exists. +# Adjacent boundary facets whose normals differ by less than this are the +# same SMOOTH WALL: the faceting of a curved boundary (angle -> 0 under +# refinement) groups into one wall, while a box-like corner (90 degrees) +# separates two. The outcrop bowl may open anywhere on a band-touched wall. +_WALL_DIHEDRAL_COS = np.cos(np.deg2rad(45.0)) + + +def _boundary_walls(dom_verts, dom_tris, region): + """Group coplanar regions into SMOOTH WALLS across low-angle creases. + + Returns ``{region: wall}`` with wall ids arbitrary but deterministic. """ - codes = np.full(len(trace_idx), -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_facets[trace_idx]][:, :, axis] - == value).all(axis=1) - codes[on] = 2 * axis + side - if (codes < 0).any(): - # TODO(DESIGN): the general carve/sew (arc-ordered splice, cap - # meshed on the boundary facets) lifts this; the trace split above - # already finds the band on any boundary. - raise NotImplementedError( - "the zone outcrops on a non-axis-aligned boundary; the " - "general carve/sew is not yet built.") - walls = set(int(c) for c in codes) - if len(walls) > 1: - raise NotImplementedError( - "the zone meets more than one domain wall; box-edge outcrop " - "bands are not built.") - return walls.pop() + tris = np.asarray(dom_tris, dtype=np.int64) + T = dom_verts[tris] + n = np.cross(T[:, 1] - T[:, 0], T[:, 2] - T[:, 0]) + n = n / np.linalg.norm(n, axis=1)[:, None] + + parent = {int(r): int(r) for r in np.unique(region)} + + def find(a): + while parent[a] != a: + parent[a] = parent[parent[a]] + a = parent[a] + return a + + edge_first = {} + for t, tri in enumerate(tris): + a, b, c = sorted(int(v) for v in tri) + for e in ((a, b), (a, c), (b, c)): + if e not in edge_first: + edge_first[e] = t + continue + s = edge_first[e] + # The gathered complex carries no consistent orientation, so + # the dihedral test must not depend on the facet normals' + # signs; |cos| also merges a near-zero fold, which a domain + # boundary does not have. + if abs(float(n[t] @ n[s])) > _WALL_DIHEDRAL_COS: + parent[find(int(region[t]))] = find(int(region[s])) + return {r: find(r) for r in parent} + + +def _outcrop_frame_3d(X, on_wall, dom_verts, dom_tris, region, + skin_xyz, band_tris, tol=1e-9): + """Which mesh boundary vertices a 3-D outcrop may DELETE, and which + count as boundary contact NEAR the outcrop. + + The 2-D frame rule one dimension up, with the two notions the third + dimension separates. NEAR — where the cavity may open into cap faces — + is by SMOOTH WALL (:func:`_boundary_walls`): on a faceted sphere the + bowl legitimately spills onto facets beside the band's own, while a + box's other walls stay refused across their sharp creases. DELETABLE + is by SHAPE: a vertex may go when the band re-provides its geometry + (it lies ON the band), when it is interior to a single band-touched + coplanar region (the collar re-mesh stays in the region's own plane), + or when it is interior to a STRAIGHT crease between two touched + regions — a straight crease's interior vertices carry no shape, and + the overlay rebuilds the line through the survivors. A vertex where + three or more regions meet (a box corner, every vertex of a faceted + sphere) is the domain's shape itself — the 2-D compressed corner — + and is never deleted off the band. + + Membership is exact: the gathered complex's vertex coordinates are the + mesh's own bytes, so a mesh vertex's regions are read off the facets + incident to its coordinate, not off a distance. + + Returns ``(deletable, near_touched, touched)`` — masks over the mesh + vertices, and the set of touched region ids. + """ + band_cen = skin_xyz[np.asarray(band_tris)].mean(axis=1) + d_cen, at = _nearest_facet(band_cen, dom_verts, dom_tris) + if (d_cen > tol).any(): + raise RuntimeError( + "the outcrop band lies off the domain boundary; the clip and " + "the boundary disagree") + touched = {int(region[k]) for k in at} + wall_of = _boundary_walls(dom_verts, dom_tris, region) + touched_walls = {wall_of[r] for r in touched} + + row_of = {tuple(v): i for i, v in enumerate(dom_verts)} + regions_at = {} + crease_dirs = {} + edge_first = {} + for t, tri in enumerate(dom_tris): + for v in tri: + regions_at.setdefault(int(v), set()).add(int(region[t])) + a, b, c = sorted(int(v) for v in tri) + for e in ((a, b), (a, c), (b, c)): + if e not in edge_first: + edge_first[e] = t + continue + if int(region[t]) != int(region[edge_first[e]]): + d = dom_verts[e[1]] - dom_verts[e[0]] + d = d / np.linalg.norm(d) + for v in e: + crease_dirs.setdefault(int(v), []).append(d) + + d_band = _sheet_distance(X, skin_xyz, np.asarray(band_tris)) + deletable = np.zeros(len(X), dtype=bool) + near_touched = np.zeros(len(X), dtype=bool) + for i in np.flatnonzero(on_wall): + dv = row_of.get(tuple(X[i])) + if dv is None: + raise RuntimeError( + "a domain-wall vertex is missing from the gathered " + "boundary complex; the wall mask and the complex disagree") + mine = regions_at[dv] + if not any(wall_of[r] in touched_walls for r in mine): + continue + near_touched[i] = True + if d_band[i] < tol: + deletable[i] = True + elif len(mine) == 1: + deletable[i] = mine <= touched + elif len(mine) == 2 and mine <= touched: + dirs = crease_dirs.get(dv, []) + deletable[i] = (len(dirs) == 2 and + float(np.linalg.norm( + np.cross(dirs[0], dirs[1]))) < 1e-9) + return deletable, near_touched, touched + + +def _outcrop_collar_3d(X, alive, cap_faces, cap_region, planes_of, + skin_xyz, band_tris, band_tri_region, band_outline, + tol=1e-9): + """Re-triangulate the outcrop bowl's wall collar, region by region. + + The collar is the bowl's wall footprint minus the band — in the box + case a flat annulus between the cavity's wall rim and the band's + outline. On a general boundary the footprint spans several COPLANAR + REGIONS, so each region's piece is meshed flat in its own plane — the + cap stays exactly on the faceted surface, which is what conserves the + domain volume — and the pieces conform along the CREASES between + regions. + + Along a crease the two sides segment the line differently: the bowl's + nodes are mesh vertices, the band's crossing nodes are assembly nodes, + so edge identity cannot match them. The 1-D overlay merges both node + sets by arc coordinate along the (straight) crease and keeps each + elementary sub-segment for the side(s) whose bowl piece covers it and + the band does not. Deleted crease vertices (``alive`` False — the + frame rule lets a straight crease's interior go) contribute geometry + but not nodes: the chain is rebuilt through the survivors. A + sub-segment covered from ONE side only must be a whole mesh edge, + because the face matching it in the fill's closed surface is a + cavity-shell face carrying the mesh segmentation — anything else there + is a refusal, not a crack. + + ``cap_faces`` are the bowl's wall triangles as rows of ``X`` with + ``alive`` the surviving-vertex mask, ``cap_region`` their coplanar + region ids, ``planes_of`` the region planes; ``band_tris`` are the + band's skin triangles (rows of ``skin_xyz``) with ``band_tri_region`` + their region ids and ``band_outline`` the band's closed outline loop + of skin nodes. + + Returns ``(mesh_nodes, skin_nodes, extra_xyz, tris)`` — the cap in the + fill's three node namespaces (kept mesh vertices, assembly skin nodes, + new in-plane points), ``tris`` indexing their concatenation in that + order, exactly as :func:`_gmsh_fill_annulus_3d` reads its payload. + """ + cap_faces = np.asarray(cap_faces, dtype=np.int64) + cap_region = np.asarray(cap_region, dtype=np.int64) + band_tris = np.asarray(band_tris, dtype=np.int64) + + # The bowl's wall edges. Within one region an edge seen once is the + # bowl rim (matched by a cavity-shell face, so it stays a mesh edge + # verbatim); an edge shared by two regions' faces is a crease edge and + # goes to the overlay; an edge seen twice within one region is + # interior to that region's piece. + edge_faces = {} + for k, tri in enumerate(cap_faces): + a, b, c = sorted(int(v) for v in tri) + for e in ((a, b), (a, c), (b, c)): + edge_faces.setdefault(e, []).append(k) + + soup = {} + crease_edges = {} + for e, ks in edge_faces.items(): + if len(ks) > 2: + raise RuntimeError( + "an outcrop bowl wall edge belongs to more than two wall " + "faces; the boundary complex is not manifold") + rs = sorted({int(cap_region[k]) for k in ks}) + if len(ks) == 1: + soup.setdefault(rs[0], []).append((('m', e[0]), ('m', e[1]))) + elif len(rs) == 2: + crease_edges.setdefault(tuple(rs), []).append(e) + + # Region soups of the coverage tests, once. + region_faces = {int(r): cap_faces[cap_region == r] + for r in np.unique(cap_region)} + + def covered(r, P): + return bool(_sheet_distance(P[None, :], X, region_faces[r])[0] + < tol) + + def in_band(P): + return bool(_sheet_distance(P[None, :], skin_xyz, band_tris)[0] + < tol) + + n_out = len(band_outline) + outline_edges = [(int(band_outline[i]), + int(band_outline[(i + 1) % n_out])) + for i in range(n_out)] + + # ------------------------------------------------- the crease overlay + for pair, edges in crease_edges.items(): + parent = {} + + def find(v): + while parent[v] != v: + parent[v] = parent[parent[v]] + v = parent[v] + return v + + for a, b in edges: + parent.setdefault(a, a) + parent.setdefault(b, b) + parent[find(a)] = find(b) + chains = {} + for a, b in edges: + chains.setdefault(find(a), []).append((a, b)) + + for chain in chains.values(): + chain_all = sorted({v for e in chain for v in e}) + nodes_m = [v for v in chain_all if alive[v]] + origin = X[chain_all[0]] + far = max(chain_all, key=lambda v: float( + np.linalg.norm(X[v] - origin))) + u = X[far] - origin + u = u / np.linalg.norm(u) + + def arc(P): + return float((P - origin) @ u) + + def off_line(P): + return float(np.linalg.norm((P - origin) + - arc(P) * u)) + + span = [min(arc(X[v]) for v in chain_all) - tol, + max(arc(X[v]) for v in chain_all) + tol] + resolved = [('m', v, arc(X[v])) for v in nodes_m] + for s in {v for e in outline_edges for v in e}: + P = skin_xyz[s] + if off_line(P) < tol and span[0] < arc(P) < span[1]: + resolved.append(('a', int(s), arc(P))) + for a, b in outline_edges: + mid = 0.5 * (skin_xyz[a] + skin_xyz[b]) + if (off_line(skin_xyz[a]) < tol + and off_line(skin_xyz[b]) < tol + and off_line(mid) < tol + and span[0] < arc(mid) < span[1]): + raise NotImplementedError( + "the outcrop band runs ALONG a domain crease; a " + "band bounded by a crease is not built.") + for kind, i, s in resolved: + if kind != 'a': + continue + if any(k2 == 'm' and abs(s2 - s) < tol + for k2, _i2, s2 in resolved): + raise RuntimeError( + "a band outline node lands within rounding of a " + "kept mesh vertex on a domain crease; the collapse " + "of boundary imprints is not built in 3-D — move " + "the patch or change the mesh spacing.") + chain_set = {(a, b) if a < b else (b, a) for a, b in chain} + resolved.sort(key=lambda n: n[2]) + for (k0, i0, s0), (k1, i1, s1) in zip(resolved, resolved[1:]): + mid = origin + (0.5 * (s0 + s1)) * u + if in_band(mid): + continue + sides = [r for r in pair if covered(r, mid)] + if not sides: + continue + if len(sides) == 1 and ( + 'a' in (k0, k1) + or ((i0, i1) if i0 < i1 else (i1, i0)) + not in chain_set): + raise RuntimeError( + "the outcrop band splits a wall edge that the " + "cavity shell keeps whole; raise `clearance` so " + "the carve covers the crease on both sides") + for r in sides: + soup.setdefault(r, []).append(((k0, i0), (k1, i1))) + + # The band outline bounds the collar from inside; each outline edge + # belongs to exactly one band triangle, whose region takes it. + band_edge_owner = {} + for k, tri in enumerate(band_tris): + a, b, c = sorted(int(v) for v in tri) + for e in ((a, b), (a, c), (b, c)): + band_edge_owner.setdefault(e, []).append(k) + for a, b in outline_edges: + e = (a, b) if a < b else (b, a) + owners = band_edge_owner.get(e, []) + if len(owners) != 1: + raise RuntimeError( + "a band outline edge does not bound exactly one band " + "facet; the outline and the band disagree") + r = int(band_tri_region[owners[0]]) + soup.setdefault(r, []).append((('a', a), ('a', b))) + + # ------------------------------------- mesh each region's piece flat + mesh_nodes, skin_nodes = [], [] + m_row, a_row = {}, {} + + def namespace(kind, i): + if kind == 'm': + if i not in m_row: + m_row[i] = len(mesh_nodes) + mesh_nodes.append(int(i)) + return ('m', m_row[i]) + if i not in a_row: + a_row[i] = len(skin_nodes) + skin_nodes.append(int(i)) + return ('a', a_row[i]) + + extra_xyz = [] + tris_out = [] + for r, segments in sorted(soup.items()): + local_of, local_nodes = {}, [] + for e in segments: + for kind, i in e: + if (kind, i) not in local_of: + local_of[(kind, i)] = len(local_nodes) + local_nodes.append((kind, i)) + coords = np.array([X[i] if kind == 'm' else skin_xyz[i] + for kind, i in local_nodes]) + m_coords = np.array([X[i] for kind, i in local_nodes + if kind == 'm']) + if any(kind == 'm' and not alive[i] for kind, i in local_nodes): + raise RuntimeError( + "an outcrop collar node was deleted by the carve; the " + "frame rule and the collar disagree") + for kind, i in local_nodes: + if kind != 'a' or not len(m_coords): + continue + if np.linalg.norm(m_coords - skin_xyz[i], axis=1).min() < tol: + raise RuntimeError( + "a band outline node lands within rounding of a kept " + "mesh vertex; the collapse of boundary imprints is " + "not built in 3-D — move the patch or change the mesh " + "spacing.") + + anchor, nrm = planes_of[int(r)] + f0 = region_faces[int(r)][0] + e0 = X[f0[1]] - X[f0[0]] + e0 = e0 / np.linalg.norm(e0) + w0 = np.cross(nrm, e0) + P2 = np.column_stack([(coords - anchor) @ e0, + (coords - anchor) @ w0]) + + loops = _skin_loops( + [(local_of[e[0]], local_of[e[1]]) for e in segments], + what="an outcrop collar piece's boundary") + + def signed_area(loop): + Q = P2[np.asarray(loop)] + return 0.5 * float(Q[:, 0] @ np.roll(Q[:, 1], -1) + - Q[:, 1] @ np.roll(Q[:, 0], -1)) + + containers = {k: [j for j, other in enumerate(loops) + if j != k and _inside_polygon( + P2[np.asarray(other)], P2[loop[0]])] + for k, loop in enumerate(loops)} + for k, loop in enumerate(loops): + if len(containers[k]) % 2: + continue # a hole; its outer takes it + holes = [j for j in range(len(loops)) + if containers[j] + and max(containers[j], + key=lambda c: len(containers[c])) == k + and len(containers[j]) == len(containers[k]) + 1] + ring = loop if signed_area(loop) > 0.0 else loop[::-1] + fill_tris, extra2 = _gmsh_fill_2d( + P2, list(ring), None, + holes=[list(loops[j]) for j in holes]) + base = len(local_nodes) + extra_row = {} + for t in fill_tris: + row = [] + for v in t: + if v < base: + row.append(namespace(*local_nodes[int(v)])) + else: + j = int(v) - base + if j not in extra_row: + x2 = extra2[j] + extra_xyz.append(anchor + x2[0] * e0 + + x2[1] * w0) + extra_row[j] = len(extra_xyz) - 1 + row.append(('x', extra_row[j])) + tris_out.append(row) + + n_m, n_a = len(mesh_nodes), len(skin_nodes) + + def flat(kind, i): + if kind == 'm': + return i + if kind == 'a': + return n_m + i + return n_m + n_a + i + + tris = np.array([[flat(*v) for v in row] for row in tris_out], + dtype=np.int64) + extra = (np.asarray(extra_xyz, dtype=float) if extra_xyz + else np.zeros((0, 3), dtype=float)) + return mesh_nodes, skin_nodes, extra, tris def _single_loop(edges, what): @@ -3290,7 +3731,8 @@ def _single_loop(edges, what): def _carve_around_volume_3d(dm, X, cells, skin_pts, skin_tris, reach_vertex, reach_cell, held_cells, on_wall, shared_chart, - seed_drop=None, open_wall=None): + seed_drop=None, open_deletable=None, + open_near=None): """Victims, dropped tets and the closed shell around a FAT object. Differs from the sheet's carve in two measured ways. The reach is a @@ -3301,6 +3743,11 @@ def _carve_around_volume_3d(dm, X, cells, skin_pts, skin_tris, reach_vertex, the drop set is GROWN at every non-manifold shell edge until the shell closes; dropping more cells only enlarges the fill (thin_volume_spike: converges in a few rounds). + + An OUTCROP passes its frame rule as two vertex masks + (:func:`_outcrop_frame_3d`): ``open_deletable`` — wall vertices the + carve may take — and ``open_near`` — wall vertices near the outcrop, + where the cavity may open into cap faces. """ d_skin = _sheet_distance(X, skin_pts, skin_tris) @@ -3308,9 +3755,8 @@ def _carve_around_volume_3d(dm, X, cells, skin_pts, skin_tris, reach_vertex, if held_cells: for c in held_cells: held_vertex[cells[c]] = True - on_open = np.zeros(len(X), dtype=bool) - if open_wall is not None: - on_open = X[:, open_wall[0]] == open_wall[1] + on_open = (open_deletable if open_deletable is not None + else np.zeros(len(X), dtype=bool)) victim = ((d_skin < reach_vertex) & (~on_wall | on_open) & ~held_vertex) @@ -3339,19 +3785,26 @@ def _carve_around_volume_3d(dm, X, cells, skin_pts, skin_tris, reach_vertex, shell, cap_faces, drop = _closed_shell_3d(dm, X, cells, drop, victim, held_cells, shared_chart, "thin volume", - open_wall=open_wall) - - # The growth can swallow the whole star of a vertex that is NOT itself a - # victim. Such a vertex is on no shell face — a shell face keeps a - # surviving cell — so it would come through the rebuild as an ISOLATED - # point and the global Euler gate reads 2, not 1 (caught by CI: the - # growth pattern follows the assembly mesh and is gmsh-version- - # dependent). Every surviving vertex must have a surviving cell. + open_vertex=open_near) + + # The carve can swallow the whole star of a vertex that is NOT itself + # a victim. An unreferenced survivor would ride through the rebuild as + # an isolated point (global Euler 2, not 1 — the class CI caught), so + # an interior one is promoted to a victim. A PROTECTED wall vertex + # near the outcrop needs no kept cell at all: it is a COLLAR node — + # the cap is re-triangulated through the surviving wall vertices, + # which is what preserves a faceted boundary's shape, and the gap + # fill's tets reference it. On a curved wall this is the common case: + # every facet vertex near the band is protected, and each cell of its + # small star holds some interior victim. referenced = np.zeros(len(X), dtype=bool) if (~drop).any(): referenced[cells[~drop].ravel()] = True orphan = ~referenced & ~victim - if orphan[on_wall & ~on_open].any(): + if open_near is not None: + collar_kept = orphan & on_wall & ~on_open & open_near + orphan &= ~collar_kept + if (orphan & on_wall & ~on_open).any(): raise RuntimeError( "the cavity would strand a domain-wall vertex; the volume must " "be interior, with clearance to spare") @@ -3369,11 +3822,11 @@ def _gmsh_fill_annulus_3d(shell_xyz, shell_tris, skin_xyz, skin_tris, OUTCROPPING zone (``cap`` given): the gap's boundary is ONE closed surface of three discrete pieces — the shell, the pre-meshed CAP over - the bowl (an annulus in the wall plane whose outer ring is the shell's - wall rim and whose HOLE is the zone's band outline), and the INTERIOR - part of the skin, which is open and rim-matched to the cap's hole. The - band itself is not part of the gap's boundary — the zone touches the - wall directly there. + the bowl (the wall collar :func:`_outcrop_collar_3d` builds — in the + box case an annulus between the shell's wall rim and the zone's band + outline), and the INTERIOR part of the skin, which is open and + rim-matched to the cap's hole. The band itself is not part of the + gap's boundary — the zone touches the wall directly there. Returns ``(points, tets, moved, skin_out, n_shell, cap_out)`` with the fill's nodes ordered shell first, skin second, cap extras third, new @@ -4370,6 +4823,9 @@ def remove_embedded(dm, label, label_value=1, clearance=0.6, verbose=False): f"the removal changed the domain volume: " f"{volume_before[0]:.12f} -> {volume_after[0]:.12f}") + # TODO(BUG): Euler 1 assumes a ball-topology domain; a spherical + # shell is Euler 2 and refuses here. place_thin_volume gates on + # CONSERVATION of the input's Euler number instead — do the same. owned = np.asarray(_owned_stratum_counts(new), dtype=np.int64) comm.Allreduce(MPI.IN_PLACE, owned, op=MPI.SUM) nv_g, ne_g, nf_g, nc_g = (int(x) for x in owned) @@ -4856,10 +5312,12 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, 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. 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. + the domain: the assembly is clipped against the mesh's own + boundary, and a clipped face left ON the boundary (an edge, in + 2-D) becomes the zone's OUTCROP BAND, carrying the skin label, the + wall's labels and the ``