From fe689869be371ebc39cffeb97215e1a707f4ac44 Mon Sep 17 00:00:00 2001 From: Hussam Date: Thu, 13 Aug 2026 14:06:06 +0100 Subject: [PATCH 01/11] Make serial submesh construction linear in the number of points Submesh(mesh, comm=COMM_SELF) cost O(P^2) in the number of mesh points P. On three ranks a UnitSquareMesh(181, 181) needed 8.9 seconds to build one serial submesh. Four times more cells cost about ten times more time. Every preconditioner that builds a subdomain matrix through firedrake.preconditioners.matis pays this cost, because local_mesh builds one serial submesh per process. The cost was in submesh_correct_entity_classes. A submesh on COMM_SELF has no neighbour, so it has no ghost points and no owned points, and every point is core. The function set that state one point at a time. DMLabelSetValue moves the stratum into a hash set and destroys its index set. The DMLabelHasPoint of the next iteration reads the hash set back, sorts it, and rebuilds the index set and the bit array. Each point therefore cost O(k log k) in the size k of the core stratum so far. Set the three strata in bulk instead. This is the same fix as the one in PR #5326 for the coarse and fine cell maps. submesh_create also built a temporary label that marked every cell, and passed it to DMPlexFilter. DMPlexFilter selects every cell by itself when it gets no label, so the label was unnecessary. Submesh now leaves label_name as None for a codim-0 submesh of the whole mesh. A codim-1 submesh still passes the "depth" label, because DMPlexFilter with no label selects cells, not facets. Measured on three ranks, for submesh_correct_entity_classes alone: cells before after 717 34 ms 0.12 ms 2787 320 ms 0.43 ms 11059 3214 ms 1.85 ms 22023 8907 ms 3.80 ms The whole submesh build at 11059 cells falls from 3256 ms to 71 ms. Both paths produce the same submesh. On UnitSquareMesh(128, 128) with three ranks, the submesh built with a label and the submesh built without one have the same chart, the same subpoint index set, and the same pyop2_core, pyop2_owned, pyop2_ghost and Face Sets strata. Co-Authored-By: Claude Opus 5 --- firedrake/cython/dmcommon.pyx | 106 ++++++++++--------- firedrake/cython/petschdr.pxi | 2 + firedrake/mesh.py | 8 +- tests/firedrake/submesh/test_submesh_comm.py | 21 ++++ 4 files changed, 85 insertions(+), 52 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 54774b456f..00c120c3e4 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4001,10 +4001,10 @@ def submesh_create(PETSc.DM dm, DMPlex representing the mesh topology subdim : int Topological dimension of the submesh - label_name : str - Name of the label - subdomain_id : int | Sequence - Values in the label + label_name : str | None + Name of the label, or `None` to select every cell + subdomain_id : int | Sequence | None + Values in the label, unused if ``label_name`` is `None` ignore_label_halo : bool If labeled points in the halo are ignored. comm : PETSc.Comm | None @@ -4018,38 +4018,44 @@ def submesh_create(PETSc.DM dm, PetscInt pStart, pEnd, p, i, stratum_size = 0, label_value = 1 const PetscInt *stratum_indices = NULL - # Cast subdomain_id into an iterable - if isinstance(subdomain_id, str) or not isinstance(subdomain_id, Sequence): - subdomain_id = (subdomain_id,) - # Take the union of the all the label values - label = dm.getLabel(label_name) - points = PETSc.IS() - for sub in subdomain_id: - if isinstance(sub, Integral): - subpoints = label.getStratumIS(sub) - elif sub == "on_boundary": - subpoints = dm.getStratumIS("exterior_facets", 1) - else: - raise ValueError(f"Submesh construction got invalid subdomain_id {sub}.") + if label_name is None: + # Every cell is wanted. DMPlexFilter selects all of them when given no + # label, so building one that marks the whole mesh only to hand it back + # would be wasted work. + temp_label = None + else: + # Cast subdomain_id into an iterable + if isinstance(subdomain_id, str) or not isinstance(subdomain_id, Sequence): + subdomain_id = (subdomain_id,) + # Take the union of the all the label values + label = dm.getLabel(label_name) + points = PETSc.IS() + for sub in subdomain_id: + if isinstance(sub, Integral): + subpoints = label.getStratumIS(sub) + elif sub == "on_boundary": + subpoints = dm.getStratumIS("exterior_facets", 1) + else: + raise ValueError(f"Submesh construction got invalid subdomain_id {sub}.") + if points: + points = points.union(subpoints) + else: + points = subpoints + # Create temp_label that contains no lower-dimensional points. + dm.createLabel(temp_label_name) + temp_label = dm.getLabel(temp_label_name) if points: - points = points.union(subpoints) - else: - points = subpoints - # Create temp_label that contains no lower-dimensional points. - dm.createLabel(temp_label_name) - temp_label = dm.getLabel(temp_label_name) - if points: - CHKERR(ISGetSize(points.iset, &stratum_size)) - if stratum_size > 0: - CHKERR(ISGetIndices(points.iset, &stratum_indices)) - CHKERR(DMPlexGetDepthStratum(dm.dm, subdim, &pStart, &pEnd)) - for i in range(stratum_size): - p = stratum_indices[i] - # Only include points on the submesh topological dimension, - # culling all lower-dimensional points. - if pStart <= p < pEnd: - CHKERR(DMLabelSetValue(temp_label.dmlabel, p, label_value)) - CHKERR(ISRestoreIndices(points.iset, &stratum_indices)) + CHKERR(ISGetSize(points.iset, &stratum_size)) + if stratum_size > 0: + CHKERR(ISGetIndices(points.iset, &stratum_indices)) + CHKERR(DMPlexGetDepthStratum(dm.dm, subdim, &pStart, &pEnd)) + for i in range(stratum_size): + p = stratum_indices[i] + # Only include points on the submesh topological dimension, + # culling all lower-dimensional points. + if pStart <= p < pEnd: + CHKERR(DMLabelSetValue(temp_label.dmlabel, p, label_value)) + CHKERR(ISRestoreIndices(points.iset, &stratum_indices)) # Make submesh using temp_label. subdm, ownership_transfer_sf = dm.filter(label=temp_label, value=label_value, @@ -4057,8 +4063,9 @@ def submesh_create(PETSc.DM dm, sanitizeSubMesh=PETSC_TRUE, comm=comm) # Destroy temp_label. - dm.removeLabel(temp_label_name) - subdm.removeLabel(temp_label_name) + if temp_label is not None: + dm.removeLabel(temp_label_name) + subdm.removeLabel(temp_label_name) submesh_update_facet_labels(dm, subdm) submesh_correct_entity_classes(dm, subdm, ownership_transfer_sf) return subdm @@ -4087,6 +4094,7 @@ def submesh_correct_entity_classes(PETSc.DM dm, const PetscInt *ilocal = NULL const PetscSFNode *iremote = NULL PETSc.IS subpoint_is + PETSc.IS all_points const PetscInt *subpoint_indices = NULL np.ndarray ownership_loss np.ndarray ownership_gain @@ -4108,18 +4116,18 @@ def submesh_correct_entity_classes(PETSc.DM dm, CHKERR(DMLabelCreateIndex(lbl_ghost, subpStart, subpEnd)) if subdm.comm.size == 1: - # Undistributed case: relabel every point as core - for subp in range(subpStart, subpEnd): - CHKERR(DMLabelHasPoint(lbl_core, subp, &has)) - if has: - continue - CHKERR(DMLabelHasPoint(lbl_ghost, subp, &has)) - if has: - CHKERR(DMLabelClearValue(lbl_ghost, subp, 1)) - CHKERR(DMLabelHasPoint(lbl_owned, subp, &has)) - if has: - CHKERR(DMLabelClearValue(lbl_owned, subp, 1)) - CHKERR(DMLabelSetValue(lbl_core, subp, 1)) + # Undistributed case: relabel every point as core. Setting the strata + # in bulk keeps this linear in the number of points: DMLabelSetValue + # invalidates the label index, which the next DMLabelHasPoint would + # rebuild and re-sort, so relabelling point by point costs O(n log n) + # each time round. + all_points = PETSc.IS().createStride(subpEnd - subpStart, + first=subpStart, step=1, + comm=PETSc.COMM_SELF) + CHKERR(DMLabelClearStratum(lbl_owned, 1)) + CHKERR(DMLabelClearStratum(lbl_ghost, 1)) + CHKERR(DMLabelSetStratumIS(lbl_core, 1, (all_points).iset)) + all_points.destroy() else: ownership_loss = np.zeros(pEnd - pStart, dtype=IntType) ownership_gain = np.zeros(pEnd - pStart, dtype=IntType) diff --git a/firedrake/cython/petschdr.pxi b/firedrake/cython/petschdr.pxi index 42ac97e24d..36a69fbd3b 100644 --- a/firedrake/cython/petschdr.pxi +++ b/firedrake/cython/petschdr.pxi @@ -95,6 +95,8 @@ cdef extern from "petscdmlabel.h" nogil: PetscErrorCode DMLabelClearValue(DMLabel, PetscInt, PetscInt) PetscErrorCode DMLabelGetStratumSize(DMLabel, PetscInt, PetscInt*) PetscErrorCode DMLabelGetStratumIS(DMLabel, PetscInt, PETSc.PetscIS*) + PetscErrorCode DMLabelSetStratumIS(DMLabel, PetscInt, PETSc.PetscIS) + PetscErrorCode DMLabelClearStratum(DMLabel, PetscInt) cdef extern from "petscdm.h" nogil: PetscErrorCode DMCreateLabel(PETSc.PetscDM,char[]) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 881cdfe523..94f3559318 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -5019,9 +5019,11 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig if subdomain_id is None: if label_name is not None: raise ValueError("subdomain_id=None requires label_name=None.") - # Select all entities - label_name = "depth" - subdomain_id = subdim + if subdim != dim: + # Select all entities of the submesh dimension. A codim-0 submesh + # is every cell, which submesh_create selects without a label. + label_name = "depth" + subdomain_id = subdim elif label_name is None: if subdim == dim: label_name = dmcommon.CELL_SETS_LABEL diff --git a/tests/firedrake/submesh/test_submesh_comm.py b/tests/firedrake/submesh/test_submesh_comm.py index 0a442179a9..ef1695f369 100644 --- a/tests/firedrake/submesh/test_submesh_comm.py +++ b/tests/firedrake/submesh/test_submesh_comm.py @@ -44,6 +44,27 @@ def test_create_submesh_comm_self(reorder, ignore_halo): assert np.allclose(submesh.coordinates.dat.data_ro, x.dat.data_ro) +@pytest.mark.parallel([1, 3]) +@pytest.mark.parametrize("ignore_halo", [False, True]) +def test_submesh_comm_self_entity_classes(ignore_halo): + """A submesh on COMM_SELF must own every point that it holds. + + The parent mesh divides its points into the pyop2 classes core, owned and + ghost. A submesh on COMM_SELF has no neighbour, so it has no ghost points + and no owned points either. Every point is core. + """ + mesh = UnitSquareMesh( + 8, 8, distribution_parameters={ + "overlap_type": (DistributedMeshOverlapType.VERTEX, 1)}) + submesh = Submesh(mesh, ignore_halo=ignore_halo, comm=COMM_SELF) + + plex = submesh.topology_dm + pStart, pEnd = plex.getChart() + assert plex.getStratumSize("pyop2_core", 1) == pEnd - pStart + assert plex.getStratumSize("pyop2_owned", 1) == 0 + assert plex.getStratumSize("pyop2_ghost", 1) == 0 + + @pytest.mark.parallel([1, 3]) @pytest.mark.parametrize("family,degree", [("DG", 0), ("CG", 1)]) @pytest.mark.parametrize("reorder", [False, True]) From 3ec9bc84396d805221a7aef054e70a3774a708a6 Mon Sep 17 00:00:00 2001 From: Hussam Date: Thu, 13 Aug 2026 16:59:03 +0100 Subject: [PATCH 02/11] Edit claude comments manually --- firedrake/cython/dmcommon.pyx | 10 ++-------- firedrake/mesh.py | 4 ++-- 2 files changed, 4 insertions(+), 10 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 00c120c3e4..7005230972 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4019,9 +4019,7 @@ def submesh_create(PETSc.DM dm, const PetscInt *stratum_indices = NULL if label_name is None: - # Every cell is wanted. DMPlexFilter selects all of them when given no - # label, so building one that marks the whole mesh only to hand it back - # would be wasted work. + # label=None covers all cells. DMPlexFilter already handles this. temp_label = None else: # Cast subdomain_id into an iterable @@ -4116,11 +4114,7 @@ def submesh_correct_entity_classes(PETSc.DM dm, CHKERR(DMLabelCreateIndex(lbl_ghost, subpStart, subpEnd)) if subdm.comm.size == 1: - # Undistributed case: relabel every point as core. Setting the strata - # in bulk keeps this linear in the number of points: DMLabelSetValue - # invalidates the label index, which the next DMLabelHasPoint would - # rebuild and re-sort, so relabelling point by point costs O(n log n) - # each time round. + # Undistributed case: relabel every point as core all_points = PETSc.IS().createStride(subpEnd - subpStart, first=subpStart, step=1, comm=PETSc.COMM_SELF) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 94f3559318..2070ddeccd 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -5020,8 +5020,8 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig if label_name is not None: raise ValueError("subdomain_id=None requires label_name=None.") if subdim != dim: - # Select all entities of the submesh dimension. A codim-0 submesh - # is every cell, which submesh_create selects without a label. + # DMPlexFilter handles label=None, codim=0. + # Take an explicit label including all entities otherwise. label_name = "depth" subdomain_id = subdim elif label_name is None: From 5cb9b41c89f50c4949b35452040b2111d462dcf1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 18 Aug 2026 14:15:15 +0100 Subject: [PATCH 03/11] review comments --- firedrake/cython/dmcommon.pyx | 73 +++++++++++++++++------------------ firedrake/mesh.py | 5 --- 2 files changed, 36 insertions(+), 42 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 7005230972..35cae685e1 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4019,41 +4019,41 @@ def submesh_create(PETSc.DM dm, const PetscInt *stratum_indices = NULL if label_name is None: - # label=None covers all cells. DMPlexFilter already handles this. - temp_label = None - else: - # Cast subdomain_id into an iterable - if isinstance(subdomain_id, str) or not isinstance(subdomain_id, Sequence): - subdomain_id = (subdomain_id,) - # Take the union of the all the label values - label = dm.getLabel(label_name) - points = PETSc.IS() - for sub in subdomain_id: - if isinstance(sub, Integral): - subpoints = label.getStratumIS(sub) - elif sub == "on_boundary": - subpoints = dm.getStratumIS("exterior_facets", 1) - else: - raise ValueError(f"Submesh construction got invalid subdomain_id {sub}.") - if points: - points = points.union(subpoints) - else: - points = subpoints - # Create temp_label that contains no lower-dimensional points. - dm.createLabel(temp_label_name) - temp_label = dm.getLabel(temp_label_name) + # Default to all entities of the given dimension. + label_name = "depth" + subdomain_id = subdim + # Cast subdomain_id into an iterable + if isinstance(subdomain_id, str) or not isinstance(subdomain_id, Sequence): + subdomain_id = (subdomain_id,) + # Take the union of the all the label values + label = dm.getLabel(label_name) + points = PETSc.IS() + for sub in subdomain_id: + if isinstance(sub, Integral): + subpoints = label.getStratumIS(sub) + elif sub == "on_boundary": + subpoints = dm.getStratumIS("exterior_facets", 1) + else: + raise ValueError(f"Submesh construction got invalid subdomain_id {sub}.") if points: - CHKERR(ISGetSize(points.iset, &stratum_size)) - if stratum_size > 0: - CHKERR(ISGetIndices(points.iset, &stratum_indices)) - CHKERR(DMPlexGetDepthStratum(dm.dm, subdim, &pStart, &pEnd)) - for i in range(stratum_size): - p = stratum_indices[i] - # Only include points on the submesh topological dimension, - # culling all lower-dimensional points. - if pStart <= p < pEnd: - CHKERR(DMLabelSetValue(temp_label.dmlabel, p, label_value)) - CHKERR(ISRestoreIndices(points.iset, &stratum_indices)) + points = points.union(subpoints) + else: + points = subpoints + # Create temp_label that contains no lower-dimensional points. + dm.createLabel(temp_label_name) + temp_label = dm.getLabel(temp_label_name) + if points: + CHKERR(ISGetSize(points.iset, &stratum_size)) + if stratum_size > 0: + CHKERR(ISGetIndices(points.iset, &stratum_indices)) + CHKERR(DMPlexGetDepthStratum(dm.dm, subdim, &pStart, &pEnd)) + for i in range(stratum_size): + p = stratum_indices[i] + # Only include points on the submesh topological dimension, + # culling all lower-dimensional points. + if pStart <= p < pEnd: + CHKERR(DMLabelSetValue(temp_label.dmlabel, p, label_value)) + CHKERR(ISRestoreIndices(points.iset, &stratum_indices)) # Make submesh using temp_label. subdm, ownership_transfer_sf = dm.filter(label=temp_label, value=label_value, @@ -4061,9 +4061,8 @@ def submesh_create(PETSc.DM dm, sanitizeSubMesh=PETSC_TRUE, comm=comm) # Destroy temp_label. - if temp_label is not None: - dm.removeLabel(temp_label_name) - subdm.removeLabel(temp_label_name) + dm.removeLabel(temp_label_name) + subdm.removeLabel(temp_label_name) submesh_update_facet_labels(dm, subdm) submesh_correct_entity_classes(dm, subdm, ownership_transfer_sf) return subdm diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 5631aa5473..5c3a0b8476 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -5023,11 +5023,6 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig if subdomain_id is None: if label_name is not None: raise ValueError("subdomain_id=None requires label_name=None.") - if subdim != dim: - # DMPlexFilter handles label=None, codim=0. - # Take an explicit label including all entities otherwise. - label_name = "depth" - subdomain_id = subdim elif label_name is None: if subdim == dim: label_name = dmcommon.CELL_SETS_LABEL From 6d3da37889634c3f41ec3455ad8af8cad412f8bb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 19 Aug 2026 09:23:11 +0100 Subject: [PATCH 04/11] restore original API --- firedrake/cython/dmcommon.pyx | 12 ++++-------- firedrake/mesh.py | 3 +++ 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 35cae685e1..b0f1e36b9f 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4001,10 +4001,10 @@ def submesh_create(PETSc.DM dm, DMPlex representing the mesh topology subdim : int Topological dimension of the submesh - label_name : str | None - Name of the label, or `None` to select every cell - subdomain_id : int | Sequence | None - Values in the label, unused if ``label_name`` is `None` + label_name : str + Name of the label + subdomain_id : int | Sequence + Values in the label ignore_label_halo : bool If labeled points in the halo are ignored. comm : PETSc.Comm | None @@ -4018,10 +4018,6 @@ def submesh_create(PETSc.DM dm, PetscInt pStart, pEnd, p, i, stratum_size = 0, label_value = 1 const PetscInt *stratum_indices = NULL - if label_name is None: - # Default to all entities of the given dimension. - label_name = "depth" - subdomain_id = subdim # Cast subdomain_id into an iterable if isinstance(subdomain_id, str) or not isinstance(subdomain_id, Sequence): subdomain_id = (subdomain_id,) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 5c3a0b8476..f30f0d4ec0 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -5023,6 +5023,9 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig if subdomain_id is None: if label_name is not None: raise ValueError("subdomain_id=None requires label_name=None.") + # Select all entities + label_name = "depth" + subdomain_id = subdim elif label_name is None: if subdim == dim: label_name = dmcommon.CELL_SETS_LABEL From db87f54fe2253d0973dec44a63758489b09cab63 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 19 Aug 2026 09:34:24 +0100 Subject: [PATCH 05/11] Move all default handling into cython wrapper --- firedrake/cython/dmcommon.pyx | 39 ++++++++++++++++++++++++++++++----- firedrake/mesh.py | 28 ------------------------- 2 files changed, 34 insertions(+), 33 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index b0f1e36b9f..215ddffb42 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -3999,12 +3999,12 @@ def submesh_create(PETSc.DM dm, ---------- dm : PETSc.DM DMPlex representing the mesh topology - subdim : int + subdim : int | None Topological dimension of the submesh - label_name : str - Name of the label - subdomain_id : int | Sequence - Values in the label + label_name : str | None + Name of the label, or `None` to select every cell + subdomain_id : int | Sequence | None + Values in the label, unused if ``label_name`` is `None` ignore_label_halo : bool If labeled points in the halo are ignored. comm : PETSc.Comm | None @@ -4018,6 +4018,35 @@ def submesh_create(PETSc.DM dm, PetscInt pStart, pEnd, p, i, stratum_size = 0, label_value = 1 const PetscInt *stratum_indices = NULL + # Parse default subdim, label_name, and subdomain_id + dim = dm.getDimension() + if subdomain_id == "on_boundary": + if subdim is None: + subdim = dim - 1 + elif subdim != dim - 1: + raise ValueError('subdomain_id="on_boundary" requires subdim=dim-1') + if label_name is None: + label_name = "exterior_facets" + elif label_name != "exterior_facets": + raise ValueError('subdomain_id="on_boundary" requires label_name="exterior_facets"') + subdomain_id = 1 + + if subdim is None: + subdim = dim + if subdim not in {dim, dim - 1}: + raise NotImplementedError(f"Found submesh dim ({subdim}) and parent dim ({dim})") + if subdomain_id is None: + if label_name is not None: + raise ValueError("subdomain_id=None requires label_name=None.") + # Select all entities + label_name = "depth" + subdomain_id = subdim + elif label_name is None: + if subdim == dim: + label_name = CELL_SETS_LABEL + elif subdim == dim - 1: + label_name = FACE_SETS_LABEL + # Cast subdomain_id into an iterable if isinstance(subdomain_id, str) or not isinstance(subdomain_id, Sequence): subdomain_id = (subdomain_id,) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index f30f0d4ec0..d5f2428237 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -5003,34 +5003,6 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig elif isinstance(mesh.topology, VertexOnlyMeshTopology): raise NotImplementedError("Can not create a submesh of a ``VertexOnlyMesh``") - if subdomain_id == "on_boundary": - if subdim is None: - subdim = mesh.topological_dimension - 1 - elif subdim != mesh.topological_dimension - 1: - raise ValueError('subdomain_id="on_boundary" requires subdim=dim-1') - if label_name is None: - label_name = "exterior_facets" - elif label_name != "exterior_facets": - raise ValueError('subdomain_id="on_boundary" requires label_name="exterior_facets"') - subdomain_id = 1 - - if subdim is None: - subdim = mesh.topological_dimension - plex = mesh.topology_dm - dim = plex.getDimension() - if subdim not in {dim, dim - 1}: - raise NotImplementedError(f"Found submesh dim ({subdim}) and parent dim ({dim})") - if subdomain_id is None: - if label_name is not None: - raise ValueError("subdomain_id=None requires label_name=None.") - # Select all entities - label_name = "depth" - subdomain_id = subdim - elif label_name is None: - if subdim == dim: - label_name = dmcommon.CELL_SETS_LABEL - elif subdim == dim - 1: - label_name = dmcommon.FACE_SETS_LABEL subplex = dmcommon.submesh_create(plex, subdim, label_name, subdomain_id, ignore_halo, comm=comm) comm = comm or mesh.comm From 26f9cec647c9c52a67c8e2c32ce0992c9ad3d65d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 19 Aug 2026 10:38:51 +0100 Subject: [PATCH 06/11] Apply suggestions from code review Co-authored-by: Pablo Brubeck --- firedrake/cython/dmcommon.pyx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 215ddffb42..715be0303f 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4000,7 +4000,8 @@ def submesh_create(PETSc.DM dm, dm : PETSc.DM DMPlex representing the mesh topology subdim : int | None - Topological dimension of the submesh + Topological dimension of the submesh, or None to be infered from other kwargs. + See :func:`~mesh.Submesh`. label_name : str | None Name of the label, or `None` to select every cell subdomain_id : int | Sequence | None @@ -4034,7 +4035,8 @@ def submesh_create(PETSc.DM dm, if subdim is None: subdim = dim if subdim not in {dim, dim - 1}: - raise NotImplementedError(f"Found submesh dim ({subdim}) and parent dim ({dim})") + raise NotImplementedError(f"Submesh construction is only implemented for codimension 0 or 1. + Found submesh dim ({subdim}) and parent dim ({dim})") if subdomain_id is None: if label_name is not None: raise ValueError("subdomain_id=None requires label_name=None.") From 6a79c91024a462729fbb8e6182d793c19f75fac9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 19 Aug 2026 10:39:46 +0100 Subject: [PATCH 07/11] Apply suggestion from @pbrubeck --- firedrake/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index d5f2428237..fcc1923206 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -5003,7 +5003,7 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig elif isinstance(mesh.topology, VertexOnlyMeshTopology): raise NotImplementedError("Can not create a submesh of a ``VertexOnlyMesh``") - subplex = dmcommon.submesh_create(plex, subdim, label_name, subdomain_id, ignore_halo, comm=comm) + subplex = dmcommon.submesh_create(mesh.topology_dm, subdim, label_name, subdomain_id, ignore_halo, comm=comm) comm = comm or mesh.comm name = name or _generate_default_submesh_name(mesh.name) From 68ba4c1c40ad3493957b16df1b04755463a60f69 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 19 Aug 2026 10:52:35 +0100 Subject: [PATCH 08/11] Fixes --- firedrake/cython/dmcommon.pyx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 715be0303f..12dab4e46d 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4000,8 +4000,8 @@ def submesh_create(PETSc.DM dm, dm : PETSc.DM DMPlex representing the mesh topology subdim : int | None - Topological dimension of the submesh, or None to be infered from other kwargs. - See :func:`~mesh.Submesh`. + Topological dimension of the submesh, or None to be inferred from other kwargs. + See :func:`~.mesh.Submesh`. label_name : str | None Name of the label, or `None` to select every cell subdomain_id : int | Sequence | None @@ -4035,8 +4035,8 @@ def submesh_create(PETSc.DM dm, if subdim is None: subdim = dim if subdim not in {dim, dim - 1}: - raise NotImplementedError(f"Submesh construction is only implemented for codimension 0 or 1. - Found submesh dim ({subdim}) and parent dim ({dim})") + raise NotImplementedError(f"Submesh construction is only implemented for codimension 0 or 1. " + "Found submesh dim ({subdim}) and parent dim ({dim})") if subdomain_id is None: if label_name is not None: raise ValueError("subdomain_id=None requires label_name=None.") From 22afea54e9b25ec3b1a2758b06f10a85c899cff6 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 19 Aug 2026 14:10:59 +0100 Subject: [PATCH 09/11] Apply suggestion from @pbrubeck --- firedrake/cython/dmcommon.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 12dab4e46d..4935fcc750 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -3988,7 +3988,7 @@ def create_halo_exchange_sf(PETSc.DM dm): @cython.boundscheck(False) @cython.wraparound(False) def submesh_create(PETSc.DM dm, - PetscInt subdim, + subdim, label_name, subdomain_id, PetscBool ignore_label_halo, From 4b5c791ba1d82efe0e3b2895c53ef4247caf8f53 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 19 Aug 2026 16:06:35 +0100 Subject: [PATCH 10/11] Apply suggestions from code review Co-authored-by: Pablo Brubeck --- firedrake/cython/dmcommon.pyx | 2 ++ firedrake/mesh.py | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 4935fcc750..34ac6907bc 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4087,6 +4087,8 @@ def submesh_create(PETSc.DM dm, ignoreHalo=ignore_label_halo, sanitizeSubMesh=PETSC_TRUE, comm=comm) + if subdm.getDimension() != subdim: + raise RuntimeError(f"Found subplex dim ({subdm.getDimension()}) != expected ({subdim})") # Destroy temp_label. dm.removeLabel(temp_label_name) subdm.removeLabel(temp_label_name) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index fcc1923206..d52423e02d 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -5009,7 +5009,6 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig name = name or _generate_default_submesh_name(mesh.name) subplex.setName(_generate_default_mesh_topology_name(name)) if subplex.getDimension() != subdim: - raise RuntimeError(f"Found subplex dim ({subplex.getDimension()}) != expected ({subdim})") if reorder is None: # Ideally we should set perm_is = mesh._dm_renumbering[label_indices] reorder = mesh._did_reordering From 0f8296675b142ac216e40e4a8ab30eb4c70bdef3 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 19 Aug 2026 16:07:17 +0100 Subject: [PATCH 11/11] Apply suggestion from @pbrubeck --- firedrake/mesh.py | 1 - 1 file changed, 1 deletion(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index d52423e02d..99bde8fe43 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -5008,7 +5008,6 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig comm = comm or mesh.comm name = name or _generate_default_submesh_name(mesh.name) subplex.setName(_generate_default_mesh_topology_name(name)) - if subplex.getDimension() != subdim: if reorder is None: # Ideally we should set perm_is = mesh._dm_renumbering[label_indices] reorder = mesh._did_reordering