From 992d538e78e6aa9a302c07623916dad0cb4e21fd Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Tue, 3 Mar 2026 17:20:32 +1100
Subject: [PATCH 1/8] =?UTF-8?q?Add=20boundary=20integral=20support=20(BdIn?=
 =?UTF-8?q?tegral)=20=E2=80=94=20closes=20#47?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Wrap PETSc's DMPlexComputeBdIntegral to provide standalone boundary/surface
integral capability. In 2D these are line integrals over boundary edges;
in 3D they are surface integrals over boundary faces (both codimension-1).

The integrand may reference the outward unit normal via mesh.Gamma / mesh.Gamma_N.
Internal boundaries (e.g. AnnulusInternalBoundary, BoxInternalBoundary) work
out of the box through the same DMLabel infrastructure.

New API:
    bd_int = uw.maths.BdIntegral(mesh, fn=1.0, boundary="Top")
    length = bd_int.evaluate()

Tests cover: constant/coordinate/sympy/mesh-variable integrands, normal-vector
integrands, perimeter checks, BoxInternalBoundary, and AnnulusInternalBoundary
circumference.

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_compat.h   |  39 +++++
 src/underworld3/cython/petsc_extras.pxi |   3 +-
 src/underworld3/cython/petsc_maths.pyx  | 179 ++++++++++++++++++++++
 src/underworld3/maths/__init__.py       |   3 +
 tests/test_0502_boundary_integrals.py   | 194 ++++++++++++++++++++++++
 5 files changed, 417 insertions(+), 1 deletion(-)
 create mode 100644 tests/test_0502_boundary_integrals.py

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 1b7cb31b..3de05619 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -138,3 +138,42 @@ PetscErrorCode UW_DMPlexSetSNESLocalFEM(DM dm, PetscBool flag, void *ctx)
     return DMPlexSetSNESLocalFEM(dm, flag, NULL);
 #endif
 }
+
+// Simplified wrapper for DMPlexComputeBdIntegral.
+// Takes a single boundary pointwise function (for field 0) instead of an Nf-element array.
+// Internally allocates the full array with NULL for all other fields.
+PetscErrorCode UW_DMPlexComputeBdIntegral(DM dm, Vec X,
+                                          DMLabel label, PetscInt numVals, const PetscInt vals[],
+                                          void (*func)(UW_SIG_F0),
+                                          PetscScalar *result,
+                                          void *ctx)
+{
+    PetscDS       ds;
+    PetscSection  section;
+    PetscInt      Nf;
+
+    PetscFunctionBeginUser;
+
+    PetscCall(DMGetLocalSection(dm, &section));
+    PetscCall(PetscSectionGetNumFields(section, &Nf));
+
+    // Allocate function pointer array (all NULL except field 0)
+    void (**funcs)(UW_SIG_F0);
+    PetscCall(PetscCalloc1(Nf, &funcs));
+    funcs[0] = func;
+
+    // Allocate output array (one value per field)
+    PetscScalar *integral;
+    PetscCall(PetscCalloc1(Nf, &integral));
+
+    // Compute the boundary integral
+    PetscCall(DMPlexComputeBdIntegral(dm, X, label, numVals, vals, funcs, integral, ctx));
+
+    // Return field 0 result
+    *result = integral[0];
+
+    PetscCall(PetscFree(funcs));
+    PetscCall(PetscFree(integral));
+
+    PetscFunctionReturn(PETSC_SUCCESS);
+}
diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
index 62bf7407..c4cad934 100644
--- a/src/underworld3/cython/petsc_extras.pxi
+++ b/src/underworld3/cython/petsc_extras.pxi
@@ -41,8 +41,9 @@ cdef extern from "petsc_compat.h":
     PetscErrorCode UW_PetscDSSetBdJacobianPreconditioner(PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, void*, void*, void*, void*)
     PetscErrorCode UW_PetscDSSetBdTerms   (PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, void*, void*, void*, void*, void*, void* )
     PetscErrorCode UW_PetscDSViewWF(PetscDS)     
-    PetscErrorCode UW_PetscDSViewBdWF(PetscDS, PetscInt)     
+    PetscErrorCode UW_PetscDSViewBdWF(PetscDS, PetscInt)
     PetscErrorCode UW_DMPlexSetSNESLocalFEM( PetscDM, PetscBool, void *)
+    PetscErrorCode UW_DMPlexComputeBdIntegral( PetscDM, PetscVec, PetscDMLabel, PetscInt, const PetscInt*, void*, PetscScalar*, void*)
 
 cdef extern from "petsc.h" nogil:
     PetscErrorCode DMPlexSNESComputeBoundaryFEM( PetscDM, void *, void *)
diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx
index eb99be7e..529f68ca 100644
--- a/src/underworld3/cython/petsc_maths.pyx
+++ b/src/underworld3/cython/petsc_maths.pyx
@@ -300,3 +300,182 @@ class CellWiseIntegral:
         rvec.destroy()
 
         return results
