From 8a9d2ff2c4c84a1a94628f6cb64ca8582397d98f Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 21 Jun 2026 07:55:35 +1000 Subject: [PATCH 1/2] fix(mesh): boundary normals and domain membership must track a deformed mesh Two related bugs where analytic/cached UNDEFORMED geometry survived a mesh.deform(), corrupting Nitsche/penalty free-slip BCs and semi-Lagrangian advection on free surfaces. (1) Projected boundary normals stayed radial on a deformed surface. mesh.Gamma's base scalars carry the C-code petsc_n[] (PETSc facet normal), defined only inside surface-integral kernels; _update_projected_normals point-evaluated it, so it fell back to the coordinate -> (x,y)/r = radial. Gamma_P1 therefore never tracked deformation and add_nitsche_bc constrained v.r_hat instead of v.n_hat_true, leaking throughflow proportional to surface tilt. Fix: new Mesh.boundary_normal(label) assembles the exact PETSc facet normals (computeCellGeometryFVM, outward, area-weighted) from ONLY that boundary's faces -> per-boundary (a corner shared with another boundary keeps its own one-sided normal, never an average across the discontinuity), deformation-tracking, smooth on a smooth boundary. add_nitsche_bc (Stokes and Scalar) now uses it; global Gamma_P1 is unchanged for back-compat. The field is refreshed on every deform so BCs that captured its .sym at setup read the new geometry. Parallel-safe (guards getStratumIS against dead label values on ranks holding no faces of a boundary). (2) points_in_domain()/function.evaluate() used the original boundary. The boundary-skeleton kd-tree was built from mesh._nav_coords, captured as a reference to the ORIGINAL coords in __init__ and never refreshed in nuke_coords_and_rebuild (only the per-cell control-point arrays were, on adapt). On a volume mesh _nav_coords stayed at the undeformed boundary, so a region that bulged past the original surface read as EXTERIOR -- stranding SL trace-back feet there and mis-locating evaluations, injecting the cold boundary value at topographic highs. Fix: invalidate the boundary kd-tree and refresh _nav_coords from the current DM coordinates on every deform. Tests: test_0056 (boundary_normal tracks a deformed annulus + a Cartesian free-surface corner) and test_0057 (deformed-domain membership + evaluation), both tier_a. test_1060 Nitsche free-slip passes serial and np2. Follow-up: ADD-reduce the unnormalised facet contributions at partition-seam boundary vertices before normalising (those few vertices currently use a this-rank-only stencil; parallel-safe and correct elsewhere). Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 22 ++- .../discretisation/discretisation_mesh.py | 175 ++++++++++++++++-- tests/test_0056_projected_normals_deform.py | 161 ++++++++++++++++ tests/test_0057_deformed_domain_membership.py | 78 ++++++++ 4 files changed, 416 insertions(+), 20 deletions(-) create mode 100644 tests/test_0056_projected_normals_deform.py create mode 100644 tests/test_0057_deformed_domain_membership.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 4848c719..29c5849e 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -2928,10 +2928,11 @@ class SNES_Vector(SolverBaseClass): # in the embedding space with an implicit tangency constraint. dim = mesh.cdim - # Surface normal components — use projected P1 normals by default. - # These are smooth, consistently oriented, and converge in 3D. - Gamma_P1 = mesh.Gamma_P1 - n = [Gamma_P1[i] for i in range(dim)] + # Surface normal components — use this boundary's own deformation- + # tracking facet normal (see Mesh.boundary_normal); the legacy global + # mesh.Gamma_P1 stays radial on a deformed surface. + bnorm = mesh.boundary_normal(boundary) + n = [bnorm[i] for i in range(dim)] # Constraint direction: defaults to surface normal if direction is not None: @@ -4771,16 +4772,21 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): mesh = self.mesh dim = mesh.dim - # Surface normal components. By default use projected P1 normals - # (smooth, consistently oriented, converges in 3D). + # Surface normal components. By default use this boundary's OWN + # deformation-tracking facet normal (assembled from only this + # boundary's faces, so a corner shared with another boundary is not + # averaged across the discontinuity, and a deformed surface is + # followed correctly). The legacy global mesh.Gamma_P1 point-evaluates + # petsc_n and stays radial on a deformed surface — see + # Mesh.boundary_normal. if normal is not None: if isinstance(normal, sympy.MatrixBase): n = [normal[i] for i in range(dim)] else: n = list(normal) else: - Gamma_P1 = mesh.Gamma_P1 - n = [Gamma_P1[i] for i in range(dim)] + bnorm = mesh.boundary_normal(boundary) + n = [bnorm[i] for i in range(dim)] # Constraint direction: defaults to surface normal if direction is not None: diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 4c003b07..585e1798 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2364,6 +2364,10 @@ def nuke_coords_and_rebuild( # Invalidate projected boundary normals (rebuilt lazily on access) self._projected_normals = None + # Per-boundary deformation-tracking normals are stale now too. The + # variables persist (we keep the name->var map); their DATA is refilled + # eagerly by Mesh.deform() after the remesh completes, so BCs that + # captured boundary_normal(...).sym at setup read the new geometry. # BUGFIX(#135): invalidate the per-cell face control-point arrays. # These are populated lazily by _get_mesh_face_control_points, sized @@ -2375,6 +2379,39 @@ def nuke_coords_and_rebuild( self.faces_outer_control_points = None self.faces_inner_control_points = None + # BUGFIX (deformed-domain membership): also invalidate the boundary- + # skeleton kd-tree used by points_in_domain() (and the SL out-of-domain + # restore). It is cached from the boundary geometry and was only rebuilt + # on adapt. After a DEFORM the surface has moved, so the stale control + # points (at the OLD boundary) wrongly flag a bulged-out region (r>r_o + # on a free surface) as EXTERIOR — stranding semi-Lagrangian trace-back + # feet there and mis-locating evaluations, which injects the cold + # boundary value at the topographic highs (upwellings). Rebuilt lazily + # from the deformed boundary on next points_in_domain() access. + self.boundary_face_control_points_kdtree = None + self.boundary_face_control_points_sign = None + self._domain_radius_squared = float("inf") + + # The navigation coords (used to build those control points and for + # point location) were captured as a reference to the ORIGINAL coords + # in __init__ and never refreshed here, so on a volume mesh they stayed + # at the undeformed geometry — the real reason a bulged-out region read + # as exterior. Re-point them at the current (deformed) DM coordinates. + if getattr(self, "_nav_dm", None) is None: + self._nav_coords = numpy.asarray( + self.dm.getCoordinatesLocal().array + ).reshape(-1, self.cdim) + else: + # manifold mesh: the nav clone carries its own (ghosted) coords; + # refresh them from the rebuilt main DM where possible. + try: + self._nav_dm.setCoordinatesLocal(self.dm.getCoordinatesLocal()) + self._nav_coords = numpy.asarray( + self._nav_dm.getCoordinatesLocal().array + ).reshape(-1, self.cdim) + except Exception: + pass + # BUGFIX(#130): recompute the DOF coordinate cache for every # already-registered variable. Variables created before this # rebuild would otherwise have their cache entry (from __init__) @@ -2420,27 +2457,23 @@ def _update_projected_normals(self): thereafter. The result is a smooth, consistently-oriented unit normal field that works well for penalty and Nitsche BCs on curved boundaries. + + NOTE: this GLOBAL field point-evaluates ``mesh.Gamma`` whose petsc_n + only exists in surface-integral kernels, so it falls back to the + coordinate (radial for a circle) and does NOT track a deformed + surface. For deformation-aware, corner-correct normals use the + per-boundary :meth:`boundary_normal` (which ``add_nitsche_bc`` now + uses). This global field is retained unchanged for back-compat. """ import underworld3 as uw Gamma = self.Gamma if not hasattr(self, '_projected_normals') or self._projected_normals is None: - # nuke_coords_and_rebuild() clears this attribute on every deform, - # but the underlying MeshVariable persists on the mesh. Recover it - # quietly rather than re-creating (which logs a name collision). existing = self.vars.get("_n_proj") if existing is not None: self._projected_normals = existing else: - # REINIT policy: _n_proj is a recomputed-on-access work - # variable. The values depend on the *current* mesh - # geometry (Gamma normals projected onto P1 nodes), so - # the Phase-1 remesh helper should NOT REMAP — the value - # at old node positions is wrong on the new mesh. - # Recomputed below from Gamma; reset on every deform - # via the _projected_normals=None clearing in - # nuke_coords_and_rebuild. self._projected_normals = uw.discretisation.MeshVariable( "_n_proj", self, self.cdim, degree=1, remesh_policy="reinit", @@ -2454,6 +2487,112 @@ def _update_projected_normals(self): nonzero = mag > 1.0e-30 n.data[nonzero] /= mag[nonzero, numpy.newaxis] + def boundary_normal(self, boundary): + """Outward unit normal of a single boundary, tracking deformation. + + Assembles the EXACT, outward, area-weighted PETSc facet normals + (``dm.computeCellGeometryFVM``) from ONLY this boundary's facets onto + its P1 vertices. Because each boundary is assembled independently, + a vertex shared by two boundaries (a sharp corner) is NOT averaged + across the discontinuity — each boundary keeps its own normal. On a + smooth boundary (e.g. a free surface) the result is the smooth + deformed normal. Cached per boundary; rebuilt lazily after a deform. + + Returns a sympy Matrix (row) of the P1 normal-field components, for + use as the constraint direction in Nitsche/penalty BCs. + + Parameters + ---------- + boundary : str or enum + Boundary label name (or a ``mesh.boundaries`` enum member). + """ + import underworld3 as uw + from scipy.spatial import cKDTree + + name = getattr(boundary, "name", str(boundary)) + if not hasattr(self, "_boundary_normal_vars") or self._boundary_normal_vars is None: + self._boundary_normal_vars = {} + var = self._boundary_normal_vars.get(name) + if var is None: + existing = self.vars.get(f"_n_bd_{name}") + var = existing if existing is not None else uw.discretisation.MeshVariable( + f"_n_bd_{name}", self, self.cdim, degree=1, remesh_policy="reinit") + self._boundary_normal_vars[name] = var + self._assemble_boundary_normal(var, name) + return var.sym + + def _assemble_boundary_normal(self, var, name): + """Fill ``var`` with the area-weighted outward facet normal assembled + from the faces of boundary ``name`` only (see :meth:`boundary_normal`).""" + import underworld3 as uw + from scipy.spatial import cKDTree + cdim = self.cdim + dm = self.dm + coords = numpy.ascontiguousarray(var.coords) + accum = numpy.zeros((coords.shape[0], cdim)) + + # faces carrying this boundary label: DM label named after the boundary, + # stratum keyed by the boundary's value (same access the BC code uses). + bvalue = None + for b in (self.boundaries or []): + if b.name == name: + bvalue = b.value + break + face_pts = [] + label = dm.getLabel(name) if dm.hasLabel(name) else None + if label is not None and bvalue is not None: + # NB: getStratumIS(value) for a value NOT in this rank's live value + # set can hard-abort the interpreter (e.g. a rank holding no faces + # of this boundary in parallel). Only query a live value. + try: + vis = label.getValueIS() + live = set(int(x) for x in vis.getIndices()) if vis.getSize() else set() + except Exception: + live = set() + if int(bvalue) in live: + pis = label.getStratumIS(bvalue) + if pis is not None and pis.getSize(): + fS, fE = dm.getHeightStratum(1) + for p in pis.getIndices(): + if fS <= int(p) < fE: + face_pts.append(int(p)) + + tree = cKDTree(coords) + nverts = cdim # P1 vertices per facet (2D edge=2, 3D tri=3) + for f in face_pts: + if dm.getSupportSize(f) != 1: + continue + area, cent, nrm = dm.computeCellGeometryFVM(f) + nrm = numpy.asarray(nrm)[:cdim] + cell = dm.getSupport(f)[0] + _, ccent, _ = dm.computeCellGeometryFVM(cell) + if numpy.dot(nrm, numpy.asarray(cent)[:cdim] + - numpy.asarray(ccent)[:cdim]) < 0: + nrm = -nrm + # Accumulate to the facet's P1 DOFs (its vertices) — found as the + # nearest DOFs to the facet centroid. This avoids indexing the local + # coordinate array by (vertex_point - vStart), which is only valid + # in serial (the parallel coordinate layout differs → out-of-range). + # Normalisation at the end makes the per-DOF weight (full vs share) + # irrelevant to the resulting direction. + _, idxs = tree.query(numpy.asarray(cent)[:cdim], k=nverts) + for idx in numpy.atleast_1d(idxs): + accum[idx] += area * nrm + + # TODO(parallel): a boundary vertex on a partition interface should + # ADD-reduce the UNnormalised facet contributions from both ranks before + # normalising (DMLocalToGlobal ADD_VALUES on the variable's section), + # so its normal is the full-stencil average rather than this rank's + # partial stencil. This is parallel-SAFE as-is (rank-interior boundary + # vertices are exact; only the handful of partition-seam surface + # vertices get a slightly-rotated unit normal). A first ADD-reduce + # attempt SEGV'd on the lazily-built work variable's global vec; deferred + # to a focused follow-up with the right vec/section plumbing. + mag = numpy.sqrt(numpy.sum(accum ** 2, axis=1)) + nonzero = mag > 1.0e-30 + accum[nonzero] /= mag[nonzero, numpy.newaxis] + var.data[...] = accum + @property def Gamma_P1(self): """Projected P1 boundary normals as a sympy Matrix. @@ -2897,9 +3036,21 @@ def deform(self, new_coords, *, dt=None, zero=None, verbose=False): def _do_move(): self._deform_mesh(_nc) - return remesh_with_field_transfer( + result = remesh_with_field_transfer( self, _do_move, dt=dt, extra_zero=zero, verbose=verbose) + # Refresh deformation-tracking per-boundary normals so Nitsche/penalty + # BCs that captured ``boundary_normal(...).sym`` at setup read the new + # geometry (the JIT reads the variable's .data, which would otherwise + # hold the setup-time normal). Re-assemble each cached boundary normal. + if getattr(self, "_boundary_normal_vars", None): + for _nm, _var in list(self._boundary_normal_vars.items()): + try: + self._assemble_boundary_normal(_var, _nm) + except Exception: + pass + return result + def _deform_mesh(self, new_coords: numpy.ndarray, verbose=False, active_vars=None): """ diff --git a/tests/test_0056_projected_normals_deform.py b/tests/test_0056_projected_normals_deform.py new file mode 100644 index 00000000..53f51222 --- /dev/null +++ b/tests/test_0056_projected_normals_deform.py @@ -0,0 +1,161 @@ +"""Per-boundary normals (mesh.boundary_normal) must track a DEFORMED surface +and must NOT be contaminated by a neighbouring boundary at a corner. + +Regression for the deformed-normal bug: the global mesh.Gamma_P1 point-evaluates +the mesh.Gamma base scalars, whose C-code maps to the PETSc facet normal +petsc_n[] — defined ONLY inside surface-integral kernels. A general point +evaluation fell back to the coordinate value, so the "normal" came out +≈ (x,y)/r = RADIAL (or the box coordinate) no matter how the mesh deformed. +Every Nitsche/penalty free-slip BC then constrained v·(stale normal) on a +deformed surface, leaking throughflow ∝ surface tilt. + +The fix is mesh.boundary_normal(label): assemble the EXACT PETSc facet normals +(dm.computeCellGeometryFVM) area-weighted onto the vertices of ONLY that +boundary's faces. Smooth boundary → smooth deformed normal; a corner shared +with another boundary keeps each boundary's own (one-sided) normal instead of +averaging across the discontinuity. + +Covers (a) a curved free surface (annulus Upper), and (b) a CARTESIAN free +surface (box Top deformed, vertical side walls) — the corner case. + +Run: pixi run python -m pytest tests/test_0056_projected_normals_deform.py -v +""" + +import pytest +import numpy as np +import underworld3 as uw + +pytestmark = [ + pytest.mark.level_1, + pytest.mark.tier_a, + # Geometric assembly is verified against a rank-local reference; run serial. + pytest.mark.skipif(uw.mpi.size > 1, reason="serial geometric-assembly check"), +] + + +def _facet_vertex_normals(mesh, label_name, label_value): + """Ground-truth: area-weighted outward vertex normals from the EXACT PETSc + facet normals of one boundary's faces, computed independently of the code + under test.""" + dm = mesh.dm + cdim = mesh.cdim + coords = np.asarray(mesh.X.coords) + accum = np.zeros_like(coords) + from scipy.spatial import cKDTree + tree = cKDTree(coords) + vS, vE = dm.getDepthStratum(0) + xc = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, cdim) + fS, fE = dm.getHeightStratum(1) + label = dm.getLabel(label_name) + pis = label.getStratumIS(label_value) + for f in pis.getIndices(): + if not (fS <= int(f) < fE) or dm.getSupportSize(int(f)) != 1: + continue + area, cent, nrm = dm.computeCellGeometryFVM(int(f)) + nrm = np.asarray(nrm)[:cdim] + cell = dm.getSupport(int(f))[0] + _, ccent, _ = dm.computeCellGeometryFVM(cell) + if np.dot(nrm, np.asarray(cent)[:cdim] - np.asarray(ccent)[:cdim]) < 0: + nrm = -nrm + for v in dm.getCone(int(f)): + if vS <= v < vE: + _, idx = tree.query(xc[v - vS], k=1) + accum[idx] += area * nrm + mag = np.sqrt((accum ** 2).sum(1)) + good = mag > 1e-30 + accum[good] /= mag[good, None] + return accum, good + + +def _angle_deg(a, b): + return np.degrees(np.arccos(np.clip(np.abs((a * b).sum(1)), -1.0, 1.0))) + + +def _bvalue(mesh, name): + for b in mesh.boundaries: + if b.name == name: + return b.value + raise KeyError(name) + + +def _eval_boundary_normal(mesh, name, pts): + bn = mesh.boundary_normal(name) + gx = np.asarray(uw.function.evaluate(bn[0], pts)).flatten() + gy = np.asarray(uw.function.evaluate(bn[1], pts)).flatten() + return np.column_stack([gx, gy]) + + +def test_boundary_normal_tracks_deformed_annulus(): + r_i, r_o, cs = 0.5, 1.0, 0.2 + mesh = uw.meshing.Annulus(radiusOuter=r_o, radiusInner=r_i, cellSize=cs, qdegree=3) + X = np.asarray(mesh.X.coords).copy() + R = np.sqrt((X ** 2).sum(1)); TH = np.arctan2(X[:, 1], X[:, 0]) + surf = np.abs(R - r_o) < 0.5 * cs + s_idx = np.where(surf)[0] + + # undeformed: boundary_normal(Upper) ≈ radial + n0 = _eval_boundary_normal(mesh, "Upper", X[s_idx]) + rhat0 = X[s_idx] / R[s_idx, None] + assert _angle_deg(n0, rhat0).max() < 3.0 + + # deform outer surface in mode-3 (12%) + Xd = X.copy() + Xd[surf] *= (1.0 + 0.12 * np.cos(3 * TH[surf]))[:, None] + mesh.deform(Xd, dt=1.0) + Xn = np.asarray(mesh.X.coords); Rn = np.sqrt((Xn ** 2).sum(1)) + rhat_d = Xn[s_idx] / Rn[s_idx, None] + + ref, good = _facet_vertex_normals(mesh, "Upper", _bvalue(mesh, "Upper")) + ref_s, good_s = ref[s_idx], good[s_idx] + tilt = _angle_deg(ref_s[good_s], rhat_d[good_s]) + assert tilt.max() > 8.0, "test setup: deformation should tilt the surface" + + nd = _eval_boundary_normal(mesh, "Upper", Xn[s_idx]) + err = _angle_deg(nd[good_s], ref_s[good_s]) + assert err.max() < 5.0, ( + f"boundary_normal(Upper) must track the deformed facet normal; " + f"max error {err.max():.1f}° (radial bug ~ {tilt.max():.1f}°)") + + +def test_boundary_normal_cartesian_free_surface_corner(): + """Cartesian free surface (box Top deformed, vertical side walls). The + Top normal must follow the deformed top and must NOT be averaged with the + side-wall normals at the corners.""" + cs = 1.0 / 8 + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cs, qdegree=3) + X = np.asarray(mesh.X.coords).copy() + top = np.abs(X[:, 1] - 1.0) < 1e-9 + t_idx = np.where(top)[0] + + # deform the top as a smooth bump h(x) = 0.1 sin(pi x); sides stay vertical + Xd = X.copy() + Xd[top, 1] += 0.1 * np.sin(np.pi * X[top, 0]) + mesh.deform(Xd, dt=1.0) + Xn = np.asarray(mesh.X.coords) + + ref, good = _facet_vertex_normals(mesh, "Top", _bvalue(mesh, "Top")) + ref_t, good_t = ref[t_idx], good[t_idx] + + nd = _eval_boundary_normal(mesh, "Top", Xn[t_idx]) + err = _angle_deg(nd[good_t], ref_t[good_t]) + assert err.max() < 5.0, ( + f"boundary_normal(Top) must track the deformed top; max {err.max():.1f}°") + + # the deformed top genuinely tilts (∂h/∂x ≠ 0 away from the crest) + yhat = np.tile([0.0, 1.0], (good_t.sum(), 1)) + assert _angle_deg(ref_t[good_t], yhat).max() > 8.0 + + # CORNER CHECK: the top-corner vertices (x=0, x=1) must carry the TOP + # normal, NOT a 45° average with the vertical walls. ∂h/∂x at x=0,1 is + # ±0.1π, so the true top normal there is tilted ~17° from vertical — and + # crucially has a SMALL... no: it must equal the one-sided top facet normal + # (the reference), i.e. err already bounds it. Assert the corner normal is + # not the 45° wall-average: + xc = Xn[t_idx, 0] + corner = (np.abs(xc) < 1e-9) | (np.abs(xc - 1.0) < 1e-9) + if corner.any(): + wall45 = np.tile([np.sqrt(0.5), np.sqrt(0.5)], (corner.sum(), 1)) + # corner top normal must be FAR from the (±0.707,0.707) wall-blend + assert _angle_deg(nd[corner], wall45).min() > 10.0, ( + "top corner normal looks averaged with the side wall (45°)") diff --git a/tests/test_0057_deformed_domain_membership.py b/tests/test_0057_deformed_domain_membership.py new file mode 100644 index 00000000..5fdba4f1 --- /dev/null +++ b/tests/test_0057_deformed_domain_membership.py @@ -0,0 +1,78 @@ +"""Domain membership + point location must use the DEFORMED geometry. + +Regression for the deformed-domain bug: points_in_domain() and +uw.function.evaluate() locate points using a boundary-skeleton kd-tree built +from mesh._nav_coords. That was captured as a reference to the ORIGINAL +coordinates in __init__ and never refreshed in nuke_coords_and_rebuild (only +the per-cell face-control-point arrays were invalidated, on adapt). So on a +volume mesh _nav_coords stayed at the undeformed boundary after a deform, and a +region that bulged OUT past the original boundary (r > r_o on a free surface) +read as EXTERIOR — stranding semi-Lagrangian trace-back feet there and +mis-locating evaluations to the old boundary, injecting the cold boundary value +at topographic highs. + +The fix invalidates the boundary kd-tree AND refreshes _nav_coords from the +current DM coordinates on every deform. + +Run: pixi run python -m pytest tests/test_0057_deformed_domain_membership.py -v +""" + +import pytest +import numpy as np +import sympy +import underworld3 as uw + +pytestmark = [ + pytest.mark.level_1, + pytest.mark.tier_a, + # Probes a rank-local bulge; the membership/eval logic is verified serial. + pytest.mark.skipif(uw.mpi.size > 1, reason="serial point-location check"), +] + + +def test_membership_and_eval_track_bulged_surface(): + r_i, r_o, cs = 0.5, 1.0, 0.15 + mesh = uw.meshing.Annulus(radiusOuter=r_o, radiusInner=r_i, cellSize=cs, qdegree=3) + + # T = radius (a field whose value identifies where a point actually is) + T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) + r_sym = sympy.sqrt(mesh.X[0] ** 2 + mesh.X[1] ** 2) + T.data[:, 0] = np.asarray(uw.function.evaluate(r_sym, T.coords)).flatten() + + # Bulge the outer surface OUTWARD in mode-3 (only the outward half), so the + # surface near theta=0 sits well beyond the original radius r_o. + X = np.asarray(mesh.X.coords).copy() + R = np.sqrt((X ** 2).sum(1)); TH = np.arctan2(X[:, 1], X[:, 0]) + surf = np.abs(R - r_o) < 0.5 * cs + Xd = X.copy() + Xd[surf] *= (1.0 + 0.12 * np.maximum(np.cos(3 * TH[surf]), 0.0))[:, None] + mesh.deform(Xd, dt=1.0) + + Rn = np.sqrt((np.asarray(mesh.X.coords) ** 2).sum(1)) + surf_now = np.abs(np.asarray(mesh.X.coords)[:, 1]) < 0.05 # near theta=0 + crest_r = Rn[surf & (np.abs(TH) < 0.2)].max() if (surf & (np.abs(TH) < 0.2)).any() else Rn.max() + assert crest_r > r_o + 0.05, "test setup: surface should bulge past r_o" + + # nav coords must now reflect the deformed geometry + nav_max = np.sqrt((np.asarray(mesh._nav_coords) ** 2).sum(1)).max() + assert nav_max > r_o + 0.05, ( + f"_nav_coords stale: max radius {nav_max:.3f} (deformed ~ {crest_r:.3f})") + + # probe points along the bulge crest (theta=0), from inside r_o out to the + # deformed surface. ALL must be inside the (deformed) domain. + radii = np.linspace(0.95, crest_r - 0.01, 10) + pts = np.column_stack([radii, np.zeros_like(radii)]) + + inside = mesh.points_in_domain(pts, strict_validation=True) + assert inside.all(), ( + f"points in the bulge (r in [{radii[0]:.2f},{radii[-1]:.2f}], surface " + f"at {crest_r:.2f}) wrongly flagged exterior: {inside}") + + # evaluate must LOCATE them (not clamp to the old boundary). T is carried + # Lagrangian-ly, so the value need not equal r; but it must NOT all collapse + # to ~r_o (the cold-clamp signature), and must vary smoothly with radius. + Tev = np.asarray(uw.function.evaluate(T.sym[0], pts)).flatten() + assert np.ptp(Tev) > 0.02, ( + f"evaluate collapsed to a single boundary value (cold-clamp): {Tev}") + assert np.all(np.diff(Tev) >= -1e-6), ( + f"evaluate not monotone along the carried-T crest: {Tev}") From e086e72a9e69b6b06263739b704b31e1da273ccb Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 22 Jun 2026 09:45:01 +1000 Subject: [PATCH 2/2] Address Copilot review on #264 - boundary_normal(): drop unused cKDTree import (tree is built in _assemble_boundary_normal). - _assemble_boundary_normal(): drop unused `import underworld3 as uw`; count P1 vertices PER FACET from its closure (vStart..vEnd) instead of assuming nverts == cdim, so non-simplex facets (3D quad/hex boundaries) are handled correctly. - add_nitsche_bc docstring: `normal` default is now the per-boundary, deformation-tracking mesh.boundary_normal(boundary), not mesh.Gamma_N. - tests: drop unused surf_now; tidy an editorial comment fragment. test_0056 + test_0057 pass (3 passed). Underworld development team with AI support from Claude Code --- src/underworld3/cython/petsc_generic_snes_solvers.pyx | 4 ++-- src/underworld3/discretisation/discretisation_mesh.py | 8 +++++--- tests/test_0056_projected_normals_deform.py | 7 +++---- tests/test_0057_deformed_domain_membership.py | 1 - 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 29c5849e..1db6079e 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -4723,8 +4723,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): a fault orientation field). normal : sympy.Matrix or list, optional Boundary unit normal used in the Nitsche consistency, symmetry, - and pressure-coupling terms. Default ``None`` uses the PETSc - boundary facet normal ``mesh.Gamma_N``. + and pressure-coupling terms. Default ``None`` uses the per-boundary, + deformation-tracking ``mesh.boundary_normal(boundary)``. gamma : float, default=10.0 Dimensionless stabilisation parameter. Typical values 5--20 for P2 elements. diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 585e1798..32b90128 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2507,7 +2507,6 @@ def boundary_normal(self, boundary): Boundary label name (or a ``mesh.boundaries`` enum member). """ import underworld3 as uw - from scipy.spatial import cKDTree name = getattr(boundary, "name", str(boundary)) if not hasattr(self, "_boundary_normal_vars") or self._boundary_normal_vars is None: @@ -2524,7 +2523,6 @@ def boundary_normal(self, boundary): def _assemble_boundary_normal(self, var, name): """Fill ``var`` with the area-weighted outward facet normal assembled from the faces of boundary ``name`` only (see :meth:`boundary_normal`).""" - import underworld3 as uw from scipy.spatial import cKDTree cdim = self.cdim dm = self.dm @@ -2558,7 +2556,9 @@ def _assemble_boundary_normal(self, var, name): face_pts.append(int(p)) tree = cKDTree(coords) - nverts = cdim # P1 vertices per facet (2D edge=2, 3D tri=3) + # P1 vertices per facet, counted from the facet's own closure so this + # works for non-simplex facets too (2D edge=2, 3D tri=3, 3D quad=4). + vStart, vEnd = dm.getDepthStratum(0) for f in face_pts: if dm.getSupportSize(f) != 1: continue @@ -2569,6 +2569,8 @@ def _assemble_boundary_normal(self, var, name): if numpy.dot(nrm, numpy.asarray(cent)[:cdim] - numpy.asarray(ccent)[:cdim]) < 0: nrm = -nrm + _clo = dm.getTransitiveClosure(f)[0] + nverts = int(numpy.count_nonzero((_clo >= vStart) & (_clo < vEnd))) or cdim # Accumulate to the facet's P1 DOFs (its vertices) — found as the # nearest DOFs to the facet centroid. This avoids indexing the local # coordinate array by (vertex_point - vStart), which is only valid diff --git a/tests/test_0056_projected_normals_deform.py b/tests/test_0056_projected_normals_deform.py index 53f51222..056999ce 100644 --- a/tests/test_0056_projected_normals_deform.py +++ b/tests/test_0056_projected_normals_deform.py @@ -148,10 +148,9 @@ def test_boundary_normal_cartesian_free_surface_corner(): # CORNER CHECK: the top-corner vertices (x=0, x=1) must carry the TOP # normal, NOT a 45° average with the vertical walls. ∂h/∂x at x=0,1 is - # ±0.1π, so the true top normal there is tilted ~17° from vertical — and - # crucially has a SMALL... no: it must equal the one-sided top facet normal - # (the reference), i.e. err already bounds it. Assert the corner normal is - # not the 45° wall-average: + # ±0.1π, so the true top normal there is tilted ~17° from vertical; it must + # equal the one-sided top facet normal (the reference), which `err` already + # bounds. Assert the corner normal is not the 45° wall-average: xc = Xn[t_idx, 0] corner = (np.abs(xc) < 1e-9) | (np.abs(xc - 1.0) < 1e-9) if corner.any(): diff --git a/tests/test_0057_deformed_domain_membership.py b/tests/test_0057_deformed_domain_membership.py index 5fdba4f1..d17a42c1 100644 --- a/tests/test_0057_deformed_domain_membership.py +++ b/tests/test_0057_deformed_domain_membership.py @@ -49,7 +49,6 @@ def test_membership_and_eval_track_bulged_surface(): mesh.deform(Xd, dt=1.0) Rn = np.sqrt((np.asarray(mesh.X.coords) ** 2).sum(1)) - surf_now = np.abs(np.asarray(mesh.X.coords)[:, 1]) < 0.05 # near theta=0 crest_r = Rn[surf & (np.abs(TH) < 0.2)].max() if (surf & (np.abs(TH) < 0.2)).any() else Rn.max() assert crest_r > r_o + 0.05, "test setup: surface should bulge past r_o"