diff --git a/src/underworld3/meshing/smoothing/graph.py b/src/underworld3/meshing/smoothing/graph.py index 27d405ef..456f2fe1 100644 --- a/src/underworld3/meshing/smoothing/graph.py +++ b/src/underworld3/meshing/smoothing/graph.py @@ -527,33 +527,76 @@ def _min_incident_edge_nd(cells, coords): # ``boundary_slip`` orchestration, and ``BoundingSurface`` facet restore. # --------------------------------------------------------------------------- -def _slip_normals(mesh, boundary_coords): - """Unit outward normals at ``boundary_coords`` from the projected - boundary-normal field. - - Re-projects ``mesh._projected_normals`` (``mesh.Gamma_P1``) first so the - normals reflect the mesh's *current* coordinates — the projected field is - stale after any deform. Returns ``(normals, valid)`` where ``normals`` is - ``(k, cdim)`` and ``valid`` is a boolean mask; ``valid`` is ``False`` for - nodes with a degenerate (zero / non-finite) normal (e.g. box corners - where opposing face normals cancel, or an occasional unlocatable vertex). - Such nodes should be pinned, not slipped. +#: Composite labels that cover every boundary at once. Summing a per-boundary +#: normal over one of these would give a contribution at every boundary node and +#: make each of them look like a corner, so they are never iterated. +_COMPOSITE_BOUNDARY_LABELS = frozenset({"All_Boundaries"}) + + +def _slip_normals(mesh, boundary_coords, boundaries=None): + """Unit outward normals at ``boundary_coords``, taken PER BOUNDARY. + + Returns ``(normals, valid)``. ``valid`` is ``False`` where no boundary claims + the node, where the normal is degenerate, and -- deliberately -- where MORE + THAN ONE boundary claims it. Those are corners and 3-D edges: they have no + single normal, so they are pinned rather than slipped along a fabricated one. + + This used to evaluate ``mesh.Gamma_P1`` once. That field is deprecated and + documented as falling back to a coordinate direction off-kernel; on a unit + box it returns the edge TANGENT on Left and Bottom -- ``(0, 1)`` where the + outward normal is ``(-1, 0)``. A node whose "normal" lies along its own edge + has its tangential motion projected out, which is exactly why + ``node_redistribution(..., slip_surfaces=True)`` left those two edges frozen + at *exactly* zero displacement while Right and Top moved (#538). Per + boundary, ``_assemble_boundary_normal`` is exact on all four edges to 1.6e-11. + + NO NEW MeshVariable IS CREATED. The per-boundary normal is assembled into + ``_n_proj`` -- the field ``Gamma_P1`` already owns -- one boundary at a time. + That is not tidiness: creating even ONE extra variable on this path + restructures the DM, and the parent's static hierarchy stops matching what it + recorded (``test_0764::test_redistribute_then_adapt_composition``). Verified + by keeping the creation and reverting to the old normals, which still failed; + the creation is the trigger, not the normals. ``_n_proj`` is rebuilt by + ``_update_projected_normals`` on every call, so overwriting it here is safe. + + ``boundaries`` names the labels to consider; ``None`` uses every + non-composite boundary. Prefer passing the slip set: each label costs one + evaluation. """ cdim = mesh.cdim - n = np.zeros((boundary_coords.shape[0], cdim)) + count = boundary_coords.shape[0] + total = np.zeros((count, cdim)) + claims = np.zeros(count, dtype=int) + + if boundaries is None: + boundaries = [b.name for b in mesh.boundaries + if b.name not in _COMPOSITE_BOUNDARY_LABELS] + try: - mesh._update_projected_normals() - n = np.asarray( - uw.function.evaluate(mesh.Gamma_P1, boundary_coords) - ).reshape(-1, cdim) + mesh._update_projected_normals() # creates/refreshes _n_proj + scratch = mesh._projected_normals except Exception: - # Projection unavailable / degenerate on this mesh — fall back to - # all-pinned boundaries (valid stays all-False below). - n = np.zeros((boundary_coords.shape[0], cdim)) - mag = np.linalg.norm(n, axis=1) - valid = np.isfinite(mag) & (mag > 0.5) - out = np.zeros_like(n) - out[valid] = n[valid] / mag[valid, None] + return np.zeros((count, cdim)), np.zeros(count, dtype=bool) + + for name in boundaries: + if name in _COMPOSITE_BOUNDARY_LABELS: + continue + try: + mesh._assemble_boundary_normal(scratch, name) + here = np.asarray( + uw.function.evaluate(scratch.sym, boundary_coords) + ).reshape(-1, cdim) + except Exception: + continue # label absent on this rank / this mesh: claims nothing + magnitude = np.linalg.norm(here, axis=1) + on_it = np.isfinite(magnitude) & (magnitude > 0.5) + total[on_it] += here[on_it] + claims[on_it] += 1 + + magnitude = np.linalg.norm(total, axis=1) + valid = (claims == 1) & np.isfinite(magnitude) & (magnitude > 0.5) + out = np.zeros_like(total) + out[valid] = total[valid] / magnitude[valid, None] return out, valid diff --git a/tests/test_0764_node_redistribution.py b/tests/test_0764_node_redistribution.py index f36c5934..f6853a1e 100644 --- a/tests/test_0764_node_redistribution.py +++ b/tests/test_0764_node_redistribution.py @@ -222,3 +222,42 @@ def n_cells(m): assert n_cells(c_default) == n_cells(c_nvb) assert np.array_equal(np.asarray(c_default.X.coords), np.asarray(c_nvb.X.coords)) + + +def test_slip_normals_are_outward_on_every_wall(): + """Every wall gets its OWN outward normal, and corners are pinned (#538). + + `_slip_normals` read the deprecated `mesh.Gamma_P1`, which off-kernel falls + back to a coordinate direction rather than a normal: on a unit box it + returned the edge TANGENT on Left and Bottom -- (0, 1) where the outward + normal is (-1, 0). A node whose "normal" lies along its own edge has its + tangential motion projected out, so `node_redistribution(slip_surfaces=True)` + froze those two walls at EXACTLY zero displacement while Right and Top moved + (left 0.00000 / right 0.18241, bottom 0.00000 / top 0.18427 on a metric band + placed symmetrically about the domain centre). + + All four walls are asserted, not one: the defect was confined to two of them, + and a test on Right or Top alone passes throughout. + """ + from underworld3.meshing.smoothing.graph import _slip_normals + + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8, qdegree=3) + along = np.linspace(0.2, 0.8, 4) + zeros, ones = np.zeros_like(along), np.ones_like(along) + + for wall, points, expected in ( + ("Left", np.column_stack([zeros, along]), (-1.0, 0.0)), + ("Bottom", np.column_stack([along, zeros]), (0.0, -1.0)), + ("Right", np.column_stack([ones, along]), (1.0, 0.0)), + ("Top", np.column_stack([along, ones]), (0.0, 1.0))): + normals, valid = _slip_normals(mesh, points) + assert valid.all(), f"{wall}: {(~valid).sum()} of {len(points)} nodes unusable" + assert np.allclose(normals, np.array(expected), atol=1.0e-9), ( + f"{wall}: got {normals[0]}, expected {expected} — a tangent here " + "means the deprecated Gamma_P1 path is back") + + # A corner belongs to two walls and has no single normal: pin, do not slip. + corners = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) + _normals, valid = _slip_normals(mesh, corners) + assert not valid.any(), "corners must be pinned, not slipped along one wall"