Fix cell-family detection on empty MPI ranks

docs/developer/subsystems/meshing.md

@@ -33,6 +33,24 @@ The meshing subsystem handles computational mesh generation and manipulation for
 - QuadBox / HexBox            # Structured meshes
 ```
 
+## Empty MPI Partitions
+
+Use `mesh.isSimplex` for the mesh-wide cell family. PETSc's
+`mesh.dm.isSimplex()` is a rank-local query and returns `False` on a rank
+with no cells, even when the distributed mesh consists of triangles or
+tetrahedra. UW3 infers the family collectively from populated ranks before
+constructing coordinate finite elements or element metadata.
+
+This is also an MPI correctness requirement: constructing simplex and
+tensor-product coordinate elements on the same communicator consumes
+different PETSc message tags. The first mesh may appear to construct
+successfully, but a later HDF5 boundary-label load can deadlock in
+`PetscSFSetUp_Basic`. This failure does not require a transport solver.
+
+`tests/parallel/test_0781_empty_rank_mesh_sequence.py` covers triangles,
+tetrahedra, quadrilaterals and hexahedra with deliberately empty partitions,
+then verifies a second mesh's volume and boundary integrals using P2 data.
+
 ## Documentation Needs
 
 ### Critical Gaps
@@ -71,4 +89,4 @@ This section needs:
 
 ---
 
-*This document serves as a placeholder for comprehensive meshing system documentation.*
+*This document serves as a placeholder for comprehensive meshing system documentation.*

src/underworld3/discretisation/discretisation_mesh.py

@@ -589,9 +589,19 @@ def __init__(
         self._setup_symbolic_coordinates(coordinate_system_type)
 
         try:
-            self.isSimplex = self.dm.isSimplex()
+            local_simplex = self.dm.isSimplex()
         except:
-            self.isSimplex = simplex
+            local_simplex = simplex
+
+        # DMPlexIsSimplex is rank-local and returns False on empty ranks.
+        # Coordinate FE construction must use one cell family everywhere;
+        # mixed simplex/tensor construction desynchronises PETSc MPI tags.
+        cell_start, cell_end = self.dm.getHeightStratum(0)
+        cell_families = self.dm.comm.tompi4py().allgather(
+            local_simplex if cell_end > cell_start else None
+        )
+        populated_families = [family for family in cell_families if family is not None]
+        self.isSimplex = all(populated_families) if populated_families else simplex
 
         # Using WeakValueDictionary to prevent circular references
         self._vars = weakref.WeakValueDictionary()
@@ -650,7 +660,7 @@ class ElementInfo:
             entities: tuple
             face_entities: tuple
 
-        if self.dm.isSimplex():
+        if self.isSimplex:
             if self.dim == 2:
                 self._element = ElementInfo("triangle", (1, 3, 3), (0, 1, 2))
             else:

tests/parallel/test_0781_empty_rank_mesh_sequence.py

@@ -0,0 +1,53 @@
+"""Empty partitions must not change the cell family or poison later mesh loads.
+
+The original failure appeared in SUPG test 1077 after the empty-partition
+test: rank-local DMPlexIsSimplex returned False on empty ranks, so the
+coordinate FE consumed different COMM_WORLD tags. The next HDF5 label load
+then hung in PetscSFSetUp_Basic. No transport solve is needed to reproduce it.
+"""
+
+import numpy as np
+import pytest
+from petsc4py import PETSc
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_2,
+              pytest.mark.tier_b, pytest.mark.timeout(120)]
+
+
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("simplex", [True, False])
+def test_cell_family_and_next_mesh_load_on_empty_ranks(dim, simplex):
+    if simplex:
+        coords = np.vstack([np.zeros(dim), np.eye(dim)])
+        cells = np.arange(dim + 1, dtype=PETSc.IntType).reshape(1, -1)
+        dm = PETSc.DMPlex().createFromCellList(dim, cells, coords)
+    else:
+        dm = PETSc.DMPlex().createBoxMesh([1] * dim, simplex=False)
+    first = uw.discretisation.Mesh(dm, simplex=simplex, qdegree=3)
+    start, end = first.dm.getHeightStratum(0)
+    counts = uw.mpi.comm.allgather(end - start)
+    assert min(counts) == 0 and sum(counts) == 1, counts
+
+    families = uw.mpi.comm.allgather(first.isSimplex)
+    assert families == [simplex] * uw.mpi.size, families
+    expected = { (2, True): "triangle", (3, True): "tetrahedron",
+                 (2, False): "quadrilateral", (3, False): "hexahedron" }
+    elements = uw.mpi.comm.allgather(first._element.type)
+    assert elements == [expected[dim, simplex]] * uw.mpi.size, elements
+
+    # This public constructor exercises the HDF5 label SF exchange that hung.
+    second = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim,
+        cellSize=0.25, qdegree=3, regular=False,
+    )
+    field = uw.discretisation.MeshVariable("T", second, 1, degree=2)
+    field.array[:, 0, 0] = 1.0
+    volume = uw.maths.Integral(second, field.sym[0]).evaluate()
+    assert np.isclose(volume, 1.0, rtol=1e-12, atol=1e-12), volume
+    for boundary in second.boundaries:
+        if boundary.name in ("Null_Boundary", "All_Boundaries"):
+            continue
+        area = uw.maths.BdIntegral(second, field.sym[0], boundary=boundary.name).evaluate()
+        assert np.isclose(area, 1.0, rtol=1e-12, atol=1e-12), (boundary.name, area)