+
+
+class BdIntegral:
+    r"""
+    The ``BdIntegral`` class constructs the boundary (surface/line) integral
+
+    .. math:: F  =  \int_{\partial \Omega} \, f(\mathbf{x}, \mathbf{n}) \, \mathrm{d} S
+
+    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).
+
+    Parameters
+    ----------
+    mesh : underworld3.discretisation.Mesh
+        The mesh whose boundary is integrated over.
+    fn : float, int, or sympy.Basic
+        Scalar function to be integrated.
+    boundary : str
+        Name of the boundary label (e.g. ``"Top"``, ``"Bottom"``, ``"Internal"``).
+        Must be present in ``mesh.boundaries``.
+
+    Example
+    -------
+    Calculate length of the top boundary of a unit box:
+
+    >>> import underworld3 as uw
+    >>> import numpy as np
+    >>> mesh = uw.discretisation.Box()
+    >>> bdIntegral = uw.maths.BdIntegral(mesh=mesh, fn=1., boundary="Top")
+    >>> np.allclose( 1., bdIntegral.evaluate(), rtol=1e-8)
+    True
+
+    See Also
+    --------
+    Integral : For volume integrals over the entire mesh domain.
+    """
+
+    @timing.routine_timer_decorator
+    def __init__( self,
+                  mesh:     underworld3.discretisation.Mesh,
+                  fn:       Union[float, int, sympy.Basic],
+                  boundary: str ):
+
+        self.mesh = mesh
+        self.fn = sympy.sympify(fn)
+        self.boundary = boundary
+
+        # Validate boundary name
+        if boundary not in [b.name for b in mesh.boundaries]:
+            available = [b.name for b in mesh.boundaries]
+            raise ValueError(
+                f"Boundary '{boundary}' not found on mesh. "
+                f"Available boundaries: {available}"
+            )
+
+        super().__init__()
+
+    @timing.routine_timer_decorator
+    def evaluate(self, verbose=False):
+        """
+        Evaluate the boundary integral and return the result.
+
+        Returns
+        -------
+        float or UWQuantity
+            The boundary integral value. If the integrand has units AND the
+            mesh coordinates have units, returns a UWQuantity with proper
+            dimensional analysis (integrand_units * coord_units^(dim-1)).
+            Otherwise returns a plain float.
+        """
+        if len(self.mesh.vars) == 0:
+            raise RuntimeError(
+                "The mesh requires at least a single variable for integration "
+                "to function correctly.\nThis is a PETSc limitation."
+            )
+
+        if isinstance(self.fn, sympy.vector.Vector):
+            raise RuntimeError("BdIntegral evaluation for Vector integrands not supported.")
+        elif isinstance(self.fn, sympy.vector.Dyadic):
+            raise RuntimeError("BdIntegral evaluation for Dyadic integrands not supported.")
+
+        mesh = self.mesh
+
+        # Compile integrand using the boundary residual slot (includes petsc_n[] in signature)
+        compiled_extns, dictionaries = getext(
+            self.mesh, [], [], [], [self.fn,], [], self.mesh.vars.values(), verbose=verbose
+        )
+        cdef PtrContainer ext = compiled_extns
+
+        # Prepare the solution vector
+        self.mesh.update_lvec()
+        a_global = mesh.dm.getGlobalVec()
+        mesh.dm.localToGlobal(self.mesh.lvec, a_global)
+
+        cdef Vec cgvec = a_global
+        cdef DM dm_c = mesh.dm
+
+        # Get the DMLabel and label value for the named boundary
+        boundary_enum = mesh.boundaries[self.boundary]
+        cdef PetscInt label_val = boundary_enum.value
+        cdef PetscInt num_vals = 1
+
+        c_label = mesh.dm.getLabel(self.boundary)
+        cdef DMLabel dmlabel = c_label
+
+        # Output value
+        cdef PetscScalar result = 0.0
+
+        # Call the boundary integral
+        ierr = UW_DMPlexComputeBdIntegral(
+            dm_c.dm, cgvec.vec,
+            dmlabel.dmlabel, num_vals, &label_val,
+            ext.fns_bd_residual[0],
+            &result, NULL
+        )
+        CHKERRQ(ierr)
+
+        mesh.dm.restoreGlobalVec(a_global)
+
+        cdef double vald = <double> result
+
+        # Unit propagation: result_units = integrand_units * coord_units^(dim-1)
+        # Boundary integral is over a codimension-1 manifold.
+        try:
+            integrand_units = underworld3.get_units(self.fn)
+            coord_units = underworld3.get_units(self.mesh.X[0])
+
+            from underworld3.scaling import units as ureg
+
+            def has_meaningful_units(u):
+                if u is None:
+                    return False
+                try:
+                    if u == ureg.dimensionless:
+                        return False
+                except:
+                    pass
+                return True
+
+            integrand_has_units = has_meaningful_units(integrand_units)
+            coord_has_units = has_meaningful_units(coord_units)
+
+            if integrand_has_units or coord_has_units:
+
+                if coord_has_units:
+                    surface_units = coord_units ** (self.mesh.dim - 1)
+                else:
+                    surface_units = None
+
+                if integrand_has_units and surface_units is not None:
+                    result_units = integrand_units * surface_units
+                elif integrand_has_units:
+                    result_units = integrand_units
+                elif surface_units is not None:
+                    result_units = surface_units
+                else:
+                    return vald
+
+                physical_value = vald
+
+                if underworld3.is_nondimensional_scaling_active() and integrand_has_units:
+                    try:
+                        orchestration_model = underworld3.get_default_model()
+                        if orchestration_model.has_units():
+                            integrand_dimensionality = (1 * integrand_units).dimensionality
+                            integrand_scale = orchestration_model.get_scale_for_dimensionality(integrand_dimensionality)
+                            integrand_scale_converted = integrand_scale.to(integrand_units)
+                            physical_value = vald * float(integrand_scale_converted.magnitude)
+                    except Exception:
+                        pass
+
+                return underworld3.quantity(physical_value, result_units)
+        except Exception:
+            pass
+
+        return vald
diff --git a/src/underworld3/maths/__init__.py b/src/underworld3/maths/__init__.py
index ffa774e3..1cecae6d 100644
--- a/src/underworld3/maths/__init__.py
+++ b/src/underworld3/maths/__init__.py
@@ -24,6 +24,8 @@
 
 **CellWiseIntegral** -- Cell-by-cell integration.
 
+**BdIntegral** -- Boundary (surface/line) integration.
+
 See Also
 --------
 underworld3.coordinates : Coordinate system definitions.
@@ -52,3 +54,4 @@
 # These could be wrapped so that they can be documented along with the math module
 from underworld3.cython.petsc_maths import Integral
 from underworld3.cython.petsc_maths import CellWiseIntegral
