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
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 8859464a..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()
 
@@ -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,12 +2713,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.
         """
@@ -2709,6 +2729,188 @@ 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.
+
+        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
 
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 5f7d323f..9e77f9e1 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,31 +159,24 @@ 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)",
-)
+# `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.
-    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()
-    Gamma = mesh_internal.Gamma
-    n_y = Gamma[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():
@@ -197,25 +190,91 @@ 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}"
 
 
-@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 with plain
+    mesh.Gamma: int_0^1 x * (+1) dx = +0.5."""
 
     mesh_internal, x_i, _ = _get_internal_mesh()
-    Gamma = mesh_internal.Gamma
-    n_y = Gamma[1]
+    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(abs(value) - 0.5) < 0.01, f"Expected |value| = 0.5, got {value}"
+    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=n_y, boundary="Internal")
+    value = bd_int.evaluate()
+
+    assert abs(value - 1.0) < 1e-6, f"Expected +1.0, got {value}"
+
+
+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 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,
+    )
+    # 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.Gamma[1]
+    val = uw.maths.BdIntegral(mesh_off, fn=n_y, boundary="Internal").evaluate()
+    assert abs(val - 1.0) < 1e-6, (
+        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():