fix(mesh): add mesh.canonical_normal accessor for partition-safe internal-boundary normals (#327)

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

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,24 +198,43 @@ 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()
 
     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."""