diff --git a/docs/advanced/fault-networks.md b/docs/advanced/fault-networks.md index 9d33a4c70..15609fa07 100644 --- a/docs/advanced/fault-networks.md +++ b/docs/advanced/fault-networks.md @@ -256,17 +256,36 @@ directly. Place is currently OPT-IN for networks: on graded solve is pathological — an open operator-health work item; on uniform bases it is healthy. +`build(width=...)` gives the 3-D network the same finite-width +contract as 2-D: the margin-expanded patches are thickened by +`±width/2` into ONE fused band (junctions free), and each un-expanded +patch is embedded in the band as a conforming mid-surface — so the +same mesh is cut and split (`realisation="split"`) or left whole for +the volumetric weak plane (`realisation="ti"`), exactly as in 2-D. +The honoured footprints are exact for planar patches (band cells +within the patch's own in-plane outline and half a width of its +plane), the weak-plane director is the patch normal, and `slips()` +reports the plane form of each gauge: the tangential pair jump for +the split, the in-plane velocity jump across the layer for the weak +plane. Junction glue in 3-D remains `damage_yield`'s tubes about the +junction segments — `junction_cells` is the 2-D ribbon rule and +refuses in 3-D. Interior networks only: a band that reaches the +domain boundary is refused loudly (its embedded mid-surface cannot +yet be clipped against the boundary). + v1 scope, refused loudly outside it: planar patches (the `rim_polygon` contract), convex rims, genuine X crossings (a near-miss — close but not crossing — is refused rather than guessed -at); parallel MULTI-fault splitting (the pairing does not yet migrate -through redistribution — single faults are parallel-validated). +at). Multi-fault networks split and solve in parallel: `split_faults` +redistributes ONCE, keyed on the union of the network's facets, and +every split then runs with serial topology. ## Limitations -- 3-D: planar convex patches, X crossings only, serial (above); the - weak-plane realisation is 2-D for now — place 3-D zones with - `place_thin_volume` directly. +- 3-D: planar convex patches, X crossings only (above). Finite-width + bands are interior-only (an outcropping band is refused), and + `junction_cells` is the 2-D ribbon rule — 3-D junction glue is + `damage_yield`'s tubes. - One damage dial per network in `damage_yield` (per-junction values: build the expression with `uw.meshing.damage_zone_yield` directly). - Time-dependent damage (wear-in/healing) is study-level for now: see diff --git a/src/underworld3/meshing/fault_network.py b/src/underworld3/meshing/fault_network.py index 26c42731e..957440316 100644 --- a/src/underworld3/meshing/fault_network.py +++ b/src/underworld3/meshing/fault_network.py @@ -96,6 +96,51 @@ def _nearest_segment_normals(P, X): return np.column_stack([-T[:, 1], T[:, 0]]) +def _patch_normal(P): + """Unit normal of a planar polygon by Newell's method — robust to + collinear leading vertices, which a clipped rim can carry (the + two-edge cross product is not).""" + P = np.asarray(P, dtype=float) + n = np.cross(P, np.roll(P, -1, axis=0)).sum(axis=0) + norm = float(np.linalg.norm(n)) + if norm == 0.0: + raise ValueError("degenerate patch: the rim spans no plane") + return n / norm + + +def _expand_convex_polygon(P, dist): + """Offset a convex planar polygon outward in its own plane: each edge + moves ``dist`` along its in-plane outward normal, adjacent edge lines + re-intersected — the 3-D tip margin (the band extends past the fault; + the mid-surface is the fault itself, honoured exactly).""" + P = np.asarray(P, dtype=float) + n_hat = _patch_normal(P) + m = len(P) + anchors, dirs = [], [] + for i in range(m): + e = P[(i + 1) % m] - P[i] + out_dir = np.cross(e, n_hat) + out_dir /= np.linalg.norm(out_dir) + anchors.append(P[i] + dist * out_dir) + dirs.append(e / np.linalg.norm(e)) + out = np.empty_like(P) + for i in range(m): + # corner i: intersect edge line i-1 with edge line i (in-plane) + a0, d0 = anchors[i - 1], dirs[i - 1] + a1, d1 = anchors[i], dirs[i] + w = np.cross(d0, d1) + denom = float(w @ w) + if denom < 1e-24: + # collinear adjacent edges (a clipped rim keeps such + # points): both edges offset onto one line, and the corner + # is simply the offset point itself + out[i] = a1 + continue + t = float(np.cross(a1 - a0, d1) @ w) / denom + out[i] = a0 + t * d0 + return out + + def _densify_polyline(E, piece, per_segment=4): """Points along an extended spine, ``per_segment`` per edge, each flagged as lying on a CUT — an edge whose two vertices belong to the @@ -156,6 +201,7 @@ def __init__(self, faults, hierarchy=None): raise ValueError(f"hierarchy names not in the network: " f"{sorted(unknown)}") self.h_near = None + self.ligament = None self.prepared = None self.junctions = None self.report = None @@ -184,6 +230,7 @@ def prepare(self, h, ligament=2.0, through=None, verbose=True): cut anywhere); the hierarchy handles everything else pairwise. """ self.h_near = float(h) + self.ligament = float(ligament) if self.dim == 3: from .fault_network_3d import prepare_fault_surfaces if through: @@ -251,13 +298,20 @@ def build(self, base=None, h_far=None, band=0.03, ramp=0.08, """ if self.prepared is None: raise RuntimeError("call prepare(h=...) first") - meshers = {2: ("network", "ladder"), 3: ("embed", "place")}[self.dim] - if mesher is None: - mesher = meshers[0] - if mesher not in meshers: - raise ValueError( - f"mesher must be one of {meshers} in {self.dim}-D, not " - f"{mesher!r}") + if self.dim == 3 and width is not None: + if mesher not in (None, "network"): + raise ValueError( + "the 3-D band has one mesher, 'network' — the fused " + f"band with embedded mid-surfaces — not {mesher!r}") + else: + meshers = {2: ("network", "ladder"), + 3: ("embed", "place")}[self.dim] + if mesher is None: + mesher = meshers[0] + if mesher not in meshers: + raise ValueError( + f"mesher must be one of {meshers} in {self.dim}-D, " + f"not {mesher!r}") if realisation not in ("split", "ti"): raise ValueError( f"realisation must be 'split' or 'ti', not {realisation!r}") @@ -269,17 +323,17 @@ def build(self, base=None, h_far=None, band=0.03, ramp=0.08, self.realisation = realisation self.width = None if width is None else float(width) if self.dim == 3: - if realisation != "split": - raise NotImplementedError( - "the 3-D network builds the split realisation only; " - "place the patches with place_thin_volume for a " - "volumetric zone.") if width is not None: - raise NotImplementedError( - "the 3-D network does not place a band: its patches " - "are meshed conforming (mesher='embed') or placed as " - "sheets (mesher='place'), both of zero thickness. For " - "a finite-width 3-D zone call place_thin_volume.") + if base is not None: + raise NotImplementedError( + "the 3-D band builds its own base box; adapting " + "a supplied base is not built yet") + return self._build_3d_band( + h_far=h_far, qdegree=qdegree, realisation=realisation, + margin_rings=margin_rings, + carve_clearance=carve_clearance) + # width=None here implies realisation == "split": the + # ti-needs-width check above already refused the other case return self._build_3d(h_far=h_far, qdegree=qdegree, mesher=mesher) from .cartesian import UnstructuredSimplexBox @@ -425,6 +479,142 @@ def _junction_margins(self, rings): out.append((ends[0], ends[1])) return out + def _build_3d_band(self, h_far=None, qdegree=2, realisation="split", + margin_rings=2, carve_clearance=0.3, + minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0)): + """The finite-width 3-D network: ONE band, both realisations. + + The prepared patches, expanded in-plane by the tip margin, are + thickened by ``±width/2`` and fused (junctions free); each + UN-expanded patch is embedded in the band as a conforming + mid-surface, so the same mesh is cut and split + (``realisation="split"``) or left whole for the weak plane + (``"ti"``) — the 2-D contract, one dimension up. Interior + networks only for now (an outcropping band cannot yet carry an + embedded mid-surface), and the MARGIN-EXPANDED band must clear + the walls by about a base cell for the carve to succeed — + shrink ``margin_rings`` or refine ``h_far`` when a patch sits + close to a boundary (the carve refuses loudly).""" + import underworld3 as uw + from enum import Enum + from underworld3.utilities.place_surface import ( + place_thin_volume, _label_mid_surface, _cell_centroids_of) + from underworld3.utilities.custom_mg import adopt_hierarchy + from underworld3.utilities.fault_split import split_faults + from .cartesian import UnstructuredSimplexBox + + h = self.h_near + h_far = 4.0 * h if h_far is None else float(h_far) + if self.junctions and float(margin_rings) > self.ligament - 0.5: + # the junction gap is PHYSICS (the intact ligament); the tip + # margin is convenience. A margin that reaches within half a + # cell of the gap welds the expanded bands at the junction + # line and the split then produces degenerate self-paired + # nodes there (measured: a bit-frozen residual floor). + raise ValueError( + f"margin_rings={float(margin_rings):g} closes the " + f"junction gap (ligament={self.ligament:g}): the " + "expanded bands weld at the junction and the split " + "degenerates. Keep margin_rings <= ligament - 0.5.") + base = UnstructuredSimplexBox( + cellSize=h_far, minCoords=minCoords, maxCoords=maxCoords, + refinement=1, qdegree=qdegree) + margin = float(margin_rings) * h + expanded = [_expand_convex_polygon(P, margin) + for _n, P in self.prepared] + # interior networks only: an outcropping band cannot yet carry + # its embedded mid-surface (the surface would need the same + # boundary clip), so refuse before any meshing happens + lo = np.asarray(minCoords, dtype=float) + hi = np.asarray(maxCoords, dtype=float) + for (name, _P), E in zip(self.prepared, expanded): + n_hat = _patch_normal(E) + slab = np.vstack([E + 0.5 * self.width * n_hat, + E - 0.5 * self.width * n_hat]) + if (slab <= lo).any() or (slab >= hi).any(): + raise NotImplementedError( + f"the band of {name!r} reaches the domain boundary; " + "an outcropping 3-D band cannot yet carry its " + "embedded mid-surface. Keep the network interior, " + "or place the zone with place_thin_volume directly " + "(no embed).") + dm, info = place_thin_volume( + base.dm, expanded, self.width, label="Band", label_value=71, + clearance=carve_clearance, size=h, mesher="network", + embed=[np.asarray(P, dtype=float) for _n, P in self.prepared]) + values = {} + for k, (name, _P) in enumerate(self.prepared): + values[name] = 41 + k + n_faces = _label_mid_surface(dm, info["embedded_nodes"][k], + name, values[name]) + if n_faces == 0: + raise RuntimeError( + f"the embedded mid-surface of {name!r} labelled no " + f"faces in the placed mesh") + members = [(b.name, b.value) for b in base.boundaries] + members += [(n, v) for n, v in values.items()] + mesh = uw.discretisation.Mesh( + dm, simplex=True, qdegree=qdegree, + coordinate_system_type=base.CoordinateSystem.coordinate_type, + boundaries=Enum("boundaries", members), verbose=False) + band = mesh.cells_labelled("Band", 71) + adopt_hierarchy(mesh, base, fac_zone=band) + if realisation == "split": + mesh = split_faults(mesh, [n for n, _P in self.prepared]) + # reduce first, then branch: the defect must raise on every + # rank together or not at all. A healthy pairing is a + # bijection between disjoint sides: a node paired with + # itself OR appearing as both a minus and a plus (a chain) + # is the same degeneracy. + n_bad = sum(len(set(pairs) & set(pairs.values())) + for pairs in mesh._fault_point_pairs.values()) + n_bad = mesh.dm.comm.tompi4py().allreduce(n_bad) + if n_bad: + raise RuntimeError( + f"the network split produced {n_bad} node(s) on " + "both sides of a pairing — the embedded " + "mid-surfaces are degenerate (a defect, not a " + "configuration error).") + mesh._custom_mg_fac_zone = None # a split fault needs no patch + band = mesh.cells_labelled("Band", 71) + # honoured footprints: band cells within the USER patch's own + # in-plane outline and half a width of its plane — planar patches + # make the rule exact; the expanded margin stays unpainted. Near a + # junction a cell can satisfy two patches' rules, so ownership is + # NEAREST PLANE among them (the 2-D nearest-spine ownership, one + # dimension up): footprints are disjoint and the weak-plane + # director is well-defined. + ids, cen = _cell_centroids_of(mesh.dm, band) + plane_dist = np.full((len(self.prepared), len(ids)), np.inf) + for j, (_name, P) in enumerate(self.prepared): + P = np.asarray(P, dtype=float) + n_hat = _patch_normal(P) + d = (cen - P[0]) @ n_hat + inside = np.abs(d) <= 0.5 * self.width + 0.35 * h + in_plane = cen - np.outer(d, n_hat) + for i in range(len(P)): # convex: inside every edge + e = P[(i + 1) % len(P)] - P[i] + out_dir = np.cross(e, n_hat) + out_dir /= np.linalg.norm(out_dir) + inside &= (in_plane - P[i]) @ out_dir <= 0.35 * h + plane_dist[j, inside] = np.abs(d[inside]) + owner = np.argmin(plane_dist, axis=0) + owned = np.isfinite(plane_dist.min(axis=0)) + footprints = {} + for j, (name, _P) in enumerate(self.prepared): + m = np.zeros_like(band) + m[ids[(owner == j) & owned]] = True + footprints[name] = m + self.info = {"n_cells": int(mesh.dm.getHeightStratum(0)[1]), + "band": band, "footprints": footprints, + "spacing": [h] * len(self.prepared), + "width": float(self.width), "mesher": "network", + "margin_rings": [(margin_rings, margin_rings)] + * len(self.prepared)} + self.mesh = mesh + self._make_surfaces() + return self.mesh + def _build_3d(self, h_far=None, qdegree=2, mesher="embed", minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), band=None, ramp=None, max_levels=2, clearance=0.8): @@ -779,7 +969,10 @@ def ti_fields(self, eta_1, eta_0=1.0, tag=""): m_ = foots[name] # there, so TI is isotropic if not m_.any(): continue - dvals[m_] = _nearest_segment_normals(P, cen[m_]) + if dim == 3: # planar patch: ONE normal + dvals[m_] = _patch_normal(P) + else: + dvals[m_] = _nearest_segment_normals(P, cen[m_]) ndir.array[...] = dvals.reshape(ndir.array.shape) self.ti = {"eta_1": eta1, "director": ndir, "footprint": foot} return eta1, ndir, foot @@ -884,6 +1077,11 @@ def junction_cells(self, ring=1): raise RuntimeError( "junction cells are the SPLIT realisation's joints; the " "weak plane has no cut to fall off") + if self.dim == 3: + raise NotImplementedError( + "junction_cells is the 2-D ribbon rule (spines and cut " + "chains); in 3-D place the junction glue with " + "damage_yield — tubes about the junction segments.") from underworld3.utilities.place_surface import _cell_centroids_of band = self.info["band"] @@ -1086,18 +1284,36 @@ def slips(self, solver): def _slips_ti(self, solver): """The weak plane's slip: the tangential velocity jump across the - band, one half-width plus a cell either side of each spine.""" + band, one half-width plus a cell either side of each spine (2-D) + or patch (3-D, where the jump is projected onto the plane).""" import underworld3 as uw if self.info is None: raise RuntimeError("no band on this mesh: build(width=...)") out = {} for k, (name, P) in enumerate(self.prepared): - P = np.asarray(P, dtype=float)[:, :2] + P = np.asarray(P, dtype=float) + skirt = 0.5 * self.width + float(self.info["spacing"][k]) + if self.dim == 3: + # planar patch, one normal: probe the patch's own points + # (corners, edge midpoints, centroid) and keep the + # in-plane part of the jump — the plane form of the 2-D + # tangential projection + n_hat = _patch_normal(P) + S = np.vstack([P, 0.5 * (P + np.roll(P, -1, axis=0)), + P.mean(axis=0)]) + vp = np.asarray(uw.function.evaluate( + solver.u.sym, S + skirt * n_hat)).reshape(len(S), -1) + vm = np.asarray(uw.function.evaluate( + solver.u.sym, S - skirt * n_hat)).reshape(len(S), -1) + dv = (vp - vm)[:, :3] + dv -= np.outer(dv @ n_hat, n_hat) + out[name] = float(np.linalg.norm(dv, axis=1).max()) + continue + P = P[:, :2] t = np.gradient(P, axis=0) t /= np.linalg.norm(t, axis=1)[:, None] n = np.column_stack([-t[:, 1], t[:, 0]]) - skirt = 0.5 * self.width + float(self.info["spacing"][k]) vp = np.asarray(uw.function.evaluate( solver.u.sym, P + skirt * n)).reshape(len(P), -1)[:, :2] vm = np.asarray(uw.function.evaluate( diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index 6dcfb9498..c83b33a62 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -3847,7 +3847,8 @@ def vid(sheet, i, j): return pts, tets -def _occ_assembly_3d(patches, width, size, domain=None, assembly="fuse"): +def _occ_assembly_3d(patches, width, size, domain=None, assembly="fuse", + embed=None): """Thicken each planar patch by ±width/2, resolve overlaps, mesh. ``assembly`` is :func:`place_thin_volume`'s: ``"fuse"`` returns the union @@ -3860,10 +3861,23 @@ def _occ_assembly_3d(patches, width, size, domain=None, assembly="fuse"): reaching a boundary leaves its clipped face exactly on the boundary's own facets (snapped onto their planes after meshing, defensively). - Returns ``(points, tets, cad_volume)`` — the assembly mesh in its own - numbering, and the CAD volume of the (clipped) pieces, against which the - meshed volume is gated (planar-faced solids mesh to their exact volume). + ``embed`` is the 2-D network lesson one dimension up: a sequence of + planar polygons (typically the patches themselves, un-thickened) + FRAGMENTED INTO the fused solid before meshing, so each becomes a + conforming interior surface of the band — the mid-surface a split + walks, at any width. Requires ``domain=None`` for now (an embedded + surface with an outcropping band would need the same clip — refused). + + Returns ``(points, tets, cad_volume, embedded)`` — the assembly mesh + in its own numbering, the CAD volume of the (clipped) pieces against + which the meshed volume is gated (planar-faced solids mesh to their + exact volume), and per ``embed`` entry the ``(m, 3)`` triangles of + its embedded faces in assembly numbering (``None`` without ``embed``). """ + if embed is not None and domain is not None: + raise NotImplementedError( + "embedded mid-surfaces with a domain clip (outcropping bands) " + "are not built yet — the surfaces would need the same clip.") import gmsh gmsh.initialize() @@ -3899,6 +3913,23 @@ def _occ_assembly_3d(patches, width, size, domain=None, assembly="fuse"): solids = [t for _d, t in gmsh.model.getEntities(3)] tool, planes = _occ_domain_3d(occ, dom_verts, dom_tris) occ.intersect([(3, t) for t in solids], [(3, tool)]) + per_embed = None + if embed is not None: + occ.synchronize() + host = [t for _d, t in gmsh.model.getEntities(3)] + mids = [] + for poly in embed: + Q = np.asarray(poly, dtype=float) + mpts = [occ.addPoint(*q) for q in Q] + mlines = [occ.addLine(mpts[i], mpts[(i + 1) % len(Q)]) + for i in range(len(Q))] + mids.append(occ.addPlaneSurface([occ.addCurveLoop(mlines)])) + _frag, frag_map = occ.fragment([(3, t) for t in host], + [(2, t) for t in mids]) + # frag_map aligns with the input: host volumes first, then + # each mid-surface's descendants + per_embed = [[t for d, t in frag_map[len(host) + k] if d == 2] + for k in range(len(mids))] occ.synchronize() vols = gmsh.model.getEntities(3) @@ -3943,7 +3974,23 @@ def _occ_assembly_3d(patches, width, size, domain=None, assembly="fuse"): # band logic and the cap's node sharing need EXACT membership, # so snap defensively. xyz = _snap_to_boundary_3d(xyz, dom_verts, dom_tris, planes) - return xyz, tets, float(cad_volume) + embedded = None + if embed is not None: + embedded = [] + for k, faces in enumerate(per_embed): + tris = [] + for sf in faces: + et, _ei, en = gmsh.model.mesh.getElements(2, sf) + for ty, nodes in zip(et, en): + if ty == 2: + tris.append(np.array( + [renum[int(x)] for x in nodes], + dtype=np.int64).reshape(-1, 3)) + if not tris: + raise RuntimeError( + f"embedded surface {k} meshed to no faces") + embedded.append(np.vstack(tris)) + return xyz, tets, float(cad_volume), embedded finally: gmsh.finalize() @@ -6673,7 +6720,7 @@ def span_labels_at(m): def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, clearance=0.7, size=None, *, assembly="fuse", - mesher=None, verbose=False): + mesher=None, embed=None, verbose=False): """Embed a THIN VOLUME of the given width around each patch, junctions free. The finite-width fault representation: each planar patch is thickened by @@ -6755,6 +6802,15 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, layers and split to tets, no remesh (:func:`_ladder_assembly_3d`). Ladder bands must lie inside the domain. + embed : sequence of array_like, keyword-only, 3-D ``mesher="network"`` + Planar polygons FRAGMENTED INTO the fused band as conforming + interior surfaces — the mid-surfaces a split walks, at any width + (the 2-D network mesher's embedded spines, one dimension up). + Typically the un-expanded fault patches while ``patches`` carry + their margin-expanded bands. Interior assemblies only (no + outcrop clip yet). ``info["embedded_nodes"]`` returns one + ``(m, 3)`` node-coordinate array per entry — the point set + :func:`_label_mid_surface` labels a placed fault from. verbose : bool Report the counts. @@ -6764,7 +6820,8 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, A new mesh with the assembly's cells embedded verbatim. info : dict Global counts: ``n_zone_cells``, ``n_skin_faces``, ``n_placed`` - (vertices added), ``n_removed`` (vertices deleted), ``min_volume``. + (vertices added), ``n_removed`` (vertices deleted), ``min_volume``; + with ``embed``, also ``embedded_nodes``. Raises ------ @@ -6784,6 +6841,11 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, if mesher not in (None, "ladder", "network"): raise ValueError( f"mesher must be None, 'ladder' or 'network', not {mesher!r}") + if embed is not None and mesher != "network": + raise ValueError( + "embed= belongs to mesher='network' (mid-surfaces fragmented " + "into the fused band); any other mesher would silently " + "ignore it.") if dm.getDimension() == 2: return _place_thin_volume_2d(dm, patches, width, label, label_value, @@ -6824,21 +6886,35 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, payload = None if comm.rank == 0: try: + embedded_nodes = None if mesher == "ladder": # The extruded band: no CAD, no remesh — the sheet's own # triangulation offset to prisms. Interior bands only # (no domain clip); a protruding ladder fails the carve. asm_pts, asm_tets = _ladder_assembly_3d( patches[0][0], patches[0][1], width) + elif mesher == "network": + # fuse + embedded mid-surfaces (interior assemblies only): + # the fault surfaces become conforming faces of the band, + # so the split walks them at any width + asm_pts, asm_tets, _cad_vol, embedded = _occ_assembly_3d( + patches, width, size, domain=None, assembly=assembly, + embed=embed) + if embedded is not None: + # coordinates, not indices: the imprint collapse below + # may renumber, but interior points do not move + embedded_nodes = [ + np.asarray(asm_pts[np.unique(t)], dtype=float) + for t in embedded] else: # The meshed-vs-CAD volume gate runs inside the assembly # builder, before its boundary snap. - asm_pts, asm_tets, _cad_vol = _occ_assembly_3d( + asm_pts, asm_tets, _cad_vol, _no_embed = _occ_assembly_3d( patches, width, size, domain=(dom_verts, dom_tris), assembly=assembly) asm_pts, asm_tets = _collapse_boundary_imprints_3d( asm_pts, asm_tets, dom_verts, dom_tris, 0.1 * size) - payload = (asm_pts, asm_tets) + payload = (asm_pts, asm_tets, embedded_nodes) # Exception, not just RuntimeError/ValueError: a raw gmsh error # (e.g. a PLC intersection) is a plain Exception, and an # uncaught raise on the surgery rank is a HANG for its peers — @@ -6849,7 +6925,7 @@ def place_thin_volume(dm, patches, width, label=ZONE_LABEL, label_value=1, real = [f for f in failures if f] if real: raise RuntimeError(f"place_thin_volume assembly failed: {real[0]}") - asm_pts, asm_tets = comm.bcast(payload, root=0) + asm_pts, asm_tets, embedded_nodes = comm.bcast(payload, root=0) skin_xyz, skin_tris, skin_node_ids = _assembly_skin(asm_pts, asm_tets) # The outcrop band: the skin's trace on the domain boundary. @@ -7237,6 +7313,8 @@ def new_id(v): "n_placed": n_placed, "n_removed": n_removed, "n_trace_facets": int(len(band_idx)), "min_volume": float(min_vol[0])} + if embedded_nodes is not None: + info["embedded_nodes"] = embedded_nodes if verbose: uw.pprint(f"[place_thin_volume {label!r}] {info['n_zone_cells']} " f"zone cells, {info['n_skin_faces']} skin faces; placed " @@ -7337,15 +7415,22 @@ def _label_mid_surface(dm, spine_points, label, value): break for f in bad: del sel[f] - if not sel: + # Gather-first placement: the placed band lives on ONE rank, so the + # other ranks legitimately select no face. The refusal is judged on + # the GLOBAL count and raised collectively — a rank-local raise on + # an empty selection is a hang for the peers. + comm = dm.getComm().tompi4py() + n_global = int(comm.allreduce(len(sel))) + if n_global == 0: raise ValueError( f"no {label!r} faces survive on the mid-surface: the inset " f"leaves too small an interior (enlarge the sheet grid or " f"reduce inset_rings).") dm.createLabel(label) + lbl = dm.getLabel(label) for f in sel: - dm.getLabel(label).setValue(f, int(value)) - return len(sel) + lbl.setValue(f, int(value)) + return n_global def place_fault_ribbon(base_mesh, sheet, width, *, normals=None, diff --git a/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py b/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py new file mode 100644 index 000000000..06749e3ce --- /dev/null +++ b/tests/parallel/ptest_0863_fault_network_3d_width_parallel.py @@ -0,0 +1,69 @@ +"""The finite-width 3-D fault network in parallel: build the band, +split every patch, contact solve at np=2 — the serial answer, on a +distributed mesh. + +The band build gathers only for the rank-0 CAD assembly +(place_thin_volume) and the network split redistributes ONCE keyed on +the union of the network's facets (split_faults); everything else runs +distributed. Run with: + mpirun -np 2 python -m pytest tests/parallel/\ +ptest_0863_fault_network_3d_width_parallel.py --with-mpi +""" +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.parallel_safe, pytest.mark.level_2, + pytest.mark.tier_b] + +H = 0.08 +WIDTH = 0.04 +P_A = np.array([[0.30, 0.50, 0.30], [0.70, 0.50, 0.30], + [0.70, 0.50, 0.70], [0.30, 0.50, 0.70]]) +P_B = np.array([[0.30, 0.62, 0.32], [0.62, 0.30, 0.32], + [0.62, 0.30, 0.68], [0.30, 0.62, 0.68]]) +# the serial run of this exact case (test_0863's end-to-end fixture): +# peak tangential pair slip per prepared piece +SERIAL = {"Main": 0.17567, "Cross_1": 0.00335, "Cross_2": 0.00295} + + +def test_network_3d_width_split_solve_np2(): + fsA = uw.meshing.FaultSurface("Main", P_A) + fsA.triangulate() + fsB = uw.meshing.FaultSurface("Cross", P_B) + fsB.triangulate() + net = uw.meshing.FaultNetwork([fsA, fsB], + hierarchy=["Main", "Cross"]) + net.prepare(h=H, ligament=1.0, verbose=False) + net.build(width=WIDTH, realisation="split", h_far=0.24, + margin_rings=0.5) + mesh = net.mesh + + # the mesh is DISTRIBUTED (the far field balanced; only the CAD + # assembly is gathered), not serial-on-one-rank + comm = mesh.dm.comm.tompi4py() + local = int(mesh.dm.getHeightStratum(0)[1]) + assert comm.allreduce(local, op=min) > 0, "a rank holds no cells" + + x, y, z = mesh.X + v = uw.discretisation.MeshVariable("v3W", mesh, 3, degree=2) + p = uw.discretisation.MeshVariable("p3W", mesh, 1, degree=0, + continuous=False) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = [0.0, 0.0, 0.0] + for wall in ("Bottom", "Top", "Left", "Right", "Front", "Back"): + stokes.add_dirichlet_bc((y - 0.5, 0.0, 0.0), wall) + net.apply(stokes) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1e-5 + info = net.solve(stokes) + assert info.get("converged") + + slips = net.slips(stokes) # rank-local pairs + for name, expected in SERIAL.items(): + peak = comm.allreduce(float(slips.get(name, 0.0)), op=max) + assert peak == pytest.approx(expected, rel=2e-2), ( + f"{name}: parallel peak {peak:.4f} vs serial {expected}") diff --git a/tests/test_0863_fault_network_3d_width.py b/tests/test_0863_fault_network_3d_width.py new file mode 100644 index 000000000..049084456 --- /dev/null +++ b/tests/test_0863_fault_network_3d_width.py @@ -0,0 +1,249 @@ +"""The finite-width 3-D fault network: ONE band, both realisations. + +test_0858's property one dimension up: ``build(width=...)`` thickens +the margin-expanded patches into one fused band with each un-expanded +patch embedded as a conforming mid-surface, so the same mesh is cut +and split (``realisation="split"``) or left whole for the volumetric +weak plane (``"ti"``). Geometry oracles are analytic: the patches are +planar, so the honoured-footprint rule and the weak-plane director +(the patch normal) are exact. + +The parallel form is ``tests/parallel/ptest_0863``. +""" +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b, + pytest.mark.skipif(uw.mpi.size > 1, + reason="serial suite; the parallel " + "form is ptest_0863")] + +H = 0.08 +WIDTH = 0.04 + +P_A = np.array([[0.30, 0.50, 0.30], [0.70, 0.50, 0.30], + [0.70, 0.50, 0.70], [0.30, 0.50, 0.70]]) +P_B = np.array([[0.30, 0.62, 0.32], [0.62, 0.30, 0.32], + [0.62, 0.30, 0.68], [0.30, 0.62, 0.68]]) + + +def _network(realisation, width=WIDTH): + fsA = uw.meshing.FaultSurface("Main", P_A) + fsA.triangulate() + fsB = uw.meshing.FaultSurface("Cross", P_B) + fsB.triangulate() + net = uw.meshing.FaultNetwork([fsA, fsB], + hierarchy=["Main", "Cross"]) + net.prepare(h=H, ligament=1.0, verbose=False) + # margin_rings=0.5: the margin must stay half a cell clear of the + # ligament=1.0 junction gap (the build refuses a welding margin), + # and these small test patches in a unit box also need the + # expanded band to clear the walls for the carve + net.build(width=width, realisation=realisation, h_far=0.24, + margin_rings=0.5) + return net + + +def test_the_two_realisations_share_one_mesh(): + split = _network("split") + ti = _network("ti") + + assert split.realisation == "split" and ti.realisation == "ti" + assert [n for n, _p in split.prepared] == \ + [n for n, _p in ti.prepared] + assert split.info["n_cells"] == ti.info["n_cells"], ( + "the realisations no longer share a mesh; the comparison " + "between them is then confounded by the discretisation") + # the ONLY difference: the split duplicated the cut's nodes + assert (split.mesh.dm.getDepthStratum(0)[1] + > ti.mesh.dm.getDepthStratum(0)[1]) + + +def test_the_weak_plane_director_is_the_patch_normal(): + ti = _network("ti") + eta1, ndir, foot = ti.ti_fields(eta_1=0.01, eta_0=1.0) + + assert foot.sum() > 0, "no footprint cells: nothing would be weak" + assert foot.sum() <= ti.info["band"].sum() + vals = eta1.array[:, 0, 0] + assert np.allclose(vals[foot], 0.01) + assert np.allclose(vals[~foot], 1.0) + + # planar patches make the director exact: y-hat on Main's plane + # (y = 0.5), (1, 1, 0)/sqrt(2) on the Cross pieces' (x + y = 0.92) + d = np.asarray(ndir.array).reshape(-1, 3) + foots = ti.footprints + assert foots["Main"].any() + assert np.allclose(np.abs(d[foots["Main"]] @ [0.0, 1.0, 0.0]), 1.0) + nB = np.array([1.0, 1.0, 0.0]) / np.sqrt(2) + for name, m in foots.items(): + if name.startswith("Cross") and m.any(): + assert np.allclose(np.abs(d[m] @ nB), 1.0) + + +def test_the_band_is_damage_material_in_either_realisation(): + for realisation in ("split", "ti"): + net = _network(realisation) + assert net.band.sum() > 0 + assert set(net.footprints) == {n for n, _p in net.prepared} + assert net.footprints["Main"].sum() > 0 + assert net.footprints["Main"].sum() <= net.band.sum() + + tau = net.band_yield(tau_y=4.0) + painted = net._band_yield_var.array[:, 0, 0] + assert np.allclose(painted[net.band], 4.0) + assert np.allclose(painted[~net.band], 1.0e8) + assert tau.free_symbols # a usable expression + + # the junction glue in 3-D is damage_yield's tubes; the 2-D + # ribbon rule refuses rather than guessing + if realisation == "split": + with pytest.raises(NotImplementedError, match="damage_yield"): + net.junction_cells() + + +def test_the_build_refusals_are_loud(): + fsA = uw.meshing.FaultSurface("Main", P_A) + fsA.triangulate() + net = uw.meshing.FaultNetwork([fsA]) + net.prepare(h=H, verbose=False) + with pytest.raises(ValueError, match="width"): + net.build(realisation="ti") + with pytest.raises(ValueError, match="realisation"): + net.build(width=WIDTH, realisation="smeared") + # the band has ONE mesher; the no-band meshers are not it + with pytest.raises(ValueError, match="network"): + net.build(width=WIDTH, mesher="embed") + with pytest.raises(NotImplementedError, match="base box"): + net.build(width=WIDTH, base="a mesh") + + # a junction network refuses a margin that welds the junction gap: + # the gap is physics (the intact ligament), the margin convenience, + # and a welded gap degenerates the split (self-paired nodes) + fsB = uw.meshing.FaultSurface("Cross", P_B) + fsB.triangulate() + fsA2 = uw.meshing.FaultSurface("Main", P_A) + fsA2.triangulate() + netj = uw.meshing.FaultNetwork([fsA2, fsB], + hierarchy=["Main", "Cross"]) + netj.prepare(h=H, ligament=1.0, verbose=False) + with pytest.raises(ValueError, match="junction gap"): + netj.build(width=WIDTH, margin_rings=1) + + +def test_an_outcropping_band_is_refused(): + """A band that reaches the domain boundary cannot yet carry its + embedded mid-surface: refused loudly, before any meshing.""" + P_out = np.array([[0.30, 0.50, 0.30], [0.70, 0.50, 0.30], + [0.70, 0.50, 0.95], [0.30, 0.50, 0.95]]) + fs = uw.meshing.FaultSurface("Daylight", P_out) + fs.triangulate() + net = uw.meshing.FaultNetwork([fs]) + net.prepare(h=H, verbose=False) + with pytest.raises(NotImplementedError, match="boundary"): + net.build(width=WIDTH) + + +def test_embed_belongs_to_the_network_mesher(): + """embed= means mid-surfaces fragmented into the fused band; any + other mesher would silently ignore it, so it is refused.""" + from underworld3.utilities.place_surface import place_thin_volume + + box = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), + cellSize=0.5, qdegree=1) + with pytest.raises(ValueError, match="network"): + place_thin_volume(box.dm, [P_A], 0.05, embed=[P_A]) + + +def test_the_expansion_survives_collinear_rim_points(): + """A clipped rim may carry collinear consecutive points (the + convexity gate passes them); the in-plane offset must not divide + by zero there — on a straight rim the corner is the edge's own + offset — and the plane normal must not come from a collinear + leading triple (Newell's method).""" + from underworld3.meshing.fault_network import _expand_convex_polygon + + P = np.array([[0.30, 0.50, 0.30], [0.50, 0.50, 0.30], + [0.70, 0.50, 0.30], + [0.70, 0.50, 0.70], [0.30, 0.50, 0.70]]) + E = _expand_convex_polygon(P, 0.08) + assert np.isfinite(E).all() + # the straight rim's midpoint moves exactly one offset outward, + # staying in the patch plane + assert E[1] == pytest.approx([0.50, 0.50, 0.22]) + assert np.allclose(E[:, 1], 0.50) + + +def test_the_weak_plane_gauges_the_jump_across_its_layer(): + """The 3-D weak plane's slip is the in-plane velocity jump across + the layer. Read on a PRESCRIBED linear shear, where the jump is + known exactly: no solve, just the gauge's own arithmetic.""" + import types + + P = np.array([[0.30, 0.30, 0.50], [0.70, 0.30, 0.50], + [0.70, 0.70, 0.50], [0.30, 0.70, 0.50]]) + fs = uw.meshing.FaultSurface("Flat", P) + fs.triangulate() + net = uw.meshing.FaultNetwork([fs]) + net.prepare(h=H, verbose=False) + net.build(width=WIDTH, realisation="ti", h_far=0.24, + margin_rings=1) + + a = 3.0 # v = (a z, 0, 0): in-plane jump = 2 a skirt + v = uw.discretisation.MeshVariable("Ug", net.mesh, 3, degree=2) + v.array[:, 0, 0] = a * np.asarray(v.coords)[:, 2] + v.array[:, 0, 1] = 0.0 + v.array[:, 0, 2] = 0.0 + + skirt = 0.5 * WIDTH + float(net.info["spacing"][0]) + got = net.slips(types.SimpleNamespace(u=v)) + assert got["Flat"] == pytest.approx(2 * a * skirt, rel=1e-6) + + +def test_the_realisations_solve_on_the_shared_band(): + """Both realisations run on the one band under the same shear + drive: the split's contact releases the drive-aligned senior, and + the weak plane's TI solve converges on the same cells.""" + slips = {} + for realisation in ("split", "ti"): + net = _network(realisation) + mesh = net.mesh + x, y, z = mesh.X + v = uw.discretisation.MeshVariable(f"v3W_{realisation}", mesh, 3, + degree=2) + p = uw.discretisation.MeshVariable(f"p3W_{realisation}", mesh, 1, + degree=0, continuous=False) + stokes = uw.systems.Stokes(mesh, velocityField=v, + pressureField=p) + stokes.bodyforce = [0.0, 0.0, 0.0] + for wall in ("Bottom", "Top", "Left", "Right", "Front", "Back"): + stokes.add_dirichlet_bc((y - 0.5, 0.0, 0.0), wall) + if realisation == "split": + stokes.constitutive_model = \ + uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + net.apply(stokes) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1e-5 + info = net.solve(stokes) + assert info.get("converged") + else: + net.apply(stokes, eta_1=0.01) + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1e-5 + stokes.solve() + slips[realisation] = net.slips(stokes) + + for realisation, got in slips.items(): + assert all(np.isfinite(s) for s in got.values()), ( + realisation, got) + # the drive-aligned senior dominates the split + s = slips["split"] + assert s["Main"] > 0.05 + assert s["Main"] > 3 * max(v for n, v in s.items() if n != "Main") + # the two gauges are different quantities on different physics, but + # both are the layer's own throughput: same order of magnitude + assert 0.1 * s["Main"] < slips["ti"]["Main"] < 10 * s["Main"]