+from underworld3.cython.petsc_maths import BdIntegral
diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py
new file mode 100644
index 00000000..3477194a
--- /dev/null
+++ b/tests/test_0502_boundary_integrals.py
@@ -0,0 +1,194 @@
+import underworld3 as uw
+import numpy as np
+import sympy
+import pytest
+
+# All tests in this module are quick core tests
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
+
+from underworld3.meshing import UnstructuredSimplexBox
+
+## Set up the mesh for tests
+mesh = UnstructuredSimplexBox(
+    minCoords=(0.0, 0.0),
+    maxCoords=(1.0, 1.0),
+    cellSize=1.0 / 32.0,
+)
+
+x, y = mesh.X
+
+# Need at least one mesh variable for PETSc integration (same as volume Integral)
+s_soln = uw.discretisation.MeshVariable("T_bd", mesh, 1, degree=2)
+
+
+def test_bd_integral_constant_top():
+    """Integrating 1.0 over the Top boundary of a unit box should give length 1.0."""
+
+    bd_int = uw.maths.BdIntegral(mesh, fn=1.0, boundary="Top")
+    value = bd_int.evaluate()
+
+    assert abs(value - 1.0) < 0.001, f"Expected 1.0, got {value}"
+
+
+def test_bd_integral_constant_bottom():
+    """Integrating 1.0 over the Bottom boundary of a unit box should give length 1.0."""
+
+    bd_int = uw.maths.BdIntegral(mesh, fn=1.0, boundary="Bottom")
+    value = bd_int.evaluate()
+
+    assert abs(value - 1.0) < 0.001, f"Expected 1.0, got {value}"
+
+
+def test_bd_integral_perimeter():
+    """Sum of integrals over all four boundaries should give perimeter = 4.0."""
+
+    total = 0.0
+    for bnd in ["Top", "Bottom", "Left", "Right"]:
+        bd_int = uw.maths.BdIntegral(mesh, fn=1.0, boundary=bnd)
+        total += bd_int.evaluate()
+
+    assert abs(total - 4.0) < 0.01, f"Expected perimeter 4.0, got {total}"
+
+
+def test_bd_integral_coordinate_fn():
+    """Integrate x along Top boundary (y=1): int_0^1 x dx = 0.5."""
+
+    bd_int = uw.maths.BdIntegral(mesh, fn=x, boundary="Top")
+    value = bd_int.evaluate()
+
+    assert abs(value - 0.5) < 0.01, f"Expected 0.5, got {value}"
+
+
+def test_bd_integral_coordinate_fn_right():
+    """Integrate y along Right boundary (x=1): int_0^1 y dy = 0.5."""
+
+    bd_int = uw.maths.BdIntegral(mesh, fn=y, boundary="Right")
+    value = bd_int.evaluate()
+
+    assert abs(value - 0.5) < 0.01, f"Expected 0.5, got {value}"
+
+
+def test_bd_integral_sympy_fn():
+    """Integrate cos(pi*x) along Top boundary: int_0^1 cos(pi*x) dx = 0."""
+
+    bd_int = uw.maths.BdIntegral(mesh, fn=sympy.cos(x * sympy.pi), boundary="Top")
+    value = bd_int.evaluate()
+
+    assert abs(value) < 0.01, f"Expected ~0, got {value}"
+
+
+def test_bd_integral_meshvar():
+    """Integrate a mesh variable (sin(pi*x)) along Top boundary."""
+
+    s_soln.data[:, 0] = np.sin(np.pi * s_soln.coords[:, 0])
+
+    bd_int = uw.maths.BdIntegral(mesh, fn=s_soln.sym[0], boundary="Top")
+    value = bd_int.evaluate()
+
+    # int_0^1 sin(pi*x) dx = 2/pi ≈ 0.6366
+    expected = 2.0 / np.pi
+    assert abs(value - expected) < 0.01, f"Expected {expected}, got {value}"
+
+
+def test_bd_integral_normal_vector():
+    """Integrand using surface normal: integrate n_y along Top boundary.
+    On Top boundary (y=1), the outward normal is (0, 1), so n_y = 1.
+    Integral should be 1.0."""
+
+    Gamma = mesh.Gamma  # Surface normal as row matrix
+    n_y = Gamma[1]
+
+    bd_int = uw.maths.BdIntegral(mesh, fn=n_y, boundary="Top")
+    value = bd_int.evaluate()
+
+    assert abs(value - 1.0) < 0.01, f"Expected 1.0, got {value}"
+
+
+def test_bd_integral_invalid_boundary():
+    """Should raise ValueError for non-existent boundary name."""
+
+    with pytest.raises(ValueError, match="not found"):
+        uw.maths.BdIntegral(mesh, fn=1.0, boundary="Nonexistent")
+
+
+# --- Internal boundary tests ---
+
+from underworld3.meshing import BoxInternalBoundary
+
+mesh_internal = BoxInternalBoundary(
+    minCoords=(0.0, 0.0),
+    maxCoords=(1.0, 1.0),
+    cellSize=1.0 / 32.0,
+    zintCoord=0.5,
+    simplex=True,
+)
+
+x_i, y_i = mesh_internal.X
+
+# Need a mesh variable for PETSc
+_v_internal = uw.discretisation.MeshVariable("T_int", mesh_internal, 1, degree=2)
+
+
+def test_bd_integral_internal_boundary_length():
+    """Internal boundary at y=0.5 across a unit box should have length 1.0."""
+
+    bd_int = uw.maths.BdIntegral(mesh_internal, fn=1.0, boundary="Internal")
+    value = bd_int.evaluate()
+
+    assert abs(value - 1.0) < 0.001, f"Expected 1.0, got {value}"
+
+
+def test_bd_integral_internal_coordinate_fn():
+    """Integrate x along internal boundary at y=0.5: int_0^1 x dx = 0.5."""
+
+    bd_int = uw.maths.BdIntegral(mesh_internal, fn=x_i, boundary="Internal")
+    value = bd_int.evaluate()
+
+    assert abs(value - 0.5) < 0.01, f"Expected 0.5, got {value}"
+
+
+def test_bd_integral_internal_does_not_affect_external():
+    """External boundaries should still work on the internal-boundary mesh."""
+
+    total = 0.0
+    for bnd in ["Top", "Bottom", "Left", "Right"]:
+        bd_int = uw.maths.BdIntegral(mesh_internal, fn=1.0, boundary=bnd)
+        total += bd_int.evaluate()
+
+    assert abs(total - 4.0) < 0.01, f"Expected perimeter 4.0, got {total}"
+
+
+# --- Annulus internal boundary tests ---
+
+from underworld3.meshing import AnnulusInternalBoundary
+
+_R_INTERNAL = 1.0
+
+mesh_annulus = AnnulusInternalBoundary(
+    radiusOuter=1.5,
+    radiusInternal=_R_INTERNAL,
+    radiusInner=0.5,
+    cellSize=0.1,
+)
+
+_v_annulus = uw.discretisation.MeshVariable("T_ann", mesh_annulus, 1, degree=2)
+
+
+def test_bd_integral_annulus_internal_circumference():
+    """Internal boundary at radius 1.0: circumference = 2*pi."""
+
+    bd_int = uw.maths.BdIntegral(mesh_annulus, fn=1.0, boundary="Internal")
+    value = bd_int.evaluate()
+    expected = 2 * np.pi * _R_INTERNAL
+
+    assert abs(value - expected) < 0.01, f"Expected {expected:.4f}, got {value}"
+
+
+def test_bd_integral_annulus_outer_circumference():
+    """Outer boundary at radius 1.5: circumference = 3*pi."""
+
+    bd_int = uw.maths.BdIntegral(mesh_annulus, fn=1.0, boundary="Upper")
+    value = bd_int.evaluate()
+    expected = 2 * np.pi * 1.5
+
+    assert abs(value - expected) < 0.01, f"Expected {expected:.4f}, got {value}"

From b67b07835a176d5e667ce0eafbe6b2e3c591c05e Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Tue, 3 Mar 2026 17:33:52 +1100
Subject: [PATCH 2/8] Clarify AI attribution convention for PRs in CLAUDE.md

Use the same team attribution wording for both commits and PR descriptions.
No emoji in PR bodies.

