Unify boundary-normal interface behind mesh.Gamma (#327) + np=4 label deadlock fix

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

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':

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
 

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()
 
@@ -2534,17 +2540,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 +2563,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 +2573,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,19 +2713,204 @@ 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.
         """
         if not hasattr(self, '_projected_normals') or self._projected_normals is None:
             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.
+
+        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``.
+
+        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
+            boundaries; may be partition-dependent on internal seams.
+        Gamma_N : normalised :attr:`Gamma`.
+        Gamma_P1 : projected P1 normals, useful for curved external
+            boundaries.
+        """
+        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):
+            return None
+        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
@@ -3453,11 +3655,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
         -------
@@ -3469,15 +3671,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