From d673594047afafc86c58989faee7c64d22474568 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Mon, 30 Mar 2026 12:13:54 +1100
Subject: [PATCH] Handle empty boundary strata in MPI BdIntegral

---
 src/underworld3/cython/petsc_compat.h         | 17 ++++--
 ...est_0765_internal_boundary_integral_mpi.py | 54 +++++++++++++++++++
 2 files changed, 68 insertions(+), 3 deletions(-)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index f2e5cc103..2ae23d60a 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -173,15 +173,26 @@ PetscErrorCode UW_DMPlexComputeBdIntegral(DM dm, Vec X,
 {
     PetscSection  section;
     PetscInt      Nf;
+    PetscInt      localCount = 0;
 
     PetscFunctionBeginUser;
 
     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) {
+    // If the label is NULL or the requested boundary has no local entities on
+    // this rank, contribute 0 but still participate in the MPI Allreduce to
+    // avoid hangs. Parallel DMPlex boundary assembly can deadlock if some
+    // ranks enter with an empty local stratum.
+    if (label) {
+        for (PetscInt i = 0; i < numVals; ++i) {
+            PetscInt stratumSize = 0;
+            PetscCall(DMLabelGetStratumSize(label, vals[i], &stratumSize));
+            localCount += stratumSize;
+        }
+    }
+
+    if (!label || localCount == 0) {
         PetscScalar zero = 0.0;
         PetscCallMPI(MPIU_Allreduce(&zero, result, 1, MPIU_SCALAR, MPIU_SUM,
                                     PetscObjectComm((PetscObject)dm)));
diff --git a/tests/parallel/test_0765_internal_boundary_integral_mpi.py b/tests/parallel/test_0765_internal_boundary_integral_mpi.py
index 17ed49e24..0d229120d 100644
--- a/tests/parallel/test_0765_internal_boundary_integral_mpi.py
+++ b/tests/parallel/test_0765_internal_boundary_integral_mpi.py
@@ -5,6 +5,7 @@
 ghost/internal facet handling in PETSc boundary assembly paths.
 """
 
+import math
 import numpy as np
 import pytest
 import underworld3 as uw
@@ -59,3 +60,56 @@ def test_outer_boundary_circumference_parallel():
 
     rel_err = abs(value - expected) / expected
     assert rel_err < 2.0e-2, f"Outer circumference rel_err={rel_err:.3e}, value={value}, expected={expected}"
+
+
+@pytest.mark.mpi(min_size=2)
+def test_deformed_spherical_shell_boundary_area_parallel():
+    """
+    Boundary integrals must remain valid after coordinate deformation in MPI.
+
+    This specifically guards the BdIntegral path when some ranks have no local
+    entities for the requested boundary after mesh deformation.
+    """
+
+    mesh = uw.meshing.SphericalShell(
+        radiusOuter=1.0,
+        radiusInner=0.5,
+        cellSize=1.0 / 4.0,
+        degree=1,
+        qdegree=2,
+    )
+    uw.discretisation.MeshVariable("T_spherical_deformed_bd", mesh, 1, degree=1, continuous=True)
+
+    coords = np.asarray(mesh.X.coords, dtype=np.float64).copy()
+    radii = np.linalg.norm(coords, axis=1)
+    thickness = 0.5
+    t = (radii - 0.5) / thickness
+    a = math.log(2.0)
+    mapped = (np.exp(a * t) - 1.0) / (math.exp(a) - 1.0)
+    new_radii = 0.5 + thickness * mapped
+    mesh._deform_mesh(coords * (new_radii / radii)[:, None])
+
+    lower = float(uw.maths.BdIntegral(mesh=mesh, fn=1.0, boundary="Lower").evaluate())
+    upper = float(uw.maths.BdIntegral(mesh=mesh, fn=1.0, boundary="Upper").evaluate())
+
+    expected_lower = 4.0 * math.pi * 0.5**2
+    expected_upper = 4.0 * math.pi * 1.0**2
+
+    rel_err_lower = abs(lower - expected_lower) / expected_lower
+    rel_err_upper = abs(upper - expected_upper) / expected_upper
+
+    assert rel_err_lower < 5.0e-2, (
+        f"Deformed lower area rel_err={rel_err_lower:.3e}, value={lower}, expected={expected_lower}"
+    )
+    assert rel_err_upper < 5.0e-2, (
+        f"Deformed upper area rel_err={rel_err_upper:.3e}, value={upper}, expected={expected_upper}"
+    )
+
+    gathered_lower = uw.mpi.comm.allgather(lower)
+    gathered_upper = uw.mpi.comm.allgather(upper)
+    assert max(gathered_lower) - min(gathered_lower) < 1.0e-12, (
+        f"Rank mismatch in lower-boundary integral values: {gathered_lower}"
+    )
+    assert max(gathered_upper) - min(gathered_upper) < 1.0e-12, (
+        f"Rank mismatch in upper-boundary integral values: {gathered_upper}"
+    )