Underworld development team with AI support from Claude Code
---
 CLAUDE.md | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/CLAUDE.md b/CLAUDE.md
index 7a3a5356..3f3d68ba 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -198,16 +198,20 @@ worktree changes, the worktree edits will not be active. Always:
 2. `./uw build`
 3. Run your code / tests from there
 
-### AI-Assisted Commit Attribution
-When committing code developed with AI assistance, end the commit message with:
+### AI-Assisted Attribution (Commits and PRs)
+When committing code or creating pull requests with AI assistance, end the
+message/body with:
 
 ```
-Underworld development team with AI support from Claude Code
+Underworld development team with AI support from [Claude Code](https://claude.com/claude-code)
 ```
 
+(In commit messages, use the plain-text form without the markdown link.)
+
 **Do NOT use**:
 - `Co-Authored-By:` with a noreply email (useless for soliciting responses)
 - Generic AI attribution without team context
+- Emoji in PR descriptions
 
 ---
 

From f2018a710527693486043535f923945551e30148 Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Tue, 3 Mar 2026 20:19:08 +1100
Subject: [PATCH 3/8] Add MPI Allreduce to boundary integral wrapper
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

DMPlexComputeBdIntegral (unlike DMPlexComputeIntegralFEM) does not
perform an MPI reduction — it returns only the local process contribution.
Add MPIU_Allreduce in UW_DMPlexComputeBdIntegral to sum across ranks.

Found by testing with mpirun -np 5..8 on AnnulusInternalBoundary.

Note: a pre-existing PETSc internal-facet orientation issue can cause
small errors (~1.5%) on internal boundaries at certain partition counts
(e.g. np=6 with the default annulus mesh). External boundaries are
unaffected. This is upstream of our wrapper.

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_compat.h | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 3de05619..285c83f4 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -166,11 +166,17 @@ PetscErrorCode UW_DMPlexComputeBdIntegral(DM dm, Vec X,
     PetscScalar *integral;
     PetscCall(PetscCalloc1(Nf, &integral));
 
-    // Compute the boundary integral
+    // Compute the boundary integral (returns local contribution only —
+    // DMPlexComputeBdIntegral is missing the MPI reduction that
+    // DMPlexComputeIntegralFEM has at plexfem.c:2633)
     PetscCall(DMPlexComputeBdIntegral(dm, X, label, numVals, vals, funcs, integral, ctx));
 
-    // Return field 0 result
-    *result = integral[0];
+    // Sum across all MPI ranks (mirrors DMPlexComputeIntegralFEM behaviour)
+    PetscScalar global_val;
+    PetscCallMPI(MPIU_Allreduce(&integral[0], &global_val, 1, MPIU_SCALAR, MPIU_SUM,
+                                PetscObjectComm((PetscObject)dm)));
+
+    *result = global_val;
 
     PetscCall(PetscFree(funcs));
     PetscCall(PetscFree(integral));

From bfc0b543fcc683ab36ba6ff68f016570276002ca Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Tue, 3 Mar 2026 21:01:04 +1100
Subject: [PATCH 4/8] Fix ghost facet double-counting in boundary integrals
 under MPI
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

DMPlexComputeBdIntegral iterates over ALL label stratum points including
ghost facets, so shared internal boundary facets get integrated on both
the owning rank and the ghost rank. External boundaries are unaffected
because external facets have only one supporting cell (no ghosting).

The fix creates a temporary DMLabel containing only owned (non-ghost)
boundary points by checking against the PetscSF leaf set, then passes
this filtered label to DMPlexComputeBdIntegral.

Also restructures tests to use lazy initialization for BoxInternalBoundary
(which has a pre-existing MPI bug) so it doesn't block other tests under
mpirun.

Verified at np=1,5,6,7,8 — all produce identical results.

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_compat.h | 64 ++++++++++++++++++++++-----
 tests/test_0502_boundary_integrals.py | 34 +++++++++-----
 2 files changed, 76 insertions(+), 22 deletions(-)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 285c83f4..649cc512 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -141,43 +141,85 @@ PetscErrorCode UW_DMPlexSetSNESLocalFEM(DM dm, PetscBool flag, void *ctx)
 
 // Simplified wrapper for DMPlexComputeBdIntegral.
 // Takes a single boundary pointwise function (for field 0) instead of an Nf-element array.
-// Internally allocates the full array with NULL for all other fields.
+//
+// Fixes two issues in DMPlexComputeBdIntegral:
+//   1. Missing MPI Allreduce (plexfem.c returns local contribution only,
+//      unlike DMPlexComputeIntegralFEM which reduces at plexfem.c:2633).
+//   2. Ghost facet double-counting: DMPlexComputeBdIntegral iterates over ALL
+//      label stratum points including ghost facets, so shared boundary facets
+//      get integrated on both the owning rank and the ghost rank. We create a
+//      temporary label containing only owned (non-ghost) boundary points.
 PetscErrorCode UW_DMPlexComputeBdIntegral(DM dm, Vec X,
                                           DMLabel label, PetscInt numVals, const PetscInt vals[],
                                           void (*func)(UW_SIG_F0),
                                           PetscScalar *result,
                                           void *ctx)
 {
-    PetscDS       ds;
     PetscSection  section;
-    PetscInt      Nf;
+    PetscInt      Nf, pStart, pEnd;
 
     PetscFunctionBeginUser;
 
     PetscCall(DMGetLocalSection(dm, &section));
     PetscCall(PetscSectionGetNumFields(section, &Nf));
 
-    // Allocate function pointer array (all NULL except field 0)
+    // --- Build a ghost-point bitset from the point SF ---
+    PetscSF   sf;
+    PetscInt  nleaves;
+    const PetscInt *ilocal;
+    PetscCall(DMGetPointSF(dm, &sf));
+    PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &ilocal, NULL));
+    PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
+
+    PetscBT ghostBT;
+    PetscCall(PetscBTCreate(pEnd - pStart, &ghostBT));
+    if (ilocal) {
+        for (PetscInt i = 0; i < nleaves; i++) {
+            PetscCall(PetscBTSet(ghostBT, ilocal[i] - pStart));
+        }
+    }
+
+    // --- Create a temporary label with only owned boundary points ---
+    DMLabel ownedLabel;
+    PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)dm), "uw_owned_bd", &ownedLabel));
+
+    for (PetscInt v = 0; v < numVals; v++) {
+        IS origIS;
+        PetscCall(DMLabelGetStratumIS(label, vals[v], &origIS));
+        if (!origIS) continue;
+
+        PetscInt        n;
+        const PetscInt *indices;
+        PetscCall(ISGetLocalSize(origIS, &n));
+        PetscCall(ISGetIndices(origIS, &indices));
+
+        for (PetscInt i = 0; i < n; i++) {
+            if (!PetscBTLookup(ghostBT, indices[i] - pStart)) {
+                PetscCall(DMLabelSetValue(ownedLabel, indices[i], vals[v]));
+            }
+        }
+        PetscCall(ISRestoreIndices(origIS, &indices));
+        PetscCall(ISDestroy(&origIS));
+    }
+    PetscCall(PetscBTDestroy(&ghostBT));
+
+    // --- Compute boundary integral over owned points only ---
     void (**funcs)(UW_SIG_F0);
     PetscCall(PetscCalloc1(Nf, &funcs));
     funcs[0] = func;
 
