The general 3-D outcrop cap: a zone may meet any boundary

tests/parallel/ptest_0854_place_sheet_parallel.py

@@ -130,3 +130,56 @@ 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_an_outcropping_sheet_traces_partition_independently():
+    """The general-boundary outcrop in parallel: the trace chain's edge
+    count is a pure function of the sheet and the gathered boundary
+    complex, never of the partition. The clip and chain run identically
+    on every rank BEFORE the gather; the labelled count is read off the
+    result and must equal them at every rank count (the #589 discipline,
+    for the sheet's chain instead of the zone's band)."""
+    comm = uw.mpi.comm
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0),
+        cellSize=0.12, regular=False, qdegree=2)
+    strike = np.array([1.0, 0.0, 0.0])
+    dip = np.array([0.0, 0.15, -1.0])
+    dip /= np.linalg.norm(dip)
+    top = np.array([0.5, 0.5, 1.1])
+    s = np.linspace(-0.22, 0.22, 5)
+    d = np.linspace(0.0, 0.5, 5)
+    pts = np.array([top + a * strike + b * dip for b in d for a in s])
+    tris, n = [], 5
+    for i in range(n - 1):
+        for j in range(n - 1):
+            a, b = i * n + j, i * n + j + 1
+            c, e = (i + 1) * n + j, (i + 1) * n + j + 1
+            tris += [(a, b, e), (a, e, c)]
+
+    dm, info = place_sheet(mesh.dm, pts, np.array(tris, dtype=np.int64),
+                           label="FltA", label_value=41)
+    assert info["n_trace_edges"] > 0, "the sheet left no trace"
+    gathered = comm.allgather(info)
+    assert all(g == gathered[0] for g in gathered)
+
+    # Owned labelled trace EDGES across all ranks == the chain's count.
+    pStart, _pEnd = dm.getChart()
+    eS, eE = dm.getDepthStratum(1)
+    leaves = np.zeros(dm.getChart()[1] - pStart, dtype=bool)
+    try:
+        _n, ilocal, _ir = dm.getPointSF().getGraph()
+        if ilocal is not None and len(ilocal):
+            leaves[np.asarray(ilocal, dtype=np.int64) - pStart] = True
+    except (ValueError, TypeError):
+        pass
+    n_local = 0
+    trace = dm.getLabel("FltA_trace")
+    if dm.hasLabel("FltA_trace") and trace.getStratumSize(41) > 0:
+        for p in trace.getStratumIS(41).getIndices():
+            if eS <= int(p) < eE and not leaves[int(p) - pStart]:
+                n_local += 1
+    n_edges = int(comm.allreduce(n_local, op=MPI.SUM))
+    assert n_edges == info["n_trace_edges"], (
+        f"{n_edges} owned trace edges for {info['n_trace_edges']} in the "
+        "chain; the trace depends on the partition")

tests/parallel/ptest_0855_place_thin_volume_parallel.py

@@ -187,6 +187,59 @@ def test_a_2d_outcrop_on_the_annulus_embeds_in_parallel():
     assert n_bare == 0, "the relabel left boundary edges without a label"
 
 
+def test_a_3d_outcrop_on_a_rotated_box_embeds_in_parallel():
+    """The general 3-D cap at np>=2: no wall axis-aligned.
+
+    The same mesh, rotation and patch as the serial
+    test_0860_outcrop_general_boundary_3d, so the trace count is a
+    partition-independence check against the serial value by suite. The
+    boundary complex is gathered from UNSHARED support-1 faces and sorted
+    canonically, so the frame masks, the collar and the info dict are
+    functions of the mesh, not of the partition.
+    """
+    comm = uw.mpi.comm
+    a, b = np.deg2rad(20.0), np.deg2rad(15.0)
+    Rz = np.array([[np.cos(a), -np.sin(a), 0.0],
+                   [np.sin(a), np.cos(a), 0.0], [0.0, 0.0, 1.0]])
+    Rx = np.array([[1.0, 0.0, 0.0], [0.0, np.cos(b), -np.sin(b)],
+                   [0.0, np.sin(b), np.cos(b)]])
+    R = Rz @ Rx
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0),
+        cellSize=0.12, regular=False, qdegree=2)
+    vec = mesh.dm.getCoordinatesLocal()
+    vec.array[:] = (vec.array.reshape(-1, 3) @ R.T).reshape(-1)
+    mesh.dm.setCoordinatesLocal(vec)
+    patch = np.array([[0.3, 0.45, 1.1], [0.7, 0.45, 1.1],
+                      [0.72, 0.55, 0.55], [0.32, 0.55, 0.55]]) @ R.T
+    new, info = place_thin_volume(mesh.dm, [patch], width=0.05,
+                                  label="Zone", label_value=5)
+    assert info["n_zone_cells"] > 0
+    assert info["n_trace_facets"] > 0, "the zone left no trace"
+    assert _owned_label_count(new, "Zone", 5) == info["n_zone_cells"]
+    assert _owned_label_count(new, "Zone_skin", 5) == info["n_skin_faces"]
+    assert _owned_label_count(new, "Zone_trace", 5) > 0
+    gathered = comm.allgather(info)
+    assert all(g == gathered[0] for g in gathered)
+
+    # Every true boundary face still carries a wall label — the relabel
+    # restored the rotated wall's labels per removed face.
+    from underworld3.utilities.place_surface import _shared_point_flags
+    pStart, _pEnd = new.getChart()
+    shared = _shared_point_flags(new).astype(bool)
+    fS, fE = new.getHeightStratum(1)
+    walls = [new.getLabel(n)
+             for n in ("Top", "Bottom", "Left", "Right", "Front", "Back")]
+    n_bare = 0
+    for f in range(fS, fE):
+        if shared[f - pStart] or len(new.getSupport(f)) != 1:
+            continue
+        if all(w.getValue(f) < 0 for w in walls):
+            n_bare += 1
+    n_bare = int(comm.allreduce(n_bare, op=MPI.SUM))
+    assert n_bare == 0, "the relabel left boundary faces without a label"
+
+
 def test_2d_refusals_are_collective():
     """A ribbon stopping short of the wall refuses IDENTICALLY everywhere."""
     comm = uw.mpi.comm

