From 89bb1a1cf340aec711feb6dcf1b1b524547413ec Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Fri, 10 Jul 2026 07:22:29 +1000 Subject: [PATCH 1/4] fix(mesh): add mesh.canonical_normal accessor for partition-safe internal-boundary normals (#327) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `mesh.Gamma[i]` JIT-substitutes `petsc_n[i]`, PETSc's per-quadrature outward-normal-from-`support[0]` in `DMPlexComputeBdIntegral`. On an internal boundary at a partition seam, different ranks disagree on which cell is `support[0]` for the one seam facet, so the outward normal flips sign there and a signed integral like `∫ n_y dS` is silently off by O(seam facets / total facets) in parallel — the exact `1 − 2/32 = 0.9375` signature #327 reports. The mesh factory Enums (`BoxInternalBoundary.boundary_normals_{2,3}D`, `AnnulusInternalBoundary.boundary_normals`, `SphericalShellInternalBoundary.boundary_normals`) already declare the analytic outward-pointing normal for each boundary, but nothing consumed them. Wire that up: `Mesh.canonical_normal(boundary_name)` returns the declared sympy Matrix (or None if not declared), giving callers a partition-independent alternative to `Gamma` for meshes that know their internal-boundary geometry analytically. The `test_0502` "internal_normal_ny" / "internal_normal_weighted" tests are un-skipped (they were serial-only with a TODO(BUG) marker citing #327) and now use the canonical accessor. A new `test_bd_integral_internal_canonical_normal_off_grid_zint` regression constructs the failing `zint=0.55` configuration that reliably picks a partition seam through the internal boundary at np=2 (from the sweep in the /remote-control session's #327 investigation). 23/23 pass in `tests/test_0502_boundary_integrals.py` at both np=1 and np=2. Not fixed: `Gamma[i]` itself remains partition-dependent on internal boundaries — users writing code that reads `Gamma[i]` on an internal boundary still get the wrong answer at np>1 unless they switch to `canonical_normal`. The full fix (canonicalise `support[0]` at construction, or expose `DMPlexOrientLabel` via a Cython wrapper) is larger and out of scope for this PR — see the #327 investigation notes. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 55 ++++++++++++++++ tests/test_0502_boundary_integrals.py | 63 ++++++++++++------- 2 files changed, 95 insertions(+), 23 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 8859464a..1c499d2b 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2709,6 +2709,61 @@ def Gamma_P1(self): self._update_projected_normals() return self._projected_normals.sym + def canonical_normal(self, boundary_name): + r"""Analytic outward-pointing normal for a boundary, or ``None`` + if no analytic normal was declared for that boundary. + + Sourced from the mesh factory's ``boundary_normals`` Enum: for + axis-aligned box boundaries this is a constant sympy Matrix, for + annulus / spherical-shell radial boundaries it is the analytic + radial unit vector, and so on. + + The primary caller is code that needs a **partition-safe** normal + on an *internal* boundary — see :issue:`327`. On an internal + boundary at a partition seam, PETSc's per-quadrature ``petsc_n[]`` + (surfacing as :attr:`Gamma` / :attr:`Gamma_N`) is derived from + ``support[0]`` of the DMPlex facet closure, which is + partition-dependent; different ranks disagree on which cell is + "support[0]" for the one seam facet, and the outward normal of + that facet flips sign. A signed integral of ``Gamma[k]`` is then + wrong by O(seam-facets / total-facets). The analytic normal + returned here is partition-independent and does not touch + ``petsc_n``, so it sidesteps the defect entirely for the mesh + classes that know their internal-boundary geometry + (:class:`BoxInternalBoundary`, :class:`AnnulusInternalBoundary`, + :class:`SphericalShellInternalBoundary`). + + Parameters + ---------- + boundary_name : str + Name of the boundary label to look up (case-sensitive; must + match one of :attr:`boundaries`). + + Returns + ------- + sympy.Matrix or None + Row matrix of length :attr:`cdim` giving the outward-pointing + normal, or ``None`` if this mesh factory did not declare + an analytic normal for ``boundary_name``. + + See Also + -------- + Gamma : the raw per-quadrature normal — use for external + boundaries; may be partition-dependent on internal seams. + Gamma_N : normalised :attr:`Gamma`. + Gamma_P1 : projected P1 normals, useful for curved external + boundaries. + """ + bn = getattr(self, "boundary_normals", None) + if bn is None: + return None + try: + member = bn[boundary_name] + except (KeyError, AttributeError): + return None + value = getattr(member, "value", member) + return value + # =================================================================== # Bounding surfaces — per-surface tangent-slip + restore. # See docs/developer/design/boundary-slip-strategy.md. SEPARATE from diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py index 5f7d323f..a67e040f 100644 --- a/tests/test_0502_boundary_integrals.py +++ b/tests/test_0502_boundary_integrals.py @@ -159,24 +159,22 @@ def test_bd_integral_internal_coordinate_fn(): assert abs(value - 0.5) < 0.01, f"Expected 0.5, got {value}" -# TODO(BUG): internal-boundary facet-normal orientation is rank-dependent at -# partition seams: at np2 one seam facet contributes with flipped sign, so -# |integral of n_y| = 1 - 2/32 = 0.9375. Scalar integrands (length, -# coordinate functions) are exact in parallel; only signed-normal integrands -# are affected. Found while unskipping after the BF-13 constructor fix -# (2026-07 audit) — separate defect, not covered by BF-13. -@pytest.mark.skipif( - uw.mpi.size > 1, - reason="Internal-boundary normal orientation is rank-dependent at partition seams (see TODO(BUG) above)", -) +# The `mesh.Gamma` normal on an *internal* boundary is derived from petsc_n[] +# which uses DMPlex support[0] — partition-dependent at seam facets, so a +# signed-normal integral like ∫ n_y dS is silently off by O(seam facets / +# total facets) in parallel (issue #327). These tests use the analytic +# `canonical_normal` accessor, which returns the mesh factory's declared +# outward normal for the boundary (a partition-independent sympy expression); +# see the accessor's docstring in discretisation_mesh.py. def test_bd_integral_internal_normal_ny(): - """Integrate n_y along internal boundary at y=0.5. - The internal boundary has normals pointing in +y or -y direction, - so integrating n_y should give +1 or -1 (length 1 boundary).""" + """Integrate n_y along internal boundary at y=0.5 using the analytic + canonical normal. The internal boundary has normals pointing in +y or + -y direction, so integrating n_y should give +1 or -1 (length 1 + boundary).""" mesh_internal, _, _ = _get_internal_mesh() - Gamma = mesh_internal.Gamma - n_y = Gamma[1] + normal = mesh_internal.canonical_normal("Internal") + n_y = normal[1] bd_int = uw.maths.BdIntegral(mesh_internal, fn=n_y, boundary="Internal") value = bd_int.evaluate() @@ -200,17 +198,14 @@ def test_bd_integral_internal_normal_nx(): assert abs(value) < 0.01, f"Expected ~0, got {value}" -@pytest.mark.skipif( - uw.mpi.size > 1, - reason="Internal-boundary normal orientation is rank-dependent at partition seams (see TODO(BUG) above)", -) def test_bd_integral_internal_normal_weighted(): - """Integrate x * n_y along internal boundary at y=0.5. - int_0^1 x * n_y dx = n_y * 0.5. Since |n_y| = 1, result should be ~0.5.""" + """Integrate x * n_y along internal boundary at y=0.5 using the analytic + canonical normal. int_0^1 x * n_y dx = n_y * 0.5. Since |n_y| = 1, + result should be ~0.5.""" mesh_internal, x_i, _ = _get_internal_mesh() - Gamma = mesh_internal.Gamma - n_y = Gamma[1] + normal = mesh_internal.canonical_normal("Internal") + n_y = normal[1] bd_int = uw.maths.BdIntegral(mesh_internal, fn=x_i * n_y, boundary="Internal") value = bd_int.evaluate() @@ -218,6 +213,28 @@ def test_bd_integral_internal_normal_weighted(): assert abs(abs(value) - 0.5) < 0.01, f"Expected |value| = 0.5, got {value}" +def test_bd_integral_internal_canonical_normal_off_grid_zint(): + """Regression for the failing partition-through-boundary case from #327. + + With ``zintCoord=0.55`` (off-grid), the mpirun -n 2 partition seam runs + through the internal boundary and one seam facet's ``petsc_n[]`` flips + sign: ``∫ Gamma[1] dS`` returns 0.9375 = 1 − 2/32 instead of 1.0. The + analytic ``canonical_normal("Internal")`` returned by the mesh factory + does not touch ``petsc_n[]`` and gives the exact value regardless of + partition.""" + mesh_off = BoxInternalBoundary( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=1.0/32.0, zintCoord=0.55, simplex=True, + ) + # Need at least one variable so BdIntegral has a section to integrate against + uw.discretisation.MeshVariable("T_off", mesh_off, 1, degree=2) + + n_y = mesh_off.canonical_normal("Internal")[1] + val = uw.maths.BdIntegral(mesh_off, fn=n_y, boundary="Internal").evaluate() + assert abs(abs(val) - 1.0) < 1e-6, ( + f"canonical_normal internal integral should be exactly ±1, got {val}") + + def test_bd_integral_internal_does_not_affect_external(): """External boundaries should still work on the internal-boundary mesh.""" From 2829c9ac57464762b1a5e78c8e4373408a2cee75 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 14 Jul 2026 10:03:52 +1000 Subject: [PATCH 2/4] fix(adaptivity): deadlock in _dm_unstack_bcs at np>=4 (collective mismatch) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit getNonEmptyStratumValuesIS() is rank-local, but the per-value loop makes collective calls (labelComplete). When a rank's partition touches no facet of some boundary — routine at np>=4 on BoxInternalBoundary — ranks make different numbers of collective calls and the run deadlocks in mesh construction (stacks: one rank in labelComplete, others in markBoundaryFaces). np=2 escaped because both halves touch every boundary. Iterate the allgathered union of stratum values so every rank makes the same sequence of collective calls; guard the local setStratumIS to values live on this rank. Verified: BoxInternalBoundary construction and the full 0502 suite now pass at np=4 (24/24, previously hung indefinitely). Underworld development team with AI support from Claude Code --- src/underworld3/adaptivity.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/underworld3/adaptivity.py b/src/underworld3/adaptivity.py index c22a1fb8..e25282ae 100644 --- a/src/underworld3/adaptivity.py +++ b/src/underworld3/adaptivity.py @@ -502,6 +502,17 @@ def _dm_unstack_bcs(dm, boundaries, stacked_bc_label_name): stacked_bc_label = dm.getLabel(stacked_bc_label_name) vals = stacked_bc_label.getNonEmptyStratumValuesIS().getIndices() + # BUGFIX: ``vals`` is rank-local (a rank only sees stratum values that + # are non-empty on its own partition), but the loop below makes + # collective calls (``labelComplete``) once per value. If ranks + # iterate different value sets — routine at np>=4, where a partition + # quadrant may touch no facet of some boundary — the collective call + # counts diverge and the run deadlocks. Iterate the global union so + # every rank makes the same sequence of collective calls; guard the + # local stratum access to values that are live on this rank. + local_vals = set(int(v) for v in vals) + vals = sorted(set().union(*uw.mpi.comm.allgather(local_vals))) + # Clear labels just in case for b in boundaries: dm.removeLabel(b.name) @@ -520,8 +531,9 @@ def _dm_unstack_bcs(dm, boundaries, stacked_bc_label_name): continue b_dmlabel = dm.getLabel(b.name) - lab_is = stacked_bc_label.getStratumIS(v) - b_dmlabel.setStratumIS(v, lab_is) + if v in local_vals: + lab_is = stacked_bc_label.getStratumIS(v) + b_dmlabel.setStratumIS(v, lab_is) # BUGFIX(#162): expand the boundary label to include its closure # (the vertices and edges bordering the labeled facets). The From 107e9f3e089db6e2dce6eb246eb308764d6a7689 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 14 Jul 2026 10:04:13 +1000 Subject: [PATCH 3/4] fix(mesh): unify boundary-normal interface behind mesh.Gamma (#327) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit mesh.Gamma is now the single user-facing boundary-normal symbol. The consumers that know their boundary (uw.maths.BdIntegral, natural BCs via add_condition) resolve it per boundary through the new Mesh._resolve_boundary_normals: - external boundaries: unchanged — Gamma compiles to PETSc's exact per-quadrature outward unit normal petsc_n[] (zero solver risk); - internal boundaries: petsc_n[] orientation follows the partition-dependent DMPlex support[0] (issue #327), so Gamma components are substituted with the mesh factory's declared analytic normal (canonical_normal), partition-independent by construction; - internal boundary with no declared normal: loud RuntimeError — there is no orientation convention for an arbitrary internal surface until DMPlexOrientLabel supports open surfaces (pending upstream). Internal-boundary detection (_boundary_is_internal): any labeled facet with support size 2, MAX-allreduced and cached per boundary. Supporting changes: - Convention: declared internal normals point from region Inner to region Outer. 2D BoxInternalBoundary Internal flipped [0,-1] -> [0,1] to match 3D (+z) and the radial (unit_e_0) annulus/spherical convention. - AnnulusInternalBoundary: attach the existing boundary_normals Enum (was defined but never attached); SphericalShellInternalBoundary: add and attach one. TODO(BUG) noted on Box external entries, which are inward-pointing. - Gamma_N docstring: deprecated alias of Gamma (numerically identical). Gamma_P1 docstring: deprecated, with diagnosis (off-kernel petsc_n fallback + sign cancellation on internal facets -> sub-unit vectors). - SNES_Stokes_Constrained.add_constraint_bc default normal now boundary_normal(boundary) (oriented-then-averaged per-boundary field) instead of the deprecated global Gamma_P1. - tests(0502): internal-normal tests use plain mesh.Gamma with tight 1e-6 assertions and the Inner->Outer (+y) sign; back-compat test for the canonical_normal accessor; #327 off-grid zint regression kept. Verified: 0502 suite 24/24 at np=1, np=2, np=4; natural-BC/Nitsche/ constrained-freeslip regression trio 23 passed. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 7 + src/underworld3/cython/petsc_maths.pyx | 14 +- .../discretisation/discretisation_mesh.py | 153 +++++++++++++++--- src/underworld3/meshing/annulus.py | 4 +- src/underworld3/meshing/cartesian.py | 8 +- src/underworld3/meshing/spherical.py | 9 +- src/underworld3/systems/solvers.py | 8 +- tests/test_0502_boundary_integrals.py | 74 +++++---- 8 files changed, 212 insertions(+), 65 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 2d6a7f84..fcb447b9 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1600,6 +1600,13 @@ class SolverBaseClass(uw_object): from collections import namedtuple if c_type == 'neumann': + # mesh.Gamma in a natural-BC expression resolves per boundary: + # external boundaries keep the exact per-quadrature petsc_n[]; + # internal boundaries substitute the declared analytic normal + # (petsc_n is orientation-ambiguous there — issue #327). + sympy_fn = sympy.Matrix( + self.mesh._resolve_boundary_normals(sympy_fn, label) + ).as_immutable() BC = namedtuple('NaturalBC', ['f_id', 'components', 'fn_f', 'fn_F', 'fn_p', 'boundary', 'boundary_label_val', 'type', 'PETScID', 'fns']) self.natural_bcs.append(BC(f_id, components, sympy_fn, None, None, label, -1, "natural", -1, {})) elif c_type == 'dirichlet': diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx index effaa951..5be87ac7 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/src/underworld3/cython/petsc_maths.pyx @@ -346,8 +346,10 @@ class BdIntegral: for some scalar function :math:`f` over a named boundary :math:`\partial \Omega` of the mesh. In 2D this is a line integral; in 3D a surface integral. - The integrand may reference the outward unit normal via ``mesh.Gamma_N`` - (components map to ``petsc_n[0], petsc_n[1], ...`` in the generated C code). + The integrand may reference the outward unit normal via ``mesh.Gamma`` + (on external boundaries the components map to ``petsc_n[0], petsc_n[1], + ...`` in the generated C code; on internal boundaries they resolve to the + mesh's declared analytic normal — see ``Mesh._resolve_boundary_normals``). Parameters ---------- @@ -421,9 +423,15 @@ class BdIntegral: mesh = self.mesh + # mesh.Gamma is resolved per boundary: external boundaries keep the + # exact per-quadrature petsc_n[]; internal boundaries substitute the + # declared analytic normal (petsc_n is orientation-ambiguous there, + # see Mesh._resolve_boundary_normals / issue #327). + fn = mesh._resolve_boundary_normals(self.fn, self.boundary) + # Compile integrand using the boundary residual slot (includes petsc_n[] in signature) _getext_result = getext( - self.mesh, JITCallbackSet(bd_residual=(self.fn,)), + self.mesh, JITCallbackSet(bd_residual=(fn,)), self.mesh.vars.values(), verbose=verbose) cdef PtrContainer ext = _getext_result.ptrobj diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 1c499d2b..3ee2d31d 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2534,17 +2534,14 @@ def boundary_normal(self, boundary): 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`).""" - from scipy.spatial import cKDTree - cdim = self.cdim - dm = self.dm - coords = numpy.ascontiguousarray(var.coords) - accum = numpy.zeros((coords.shape[0], cdim)) + def _boundary_facets(self, name): + """This rank's facet (height-1) points carrying boundary label ``name``. - # faces carrying this boundary label: DM label named after the boundary, - # stratum keyed by the boundary's value (same access the BC code uses). + Uses the DM label named after the boundary with the stratum keyed by + the boundary's value (the same access the BC code uses). Returns an + empty list when this rank owns no facets of the boundary. + """ + dm = self.dm bvalue = None for b in (self.boundaries or []): if b.name == name: @@ -2560,6 +2557,8 @@ def _assemble_boundary_normal(self, var, name): vis = label.getValueIS() live = set(int(x) for x in vis.getIndices()) if vis.getSize() else set() except Exception: + # Sanctioned: no value IS on this rank means no facets here — + # fall through to the empty list. live = set() if int(bvalue) in live: pis = label.getStratumIS(bvalue) @@ -2568,6 +2567,18 @@ def _assemble_boundary_normal(self, var, name): for p in pis.getIndices(): if fS <= int(p) < fE: face_pts.append(int(p)) + return face_pts + + 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`).""" + from scipy.spatial import cKDTree + cdim = self.cdim + dm = self.dm + coords = numpy.ascontiguousarray(var.coords) + accum = numpy.zeros((coords.shape[0], cdim)) + + face_pts = self._boundary_facets(name) tree = cKDTree(coords) # P1 vertices per facet, counted from the facet's own closure so this @@ -2696,12 +2707,15 @@ def _assemble_cell_size(self, var): @property def Gamma_P1(self): - """Projected P1 boundary normals as a sympy Matrix. + """Deprecated — use :attr:`Gamma` in integrands and BCs, or + :meth:`boundary_normal` for a per-boundary P1 normal field. - Returns the normalised, vertex-averaged PETSc face normals - as a smooth P1 field. Preferred over :attr:`Gamma_N` for - penalty and Nitsche BCs on curved boundaries — gives - consistent orientation and better convergence in 3D. + This global field point-evaluates :attr:`Gamma` at every mesh + vertex, but the underlying ``petsc_n`` only exists inside + surface-integral kernels; off-kernel evaluation falls back to a + coordinate-based direction and, on internal boundaries, averages + oppositely-oriented facet normals into sub-unit vectors. Retained + unchanged for back-compat (mesh-smoothing internals). Automatically updated when the mesh deforms. """ @@ -2764,6 +2778,90 @@ def canonical_normal(self, boundary_name): value = getattr(member, "value", member) return value + def _boundary_is_internal(self, boundary_name): + """True if boundary ``boundary_name`` is an internal surface — its + facets have a cell on BOTH sides (support size 2) — rather than part + of the mesh exterior (support size 1). + + Collective on first call per boundary (the local answer is + MAX-reduced so every rank agrees even when it owns no facets of the + boundary); cached thereafter. + """ + import underworld3 as uw + from mpi4py import MPI + + if not hasattr(self, "_internal_boundary_cache") or \ + self._internal_boundary_cache is None: + self._internal_boundary_cache = {} + cached = self._internal_boundary_cache.get(boundary_name) + if cached is not None: + return cached + + dm = self.dm + local = 0 + for f in self._boundary_facets(boundary_name): + if dm.getSupportSize(f) == 2: + local = 1 + break + result = bool(uw.mpi.comm.allreduce(local, op=MPI.MAX)) + self._internal_boundary_cache[boundary_name] = result + return result + + def _resolve_boundary_normals(self, fn, boundary_name): + r"""Return ``fn`` with :attr:`Gamma` resolved for ``boundary_name``. + + ``mesh.Gamma`` is the single user-facing boundary-normal symbol. On + an EXTERNAL boundary its components compile directly to PETSc's + exact per-quadrature outward unit normal (``petsc_n[]``) and ``fn`` + is returned unchanged. On an INTERNAL boundary ``petsc_n[]`` is + orientation-ambiguous — the facet has two support cells and the sign + follows the partition-dependent ``support[0]`` — so a signed integral + of ``Gamma`` components is partition-dependent (issue #327). Here the + ``Gamma`` components are substituted with the mesh factory's declared + analytic normal (:meth:`canonical_normal`), which is + partition-independent by construction. + + Called by every boundary-integrand consumer that knows its boundary + (``uw.maths.BdIntegral``, natural boundary conditions), so users + write ``mesh.Gamma`` everywhere and never choose an implementation. + + Raises + ------ + RuntimeError + If ``boundary_name`` is internal and the mesh declares no + analytic normal for it. There is no orientation convention for + an arbitrary internal surface until the surface itself is + oriented (PETSc ``DMPlexOrientLabel``, open-surface support + pending upstream), so failing loudly beats integrating a + partition-dependent sign. + """ + import underworld3 as uw + + gamma_syms = tuple(self._Gamma.base_scalars()[0:self.cdim]) + expr = uw.function.expressions.unwrap( + fn, keep_constants=False, return_self=False) + if not isinstance(expr, sympy.MatrixBase): + expr = sympy.sympify(expr) + free = expr.free_symbols + if not any(g in free for g in gamma_syms): + return fn + if not self._boundary_is_internal(boundary_name): + return fn + + normal = self.canonical_normal(boundary_name) + if normal is None: + raise RuntimeError( + f"The integrand references mesh.Gamma on the internal " + f"boundary '{boundary_name}', but this mesh declares no " + f"analytic normal for it. On an internal surface the " + f"discrete facet normal has no well-defined orientation " + f"(issue #327). Declare one in the mesh factory's " + f"'boundary_normals' Enum, or use an explicit normal " + f"expression in place of mesh.Gamma." + ) + subs = {g: normal[i] for i, g in enumerate(gamma_syms)} + return expr.xreplace(subs) + # =================================================================== # Bounding surfaces — per-surface tangent-slip + restore. # See docs/developer/design/boundary-slip-strategy.md. SEPARATE from @@ -3508,11 +3606,11 @@ def N(self) -> sympy.vector.CoordSys3D: @property def Gamma_N(self) -> sympy.Matrix: - r"""Normalised boundary/surface normal as a row matrix. + r"""Deprecated alias — use :attr:`Gamma`. - Returns ``Gamma / |Gamma|`` so that the result is a unit normal - regardless of element size. Use this for penalty and Nitsche BCs - where mesh-independent scaling is required. + Retained for back-compatibility: returns ``Gamma / |Gamma|``, which + is numerically identical to :attr:`Gamma` (the quadrature-point + normal is already unit length). Returns ------- @@ -3524,15 +3622,24 @@ def Gamma_N(self) -> sympy.Matrix: @property def Gamma(self) -> sympy.Matrix: - r"""Raw (un-normalised) boundary coordinate scalars as a row matrix. + r"""The boundary unit normal, as a symbolic row matrix. + + This is the single user-facing normal symbol: use it in boundary + integrands (:class:`uw.maths.BdIntegral`) and natural boundary + conditions on any boundary, external or internal. - The magnitude scales with face edge length (2D) or face area (3D). - For a unit normal, use :attr:`Gamma_N` instead. + The consumer that compiles the integrand resolves the symbol per + boundary: on an external boundary the components map to PETSc's + exact per-quadrature outward unit normal (``petsc_n[]``); on an + internal boundary — where the discrete facet normal has no + well-defined orientation (issue #327) — they are substituted with + the mesh factory's declared analytic normal (see + :meth:`canonical_normal` and :meth:`_resolve_boundary_normals`). Returns ------- sympy.Matrix - Row matrix of boundary coordinate scalars. + Row matrix of boundary normal components. """ return sympy.Matrix(self._Gamma.base_scalars()[0 : self.cdim]).T diff --git a/src/underworld3/meshing/annulus.py b/src/underworld3/meshing/annulus.py index a1f47080..94e7ba92 100644 --- a/src/underworld3/meshing/annulus.py +++ b/src/underworld3/meshing/annulus.py @@ -1414,7 +1414,9 @@ class boundary_normals(Enum): Internal = new_mesh.CoordinateSystem.unit_e_0 Centre = None - # boundary_normals deprecated — use mesh.Gamma_P1 for boundary normals + # Consumed by Mesh.canonical_normal — mesh.Gamma on the Internal + # boundary resolves to this analytic radial normal (issue #327). + new_mesh.boundary_normals = boundary_normals new_mesh.regions = regions # Full annulus with internal boundary: rigid rotation about z-axis diff --git a/src/underworld3/meshing/cartesian.py b/src/underworld3/meshing/cartesian.py index 806cd495..4f965976 100644 --- a/src/underworld3/meshing/cartesian.py +++ b/src/underworld3/meshing/cartesian.py @@ -492,12 +492,18 @@ class regions_2D(Enum): Inner = 101 # Below internal boundary Outer = 102 # Above internal boundary + # TODO(BUG): the external entries below are INWARD-pointing while + # Mesh.canonical_normal documents the declared normal as outward — + # audit consumers before flipping them. The Internal entry follows the + # internal-boundary convention: it points from region Inner to region + # Outer (+y here; +z in 3D; radially outward, unit_e_0, on annulus and + # spherical shells). class boundary_normals_2D(Enum): Bottom = sympy.Matrix([0, 1]) Top = sympy.Matrix([0, -1]) Right = sympy.Matrix([-1, 0]) Left = sympy.Matrix([1, 0]) - Internal = sympy.Matrix([0, -1]) + Internal = sympy.Matrix([0, 1]) class boundaries_3D(Enum): Bottom = 11 diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py index 36cbd6ff..faaea74d 100644 --- a/src/underworld3/meshing/spherical.py +++ b/src/underworld3/meshing/spherical.py @@ -746,7 +746,14 @@ def spherical_mesh_refinement_callback(dm): verbose=verbose, ) - # boundary_normals deprecated — use mesh.Gamma_P1 for boundary normals + class boundary_normals(Enum): + Lower = new_mesh.CoordinateSystem.unit_e_0 + Internal = new_mesh.CoordinateSystem.unit_e_0 + Upper = new_mesh.CoordinateSystem.unit_e_0 + + # Consumed by Mesh.canonical_normal — mesh.Gamma on the Internal + # boundary resolves to this analytic radial normal (issue #327). + new_mesh.boundary_normals = boundary_normals new_mesh.regions = regions diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 01bbb144..27f6877f 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -2371,7 +2371,8 @@ def add_constraint_bc(self, conds=None, boundary=None, normal=None, screening=No boundary : str Mesh boundary label (e.g. ``"Upper"``). normal : sympy matrix, optional - Row-vector constraint normal. Defaults to ``mesh.Gamma_P1``. + Row-vector constraint normal. Defaults to + ``mesh.boundary_normal(boundary)``. screening : float or sympy expression, optional Interior screening coefficient :math:`\varepsilon` (de-singularises the interior multiplier DOFs). Defaults to ``1e-6``. @@ -2429,7 +2430,10 @@ def add_constraint_bc(self, conds=None, boundary=None, normal=None, screening=No ) if normal is None: - normal = self.mesh.Gamma_P1 + # Per-boundary assembled normal (facet normals oriented before + # averaging) — the legacy global mesh.Gamma_P1 point-evaluates + # petsc_n off-kernel and is deprecated (see Mesh.Gamma_P1). + normal = self.mesh.boundary_normal(boundary) normal = sympy.Matrix(normal) if normal.shape[0] != 1: normal = normal.reshape(1, self.mesh.dim) diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py index a67e040f..b32cca7a 100644 --- a/tests/test_0502_boundary_integrals.py +++ b/tests/test_0502_boundary_integrals.py @@ -113,9 +113,9 @@ def test_bd_integral_invalid_boundary(): # --- Internal boundary tests (BoxInternalBoundary) --- # These run in serial and parallel: the BoxInternalBoundary rank>0 -# UnboundLocalError (2026-07 audit, BF-13) is fixed. Two signed-normal -# tests remain serial-only — see the TODO(BUG) on -# test_bd_integral_internal_normal_ny. +# UnboundLocalError (2026-07 audit, BF-13) is fixed, and signed-normal +# integrands written with plain mesh.Gamma are partition-safe (resolved +# to the declared analytic normal — issue #327). from underworld3.meshing import BoxInternalBoundary @@ -159,29 +159,24 @@ def test_bd_integral_internal_coordinate_fn(): assert abs(value - 0.5) < 0.01, f"Expected 0.5, got {value}" -# The `mesh.Gamma` normal on an *internal* boundary is derived from petsc_n[] -# which uses DMPlex support[0] — partition-dependent at seam facets, so a -# signed-normal integral like ∫ n_y dS is silently off by O(seam facets / -# total facets) in parallel (issue #327). These tests use the analytic -# `canonical_normal` accessor, which returns the mesh factory's declared -# outward normal for the boundary (a partition-independent sympy expression); -# see the accessor's docstring in discretisation_mesh.py. +# `mesh.Gamma` is the single user-facing normal symbol on any boundary. +# On an internal boundary the raw petsc_n[] is orientation-ambiguous +# (DMPlex support[0] is partition-dependent at seam facets — issue #327), +# so BdIntegral resolves the Gamma components to the mesh factory's +# declared analytic normal (Mesh._resolve_boundary_normals). The declared +# internal normal points from region Inner to region Outer (+y here). def test_bd_integral_internal_normal_ny(): - """Integrate n_y along internal boundary at y=0.5 using the analytic - canonical normal. The internal boundary has normals pointing in +y or - -y direction, so integrating n_y should give +1 or -1 (length 1 - boundary).""" + """Integrate n_y along internal boundary at y=0.5 with plain mesh.Gamma. + The declared internal normal is +y, so the integral is exactly +1 + (length-1 boundary).""" mesh_internal, _, _ = _get_internal_mesh() - normal = mesh_internal.canonical_normal("Internal") - n_y = normal[1] + n_y = mesh_internal.Gamma[1] bd_int = uw.maths.BdIntegral(mesh_internal, fn=n_y, boundary="Internal") value = bd_int.evaluate() - # Normal orientation is consistent but direction depends on mesh; - # absolute value should be 1.0 - assert abs(abs(value) - 1.0) < 0.01, f"Expected |n_y integral| = 1.0, got {value}" + assert abs(value - 1.0) < 1e-6, f"Expected +1.0, got {value}" def test_bd_integral_internal_normal_nx(): @@ -195,33 +190,44 @@ def test_bd_integral_internal_normal_nx(): bd_int = uw.maths.BdIntegral(mesh_internal, fn=n_x, boundary="Internal") value = bd_int.evaluate() - assert abs(value) < 0.01, f"Expected ~0, got {value}" + assert abs(value) < 1e-6, f"Expected ~0, got {value}" def test_bd_integral_internal_normal_weighted(): - """Integrate x * n_y along internal boundary at y=0.5 using the analytic - canonical normal. int_0^1 x * n_y dx = n_y * 0.5. Since |n_y| = 1, - result should be ~0.5.""" + """Integrate x * n_y along internal boundary at y=0.5 with plain + mesh.Gamma: int_0^1 x * (+1) dx = +0.5.""" mesh_internal, x_i, _ = _get_internal_mesh() + n_y = mesh_internal.Gamma[1] + + bd_int = uw.maths.BdIntegral(mesh_internal, fn=x_i * n_y, boundary="Internal") + value = bd_int.evaluate() + + assert abs(value - 0.5) < 1e-6, f"Expected +0.5, got {value}" + + +def test_bd_integral_internal_canonical_normal_accessor(): + """The canonical_normal accessor remains available and agrees with the + normal that mesh.Gamma resolves to on the internal boundary.""" + + mesh_internal, _, _ = _get_internal_mesh() normal = mesh_internal.canonical_normal("Internal") n_y = normal[1] - bd_int = uw.maths.BdIntegral(mesh_internal, fn=x_i * n_y, boundary="Internal") + bd_int = uw.maths.BdIntegral(mesh_internal, fn=n_y, boundary="Internal") value = bd_int.evaluate() - assert abs(abs(value) - 0.5) < 0.01, f"Expected |value| = 0.5, got {value}" + assert abs(value - 1.0) < 1e-6, f"Expected +1.0, got {value}" -def test_bd_integral_internal_canonical_normal_off_grid_zint(): +def test_bd_integral_internal_gamma_off_grid_zint(): """Regression for the failing partition-through-boundary case from #327. With ``zintCoord=0.55`` (off-grid), the mpirun -n 2 partition seam runs - through the internal boundary and one seam facet's ``petsc_n[]`` flips - sign: ``∫ Gamma[1] dS`` returns 0.9375 = 1 − 2/32 instead of 1.0. The - analytic ``canonical_normal("Internal")`` returned by the mesh factory - does not touch ``petsc_n[]`` and gives the exact value regardless of - partition.""" + through the internal boundary and one seam facet's raw ``petsc_n[]`` + flips sign: the unresolved integral returned 0.9375 = 1 − 2/32 instead + of 1.0. With plain ``mesh.Gamma`` now resolved to the declared analytic + normal, the value is exact regardless of partition.""" mesh_off = BoxInternalBoundary( minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0/32.0, zintCoord=0.55, simplex=True, @@ -229,10 +235,10 @@ def test_bd_integral_internal_canonical_normal_off_grid_zint(): # Need at least one variable so BdIntegral has a section to integrate against uw.discretisation.MeshVariable("T_off", mesh_off, 1, degree=2) - n_y = mesh_off.canonical_normal("Internal")[1] + n_y = mesh_off.Gamma[1] val = uw.maths.BdIntegral(mesh_off, fn=n_y, boundary="Internal").evaluate() - assert abs(abs(val) - 1.0) < 1e-6, ( - f"canonical_normal internal integral should be exactly ±1, got {val}") + assert abs(val - 1.0) < 1e-6, ( + f"mesh.Gamma internal integral should be exactly +1, got {val}") def test_bd_integral_internal_does_not_affect_external(): From 2f98dec760190796161fd52f3d4aab8decf690a3 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 14 Jul 2026 11:16:19 +1000 Subject: [PATCH 4/4] fix(mesh): invalidate declared analytic boundary normals on deformation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The factory-declared analytic normal (canonical_normal, consumed when mesh.Gamma resolves on an internal boundary) describes the ORIGINAL factory geometry. After mesh.deform() / adapt() the internal surface may have warped, and silently integrating the stale declaration would be wrong. - nuke_coords_and_rebuild bumps a _geometry_version counter (0 at construction; every later call is a genuine coordinate change — deform, adapt, re-extract and direct coordinate writes all funnel through it). - boundary_normals is now a property whose setter stamps the declaration with the current geometry version. - canonical_normal raises RuntimeError when the declaration predates the current geometry, telling the user to re-declare (mesh.boundary_normals = mesh.boundary_normals when the expression is still valid, e.g. radial normals after a radius-preserving remesh) or to use an explicit normal expression. External boundaries are unaffected: Gamma -> petsc_n[] is exact on the deformed geometry and never consults the declaration. Test: deform a BoxInternalBoundary with a displacement vanishing on the internal plane; Gamma resolution raises; re-declaration restores the exact +1 integral. Serial-only until #360 (deform crashes at np>1). Verified: 0502 25 passed serial, 24 passed + 1 skip at np=2 and np=4; natural-BC/Nitsche/constrained regression trio 23 passed. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 51 ++++++++++++++++++- tests/test_0502_boundary_integrals.py | 36 +++++++++++++ 2 files changed, 86 insertions(+), 1 deletion(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 3ee2d31d..c6ee02b8 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2264,6 +2264,12 @@ def nuke_coords_and_rebuild( # # + # Geometry generation counter. First call (construction) -> 0; every + # later call means the coordinates changed (deform / adapt / + # re-extract), which invalidates factory-declared analytic boundary + # normals (see the boundary_normals setter and canonical_normal). + self._geometry_version = getattr(self, "_geometry_version", -1) + 1 + self.dm.clearDS() self.dm.createDS() @@ -2723,6 +2729,25 @@ def Gamma_P1(self): self._update_projected_normals() return self._projected_normals.sym + @property + def boundary_normals(self): + """Declared analytic boundary normals (Enum or mapping), or ``None``. + + Assigning stamps the declaration against the mesh's current + geometry: any later coordinate change (``deform``, ``adapt``, + direct coordinate writes) marks it stale, and + :meth:`canonical_normal` then refuses to serve it. Re-assign after + a deformation to re-declare normals that are valid for the new + geometry (e.g. a radial normal after a radius-preserving remesh). + """ + return getattr(self, "_boundary_normals", None) + + @boundary_normals.setter + def boundary_normals(self, value): + self._boundary_normals = value + self._boundary_normals_geometry_version = getattr( + self, "_geometry_version", 0) + def canonical_normal(self, boundary_name): r"""Analytic outward-pointing normal for a boundary, or ``None`` if no analytic normal was declared for that boundary. @@ -2760,6 +2785,15 @@ def canonical_normal(self, boundary_name): normal, or ``None`` if this mesh factory did not declare an analytic normal for ``boundary_name``. + Raises + ------ + RuntimeError + If the mesh coordinates have changed since the normals were + declared (``deform`` / ``adapt``): the declaration describes + the original geometry and may no longer match the deformed + surface. Re-assign :attr:`boundary_normals` to re-declare + normals valid for the current geometry. + See Also -------- Gamma : the raw per-quadrature normal — use for external @@ -2768,9 +2802,24 @@ def canonical_normal(self, boundary_name): Gamma_P1 : projected P1 normals, useful for curved external boundaries. """ - bn = getattr(self, "boundary_normals", None) + bn = self.boundary_normals if bn is None: return None + declared_at = getattr(self, "_boundary_normals_geometry_version", 0) + if declared_at != getattr(self, "_geometry_version", 0): + raise RuntimeError( + f"The declared analytic boundary normals for this mesh " + f"describe its original (factory) geometry, but the mesh " + f"coordinates have changed since (deform / adapt). The " + f"declaration for '{boundary_name}' may no longer match the " + f"deformed surface. If the normals are still valid for the " + f"new geometry (e.g. a radial normal after a " + f"radius-preserving remesh), re-declare them:\n" + f" mesh.boundary_normals = mesh.boundary_normals\n" + f"or assign a new Enum / dict of normal expressions. " + f"Otherwise use an explicit normal expression in place of " + f"mesh.Gamma / canonical_normal on this boundary." + ) try: member = bn[boundary_name] except (KeyError, AttributeError): diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py index b32cca7a..9e77f9e1 100644 --- a/tests/test_0502_boundary_integrals.py +++ b/tests/test_0502_boundary_integrals.py @@ -241,6 +241,42 @@ def test_bd_integral_internal_gamma_off_grid_zint(): f"mesh.Gamma internal integral should be exactly +1, got {val}") +@pytest.mark.skipif( + uw.mpi.size > 1, + reason="mesh.deform crashes at np>1 (issue #360, kd-tree index rebuild)", +) +def test_bd_integral_internal_gamma_stale_after_deform(): + """Deformation invalidates the factory-declared analytic normal. + + The declaration describes the original geometry; after mesh.deform() + resolving mesh.Gamma on the internal boundary must fail loudly rather + than integrate a stale normal. Re-assigning mesh.boundary_normals + re-declares it for the new geometry. The deformation used here + vanishes on y=0.5 (sin(2*pi*y) = 0), so the internal boundary is + unmoved and the re-declared +y normal gives exactly +1 again.""" + + mesh_d = BoxInternalBoundary( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=1.0/16.0, zintCoord=0.5, simplex=True, + ) + uw.discretisation.MeshVariable("T_deform", mesh_d, 1, degree=2) + n_y = mesh_d.Gamma[1] + + coords = np.array(mesh_d.X.coords) + new_coords = coords.copy() + new_coords[:, 1] += ( + 0.01 * np.sin(np.pi * coords[:, 0]) * np.sin(2.0 * np.pi * coords[:, 1]) + ) + mesh_d.deform(new_coords) + + with pytest.raises(RuntimeError, match="coordinates have changed"): + uw.maths.BdIntegral(mesh_d, fn=n_y, boundary="Internal").evaluate() + + mesh_d.boundary_normals = mesh_d.boundary_normals + val = uw.maths.BdIntegral(mesh_d, fn=n_y, boundary="Internal").evaluate() + assert abs(val - 1.0) < 1e-6, f"Expected +1.0 after re-declaration, got {val}" + + def test_bd_integral_internal_does_not_affect_external(): """External boundaries should still work on the internal-boundary mesh."""