-    // Allocate output array (one value per field)
     PetscScalar *integral;
     PetscCall(PetscCalloc1(Nf, &integral));
 
-    // Compute the boundary integral (returns local contribution only —
-    // DMPlexComputeBdIntegral is missing the MPI reduction that
-    // DMPlexComputeIntegralFEM has at plexfem.c:2633)
-    PetscCall(DMPlexComputeBdIntegral(dm, X, label, numVals, vals, funcs, integral, ctx));
+    PetscCall(DMPlexComputeBdIntegral(dm, X, ownedLabel, numVals, vals, funcs, integral, ctx));
 
-    // Sum across all MPI ranks (mirrors DMPlexComputeIntegralFEM behaviour)
+    // --- MPI reduction (sum local owned contributions across all ranks) ---
     PetscScalar global_val;
     PetscCallMPI(MPIU_Allreduce(&integral[0], &global_val, 1, MPIU_SCALAR, MPIU_SUM,
                                 PetscObjectComm((PetscObject)dm)));
-
     *result = global_val;
 
+    PetscCall(DMLabelDestroy(&ownedLabel));
     PetscCall(PetscFree(funcs));
     PetscCall(PetscFree(integral));
 
diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py
index 3477194a..39595e2e 100644
--- a/tests/test_0502_boundary_integrals.py
+++ b/tests/test_0502_boundary_integrals.py
@@ -111,27 +111,37 @@ def test_bd_integral_invalid_boundary():
         uw.maths.BdIntegral(mesh, fn=1.0, boundary="Nonexistent")
 
 
-# --- Internal boundary tests ---
+# --- Internal boundary tests (BoxInternalBoundary) ---
+# These use lazy initialization because BoxInternalBoundary has a pre-existing
+# MPI bug (UnboundLocalError at cartesian.py:881) that would prevent the entire
+# test module from loading under mpirun.
 
 from underworld3.meshing import BoxInternalBoundary
 
-mesh_internal = BoxInternalBoundary(
-    minCoords=(0.0, 0.0),
-    maxCoords=(1.0, 1.0),
-    cellSize=1.0 / 32.0,
-    zintCoord=0.5,
-    simplex=True,
-)
+_mesh_internal = None
+_x_i = None
+_y_i = None
 
-x_i, y_i = mesh_internal.X
 
-# Need a mesh variable for PETSc
-_v_internal = uw.discretisation.MeshVariable("T_int", mesh_internal, 1, degree=2)
+def _get_internal_mesh():
+    global _mesh_internal, _x_i, _y_i
+    if _mesh_internal is None:
+        _mesh_internal = BoxInternalBoundary(
+            minCoords=(0.0, 0.0),
+            maxCoords=(1.0, 1.0),
+            cellSize=1.0 / 32.0,
+            zintCoord=0.5,
+            simplex=True,
+        )
+        _x_i, _y_i = _mesh_internal.X
+        uw.discretisation.MeshVariable("T_int", _mesh_internal, 1, degree=2)
+    return _mesh_internal, _x_i, _y_i
 
 
 def test_bd_integral_internal_boundary_length():
     """Internal boundary at y=0.5 across a unit box should have length 1.0."""
 
+    mesh_internal, _, _ = _get_internal_mesh()
     bd_int = uw.maths.BdIntegral(mesh_internal, fn=1.0, boundary="Internal")
     value = bd_int.evaluate()
 
@@ -141,6 +151,7 @@ def test_bd_integral_internal_boundary_length():
 def test_bd_integral_internal_coordinate_fn():
     """Integrate x along internal boundary at y=0.5: int_0^1 x dx = 0.5."""
 
+    mesh_internal, x_i, _ = _get_internal_mesh()
     bd_int = uw.maths.BdIntegral(mesh_internal, fn=x_i, boundary="Internal")
     value = bd_int.evaluate()
 
@@ -150,6 +161,7 @@ def test_bd_integral_internal_coordinate_fn():
 def test_bd_integral_internal_does_not_affect_external():
     """External boundaries should still work on the internal-boundary mesh."""
 
+    mesh_internal, _, _ = _get_internal_mesh()
     total = 0.0
     for bnd in ["Top", "Bottom", "Left", "Right"]:
         bd_int = uw.maths.BdIntegral(mesh_internal, fn=1.0, boundary=bnd)

From 4f7637577b422136e654ac6a592d1b8aaa069287 Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Wed, 4 Mar 2026 18:30:54 +1100
Subject: [PATCH 5/8] Add PETSc patch for ghost facet double-counting in
 boundary assembly

DMPlexComputeBdResidual_Internal and DMPlexComputeBdIntegral do not
exclude ghost facets (SF leaves) from the facet IS, unlike the volume
residual code which checks the ghost label. For internal boundaries
at partition junctions, this causes shared facets to be integrated on
multiple ranks, producing O(1) errors after MPI summation.

The patch filters ghost facets using ISDifference, following the same
pattern used by DMPlexComputeBdFaceGeomFVM for FVM face flux. It
applies cleanly to PETSc v3.18.0 through v3.24.5.

Includes:
- petsc-custom/patches/plexfem-ghost-facet-fix.patch
- build-petsc.sh integration (auto-applies after clone)
- MR description for upstream PETSc submission

Upstream MR branch pushed to gitlab.com/lmoresi/petsc
(fix/plexfem-ghost-facet-boundary-residual)

Underworld development team with AI support from Claude Code
---
 petsc-custom/build-petsc.sh                   | 26 ++++++
 petsc-custom/patches/PETSC_MR_DESCRIPTION.md  | 79 +++++++++++++++++++
 .../patches/plexfem-ghost-facet-fix.patch     | 59 ++++++++++++++
 3 files changed, 164 insertions(+)
 create mode 100644 petsc-custom/patches/PETSC_MR_DESCRIPTION.md
 create mode 100644 petsc-custom/patches/plexfem-ghost-facet-fix.patch