tests/test_0854_place_sheet.py

@@ -267,3 +267,39 @@ def test_an_outcropping_sheet_meets_the_top_surface():
     err = np.abs(np.asarray(t.data[:, 0])
                  - (X[:, 0]**2 + X[:, 1]**2 + X[:, 2]**2))
     assert float(err.max()) < 1e-8
+
+def test_a_sheet_embeds_in_a_spherical_shell_and_removes_again():
+    """The domain's topology is conserved, not assumed.
+
+    A spherical shell is S^2 x I, global Euler number 2, and the old gate
+    demanded 1 — refusing the domain for its topology rather than for any
+    defect of the surgery. Placement and removal both run their own
+    volume, conformity and Euler-conservation gates; asserting the info
+    here proves they ran and agreed. Serial guard matches the suite.
+    """
+    if uw.mpi.size > 1:
+        pytest.skip("serial suite; the parallel form is ptest_0854")
+    from underworld3.utilities.place_surface import remove_embedded
+
+    shell = uw.meshing.SphericalShell(radiusInner=0.25, radiusOuter=1.0,
+                                      cellSize=0.13, qdegree=2)
+    d1 = np.array([0.0, 1.0, 0.0])
+    d2 = np.array([-1.0, 0.0, 1.0]) / np.sqrt(2.0)
+    C = np.array([0.45, 0.0, 0.45])
+    n, s = 5, np.linspace(-0.1, 0.1, 5)
+    pts = np.array([C + a * d1 + b * d2 for a in s for b in s])
+    tris = []
+    for i in range(n - 1):
+        for j in range(n - 1):
+            a, b = i * n + j, i * n + j + 1
+            c, d = (i + 1) * n + j, (i + 1) * n + j + 1
+            tris += [(a, b, d), (a, d, c)]
+
+    new, info = place_sheet(shell.dm, pts, np.array(tris, dtype=np.int64),
+                            label="Sheet", label_value=3)
+    assert info["n_surface_facets"] == len(tris)
+    assert new.getLabel("Sheet").getStratumSize(3) >= len(tris)
+
+    back, rinfo = remove_embedded(new, "Sheet", label_value=3)
+    assert rinfo["n_removed_cells"] > 0
+    assert back.getLabel("Sheet").getStratumSize(3) == 0

tests/test_0859_domain_boundary_tool.py

@@ -140,12 +140,15 @@ def test_the_3d_snap_restores_wall_values_exactly():
                       [0.0, 1.0, 0.0]])
     tris = np.array([[0, 1, 2], [0, 2, 3]])
     planes = [(verts[0], np.array([0.0, 0.0, 1.0]))]
+    # The control sits ABOVE the snap tolerance (1e-6 since the snap was
+    # widened to cover OCC's ~4e-7 boolean placement noise) and below any
+    # real feature scale — it must come through untouched.
     xyz = np.array([[0.3, 0.4, 1e-10],
                     [0.6, 0.2, -1e-10],
-                    [0.5, 0.5, 2e-9],
+                    [0.5, 0.5, 2e-5],
                     [0.5, 0.5, 0.5]])
     out = _snap_to_boundary_3d(xyz.copy(), verts, tris, planes)
     assert out[0, 2] == 0.0 and out[0, 0] == 0.3 and out[0, 1] == 0.4
     assert out[1, 2] == 0.0
-    assert out[2, 2] == 2e-9              # untouched: the control
+    assert out[2, 2] == 2e-5              # untouched: the control
     assert np.array_equal(out[3], xyz[3])