Add region labels to BoxInternalBoundary via centroid classification

BoxInternalBoundary uses a different label import path (Face Sets +
_dm_unstack_bcs) that conflicts with useRegions=True. Instead of
changing the import mechanism, classify cells by centroid position
after mesh construction: cells below zintCoord = Inner, above = Outer.

Works for both 2D and 3D, simplex and structured meshes.
extract_region("Inner") verified on 2D box.
All 21 boundary integral tests pass.

Underworld development team with AI support from Claude Code

Author: lmoresi

src/underworld3/meshing/cartesian.py

@@ -906,12 +906,49 @@ def box_return_coords_to_bounds(coords):
     )
     uw.adaptivity._dm_unstack_bcs(new_mesh.dm, new_mesh.boundaries, "Face Sets")
 
-    # Assign regions
+    # Create region labels by classifying cells based on centroid position
+    # relative to the internal boundary coordinate
     if dim == 2:
         new_mesh.regions = regions_2D
     else:
         new_mesh.regions = regions_3D
 
+    dm = new_mesh.dm
+    depth_label = dm.getLabel("depth")
+    cell_is = depth_label.getStratumIS(dim)
+
+    if cell_is:
+        cells = cell_is.getIndices()
+        coord_sec = dm.getCoordinateSection()
+        coord_vec = dm.getCoordinatesLocal()
+        coord_arr = coord_vec.array
+
+        for region in new_mesh.regions:
+            dm.createLabel(region.name)
+
+        inner_label = dm.getLabel(new_mesh.regions.Inner.name)
+        outer_label = dm.getLabel(new_mesh.regions.Outer.name)
+
+        # z-coordinate index: 1 for 2D (y), 2 for 3D (z)
+        z_idx = dim - 1
+
+        for cell in cells:
+            # Compute centroid from cell vertex coordinates
+            closure = dm.getTransitiveClosure(cell)[0]
+            vert_coords = []
+            for pt in closure:
+                ndof = coord_sec.getDof(pt)
+                if ndof > 0:
+                    off = coord_sec.getOffset(pt)
+                    vert_coords.append(coord_arr[off + z_idx])
+
+            if vert_coords:
+                centroid_z = sum(vert_coords) / len(vert_coords)
+                if centroid_z < zintCoord:
+                    inner_label.setValue(cell, new_mesh.regions.Inner.value)
+                else:
+                    outer_label.setValue(cell, new_mesh.regions.Outer.value)
+
     return new_mesh