diff --git a/petsc-custom/build-petsc.sh b/petsc-custom/build-petsc.sh
index ee246c8d..ff8c94f7 100755
--- a/petsc-custom/build-petsc.sh
+++ b/petsc-custom/build-petsc.sh
@@ -57,6 +57,27 @@ clone_petsc() {
     echo "Clone complete."
 }
 
+apply_patches() {
+    echo "Applying UW3 patches to PETSc..."
+    cd "$PETSC_DIR"
+
+    # Fix ghost facet double-counting in boundary residual/integral assembly.
+    # Without this, internal boundary natural BCs and BdIntegral produce
+    # incorrect results in parallel (shared facets integrated on multiple ranks).
+    # Submitted upstream: [TODO: add PETSc MR link when available]
+    local patch="${SCRIPT_DIR}/patches/plexfem-ghost-facet-fix.patch"
+    if [ -f "$patch" ]; then
+        if git apply --check "$patch" 2>/dev/null; then
+            git apply "$patch"
+            echo "  Applied: plexfem-ghost-facet-fix.patch"
+        else
+            echo "  Skipped: plexfem-ghost-facet-fix.patch (already applied or conflict)"
+        fi
+    fi
+
+    echo "Patches complete."
+}
+
 configure_petsc() {
     echo "Configuring PETSc with AMR tools..."
     cd "$PETSC_DIR"
@@ -147,6 +168,7 @@ show_help() {
     echo "  build     Build PETSc"
     echo "  test      Run PETSc tests"
     echo "  petsc4py  Build and install petsc4py"
+    echo "  patch     Apply UW3 patches to PETSc source"
     echo "  clean     Remove PETSc directory"
     echo "  help      Show this help"
 }
@@ -155,6 +177,7 @@ show_help() {
 case "${1:-all}" in
     all)
         clone_petsc
+        apply_patches
         configure_petsc
         build_petsc
         build_petsc4py
@@ -175,6 +198,9 @@ case "${1:-all}" in
     build)
         build_petsc
         ;;
+    patch)
+        apply_patches
+        ;;
     test)
         test_petsc
         ;;
diff --git a/petsc-custom/patches/PETSC_MR_DESCRIPTION.md b/petsc-custom/patches/PETSC_MR_DESCRIPTION.md
new file mode 100644
index 00000000..bc796bae
--- /dev/null
+++ b/petsc-custom/patches/PETSC_MR_DESCRIPTION.md
@@ -0,0 +1,79 @@
+# PETSc Merge Request: Ghost Facet Fix
+
+## Submission details
+
+- **Fork**: https://gitlab.com/lmoresi/petsc
+- **Branch**: `fix/plexfem-ghost-facet-boundary-residual`
+- **Target**: `petsc/petsc` → `release`
+- **Create MR**: https://gitlab.com/lmoresi/petsc/-/merge_requests/new?merge_request%5Bsource_branch%5D=fix%2Fplexfem-ghost-facet-boundary-residual
+
+## Title
+
+DMPlex: filter ghost facets from boundary residual and integral assembly
+
+## Description
+
+### Problem
+
+`DMPlexComputeBdResidual_Internal` and `DMPlexComputeBdIntegral` construct
+`facetIS` as all facets at depth `dim-1` but do not exclude ghost facets
+(point SF leaves). The volume residual code in `DMPlexComputeResidualByKey`
+checks the `"ghost"` label before calling `DMPlexVecSetClosure` (line 5355),
+but the boundary residual path has no equivalent check.
+
+For **external** boundaries this is benign: ghost facets on the domain
+boundary still have `support[0]` pointing to a valid local cell, and
+`DMLocalToGlobal(ADD_VALUES)` correctly resolves ownership. But for
+**internal** boundaries — facets in the mesh interior that carry a
+`DMLabel` value and a `DM_BC_NATURAL` condition — facets at partition
+junctions are present on multiple ranks. Each rank independently
+integrates the boundary flux via `PetscFEIntegrateBdResidual`, and
+`DMLocalToGlobal(ADD_VALUES)` sums all contributions, causing
+double-counting.
+
+The same issue affects `DMPlexComputeBdIntegral`: ghost boundary facets
+are included in the per-face integral sum, which is then `MPI_Allreduce`d,
+again double-counting shared facets.
+
+### Use case: internal boundaries in geodynamics
+
+[Underworld3](https://github.com/underworldcode/underworld3) is a
+Python/PETSc geodynamics framework that uses DMPlex FEM throughout.
+Geophysical models commonly require traction (natural) boundary conditions
+on **internal surfaces** — for example, a fault plane or a material
+interface embedded within the mesh. These are set up by labeling interior
+facets with a `DMLabel` and registering `DM_BC_NATURAL` via
+`PetscDSAddBoundary` / `PetscWeakFormAddBdResidual`.
+
+This works correctly in serial but produces O(1) errors whenever the
+labeled internal boundary is split across partition boundaries, because
+shared facets are integrated on multiple ranks.
+
+### Fix
+
+After obtaining `facetIS` from `DMLabelGetStratumIS(depthLabel, dim-1, ...)`,
+filter out SF leaves using `ISDifference`. This ensures each boundary facet
+is processed by exactly one rank (the owner).
+
+This follows the same pattern already used in `DMPlexComputeBdFaceGeomFVM`
+(plexfem.c, around line 1087) which calls `PetscFindInt(face, nleaves, leaves, &loc)`
+to skip ghost faces during FVM face flux computation.
+
+### Testing
+
+Tested with Underworld3 Stokes solver using `DM_BC_NATURAL` on an internal
+boundary (annular mesh with labeled interior circle). Results match serial
+to 10+ significant figures across all processor counts tested (1, 2, 4, 8).
+Before the fix, any partition that splits the internal boundary produces
+O(1) errors from double-counted boundary contributions.
+
+The patch applies cleanly to all PETSc releases from v3.18.0 through
+v3.24.5 — the affected code has been structurally unchanged across these
+versions.
+
+## Reference examples
+
+- Stokes Green's function with internal boundary natural BCs:
+  https://github.com/underworldcode/underworld3/blob/development/docs/examples/fluid_mechanics/advanced/Ex_Stokes_Annulus_Kernels.py
+- Cartesian variant:
+  https://github.com/underworldcode/underworld3/blob/development/docs/examples/fluid_mechanics/advanced/Ex_Stokes_Cartesian_Kernels.py
diff --git a/petsc-custom/patches/plexfem-ghost-facet-fix.patch b/petsc-custom/patches/plexfem-ghost-facet-fix.patch
new file mode 100644
index 00000000..e1f7c8aa
--- /dev/null
+++ b/petsc-custom/patches/plexfem-ghost-facet-fix.patch
@@ -0,0 +1,59 @@
+diff --git a/src/dm/impls/plex/plexfem.c b/src/dm/impls/plex/plexfem.c
+index 0299f569070..26e6432585d 100644
+--- a/src/dm/impls/plex/plexfem.c
++++ b/src/dm/impls/plex/plexfem.c
+@@ -2874,6 +2874,26 @@ PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt num
+   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
+   PetscCall(DMGetDimension(dm, &dim));
+   PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS));
++  /* Filter out ghost facets (SF leaves) so that each boundary facet is only
++     counted on one rank. Without this, shared facets at partition boundaries
++     are integrated on multiple ranks, causing double-counting after MPI sum. */
++  if (facetIS) {
++    PetscSF  sf;
++    PetscInt nleaves;
++    const PetscInt *leaves;
++
++    PetscCall(DMGetPointSF(dm, &sf));
++    PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
++    if (nleaves > 0 && leaves) {
++      IS leafIS, ownedFacetIS;
++
++      PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS));
++      PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS));
++      PetscCall(ISDestroy(&leafIS));
++      PetscCall(ISDestroy(&facetIS));
++      facetIS = ownedFacetIS;
++    }
++  }
+   PetscCall(DMGetLocalSection(dm, &section));
+   PetscCall(PetscSectionGetNumFields(section, &Nf));
+   /* Get local solution with boundary values */
+@@ -5113,6 +5133,27 @@ static PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX
+   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
+   PetscCall(DMGetDimension(dm, &dim));
+   PetscCall(DMLabelGetStratumIS(depthLabel, dim - 1, &facetIS));
++  /* Filter out ghost facets (SF leaves) so that boundary residual contributions
++     from shared facets are only assembled on the owning rank. Without this,
++     internal boundary natural BCs at partition junctions get double-counted
++     because LocalToGlobal with ADD_VALUES sums contributions from all ranks. */
++  if (facetIS) {
++    PetscSF  sf;
++    PetscInt nleaves;
++    const PetscInt *leaves;
++
++    PetscCall(DMGetPointSF(dm, &sf));
++    PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL));
++    if (nleaves > 0 && leaves) {
++      IS leafIS, ownedFacetIS;
++
++      PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nleaves, leaves, PETSC_USE_POINTER, &leafIS));
++      PetscCall(ISDifference(facetIS, leafIS, &ownedFacetIS));
++      PetscCall(ISDestroy(&leafIS));
++      PetscCall(ISDestroy(&facetIS));
++      facetIS = ownedFacetIS;
++    }
++  }
+   PetscCall(PetscDSGetNumBoundary(prob, &numBd));
+   for (bd = 0; bd < numBd; ++bd) {
+     PetscWeakForm           wf;

From f4fb0378fddbecd1ead8917deda0b7c5f2a6d182 Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Wed, 4 Mar 2026 21:24:30 +1100
Subject: [PATCH 6/8] Fix DMPlexComputeBdIntegral signature for PETSc < 3.22

PETSc changed DMPlexComputeBdIntegral from void (*func)(...) to
void (**funcs)(...) in v3.22.0. CI uses PETSc 3.21.5 (pinned in
environment.yaml), so the wrapper must handle both signatures.

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_compat.h | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 649cc512..78989dbf 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -211,7 +211,14 @@ PetscErrorCode UW_DMPlexComputeBdIntegral(DM dm, Vec X,
     PetscScalar *integral;
     PetscCall(PetscCalloc1(Nf, &integral));
 
+    // PETSc changed DMPlexComputeBdIntegral signature in v3.22.0:
+    //   <= 3.21.x: void (*func)(...)     — single function pointer
+    //   >= 3.22.0: void (**funcs)(...)    — array of Nf function pointers
+#if PETSC_VERSION_GE(3, 22, 0)
     PetscCall(DMPlexComputeBdIntegral(dm, X, ownedLabel, numVals, vals, funcs, integral, ctx));
+#else
+    PetscCall(DMPlexComputeBdIntegral(dm, X, ownedLabel, numVals, vals, funcs[0], integral, ctx));
+#endif
 
     // --- MPI reduction (sum local owned contributions across all ranks) ---
     PetscScalar global_val;

From 2881718dc529f1632826e5c54fa87deed430b9f6 Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Wed, 4 Mar 2026 21:38:18 +1100
Subject: [PATCH 7/8] Skip BoxInternalBoundary tests under MPI (pre-existing
 mesh constructor bug)

Underworld development team with AI support from Claude Code
---
 tests/test_0502_boundary_integrals.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py
index 39595e2e..e59437ef 100644
--- a/tests/test_0502_boundary_integrals.py
+++ b/tests/test_0502_boundary_integrals.py
@@ -112,9 +112,9 @@ def test_bd_integral_invalid_boundary():
 
 
 # --- Internal boundary tests (BoxInternalBoundary) ---
-# These use lazy initialization because BoxInternalBoundary has a pre-existing
-# MPI bug (UnboundLocalError at cartesian.py:881) that would prevent the entire
-# test module from loading under mpirun.
+# BoxInternalBoundary has a pre-existing MPI bug (UnboundLocalError in the mesh
+# constructor) so these tests are skipped under MPI. They use lazy initialization
+# to avoid crashing the entire module if the mesh constructor fails.
 
 from underworld3.meshing import BoxInternalBoundary
 
@@ -138,6 +138,7 @@ def _get_internal_mesh():
     return _mesh_internal, _x_i, _y_i
 
 
+@pytest.mark.skipif(uw.mpi.size > 1, reason="BoxInternalBoundary has pre-existing MPI bug")
 def test_bd_integral_internal_boundary_length():
     """Internal boundary at y=0.5 across a unit box should have length 1.0."""
 
@@ -148,6 +149,7 @@ def test_bd_integral_internal_boundary_length():
     assert abs(value - 1.0) < 0.001, f"Expected 1.0, got {value}"
 
 
+@pytest.mark.skipif(uw.mpi.size > 1, reason="BoxInternalBoundary has pre-existing MPI bug")
 def test_bd_integral_internal_coordinate_fn():
     """Integrate x along internal boundary at y=0.5: int_0^1 x dx = 0.5."""
 
@@ -158,6 +160,7 @@ def test_bd_integral_internal_coordinate_fn():
     assert abs(value - 0.5) < 0.01, f"Expected 0.5, got {value}"
 
 
+@pytest.mark.skipif(uw.mpi.size > 1, reason="BoxInternalBoundary has pre-existing MPI bug")
 def test_bd_integral_internal_does_not_affect_external():
     """External boundaries should still work on the internal-boundary mesh."""
 

From 220fd1cc67626e477eac4e99a7c8ae081504e72e Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Fri, 6 Mar 2026 21:28:11 +1100
Subject: [PATCH 8/8] Guard against NULL boundary labels and add internal
 normal vector tests

Address PR #70 review feedback:
- C wrapper handles NULL DMLabel gracefully (contributes 0 to MPI reduction)
- Cython passes NULL instead of segfaulting when getLabel returns None
- Add 5 new tests for normal vector integrals on internal boundaries
  (box: n_y, n_x, x*n_y; annulus: |n|^2, n.t_hat)

Underworld development team with AI support from Claude Code
---
 src/underworld3/cython/petsc_compat.h  |  9 +++
 src/underworld3/cython/petsc_maths.pyx | 13 +++-
 tests/test_0502_boundary_integrals.py  | 84 ++++++++++++++++++++++++++
 3 files changed, 104 insertions(+), 2 deletions(-)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 78989dbf..bd6e49e3 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -163,6 +163,15 @@ PetscErrorCode UW_DMPlexComputeBdIntegral(DM dm, Vec X,
     PetscCall(DMGetLocalSection(dm, &section));
     PetscCall(PetscSectionGetNumFields(section, &Nf));
 
+    // If label is NULL (boundary not present on this rank), contribute 0
+    // but still participate in the MPI Allreduce to avoid hangs.
+    if (!label) {
+        PetscScalar zero = 0.0;
+        PetscCallMPI(MPIU_Allreduce(&zero, result, 1, MPIU_SCALAR, MPIU_SUM,
+                                    PetscObjectComm((PetscObject)dm)));
+        PetscFunctionReturn(PETSC_SUCCESS);
+    }
+
     // --- Build a ghost-point bitset from the point SF ---
     PetscSF   sf;
     PetscInt  nleaves;
diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx
index 529f68ca..da12acdf 100644
--- a/src/underworld3/cython/petsc_maths.pyx
+++ b/src/underworld3/cython/petsc_maths.pyx
@@ -406,15 +406,24 @@ class BdIntegral:
         cdef PetscInt num_vals = 1
 
         c_label = mesh.dm.getLabel(self.boundary)
-        cdef DMLabel dmlabel = c_label
 
         # Output value
         cdef PetscScalar result = 0.0
 
+        # Label can be None if this boundary has no facets on this process
+        # (normal in parallel) or if the DM doesn't carry the label at all.
+        # The C wrapper handles NULL labels gracefully (contributes 0 to the
+        # MPI reduction so all ranks still participate).
+        cdef PetscDMLabel c_dmlabel = NULL
+        cdef DMLabel dmlabel
+        if c_label:
+            dmlabel = c_label
+            c_dmlabel = dmlabel.dmlabel
+
         # Call the boundary integral
         ierr = UW_DMPlexComputeBdIntegral(
             dm_c.dm, cgvec.vec,
-            dmlabel.dmlabel, num_vals, &label_val,
+            c_dmlabel, num_vals, &label_val,
             ext.fns_bd_residual[0],
             &result, NULL
         )
diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py
index e59437ef..e7a6f05a 100644
--- a/tests/test_0502_boundary_integrals.py
+++ b/tests/test_0502_boundary_integrals.py
@@ -160,6 +160,54 @@ def test_bd_integral_internal_coordinate_fn():
     assert abs(value - 0.5) < 0.01, f"Expected 0.5, got {value}"
 
 
+@pytest.mark.skipif(uw.mpi.size > 1, reason="BoxInternalBoundary has pre-existing MPI bug")
+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)."""
+
+    mesh_internal, _, _ = _get_internal_mesh()
+    Gamma = mesh_internal.Gamma
+    n_y = 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}"
+
+
+@pytest.mark.skipif(uw.mpi.size > 1, reason="BoxInternalBoundary has pre-existing MPI bug")
+def test_bd_integral_internal_normal_nx():
+    """Integrate n_x along internal boundary at y=0.5.
+    The internal boundary is horizontal, so n_x should be ~0."""
+
+    mesh_internal, _, _ = _get_internal_mesh()
+    Gamma = mesh_internal.Gamma
+    n_x = Gamma[0]
+
+    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}"
+
+
+@pytest.mark.skipif(uw.mpi.size > 1, reason="BoxInternalBoundary has pre-existing MPI bug")
+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."""
+
+    mesh_internal, x_i, _ = _get_internal_mesh()
+    Gamma = mesh_internal.Gamma
+    n_y = 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}"
+
+
 @pytest.mark.skipif(uw.mpi.size > 1, reason="BoxInternalBoundary has pre-existing MPI bug")
 def test_bd_integral_internal_does_not_affect_external():
     """External boundaries should still work on the internal-boundary mesh."""
@@ -207,3 +255,39 @@ def test_bd_integral_annulus_outer_circumference():
     expected = 2 * np.pi * 1.5
 
     assert abs(value - expected) < 0.01, f"Expected {expected:.4f}, got {value}"
+
+
+def test_bd_integral_annulus_internal_normal_squared():
+    """Normal magnitude squared integrated over internal circle at r=1.0.
+
+    Since n is a unit normal, n_x^2 + n_y^2 = 1 everywhere on the boundary.
+    Integrating 1 over the circle gives the circumference 2*pi*R.
+    (We use n_x^2 + n_y^2 rather than n dot r_hat because internal boundary
+    normals have arbitrary per-facet orientation that may cancel.)"""
+
+    Gamma = mesh_annulus.Gamma
+    n_sq = Gamma[0]**2 + Gamma[1]**2
+
+    bd_int = uw.maths.BdIntegral(mesh_annulus, fn=n_sq, boundary="Internal")
+    value = bd_int.evaluate()
+    expected = 2 * np.pi * _R_INTERNAL
+
+    assert abs(value - expected) < 0.05, f"Expected {expected:.4f}, got {value}"
+
+
+def test_bd_integral_annulus_internal_normal_tangential():
+    """Tangential component of normal integrated over internal circle should be ~0.
+
+    The tangential direction is t = (-y/r, x/r). Since n is radial,
+    n.t should integrate to zero."""
+
+    x_a, y_a = mesh_annulus.X
+    Gamma = mesh_annulus.Gamma
+    r = sympy.sqrt(x_a**2 + y_a**2)
+    # n dot t_hat = (-n_x * y + n_y * x) / r
+    n_dot_that = (-Gamma[0] * y_a + Gamma[1] * x_a) / r
+
+    bd_int = uw.maths.BdIntegral(mesh_annulus, fn=n_dot_that, boundary="Internal")
+    value = bd_int.evaluate()
+
+    assert abs(value) < 0.05, f"Expected ~0, got {value}"