From afde9fafb7899f78cc20f4065ef6c0bc706d9b69 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Mon, 9 Mar 2026 15:42:36 +0100 Subject: [PATCH 01/15] fix submesh labels --- firedrake/cython/dmcommon.pyx | 70 ++++++++++++++++++++++++++++++++--- 1 file changed, 64 insertions(+), 6 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index ec7ab631af..ef5fb1782c 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -3951,11 +3951,18 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): Notes ----- - This function marks the new exterior facets with current max label value + 1 in "Face Sets". + For same-dimension submeshes (subdim == dim), marks new exterior facets + with current max label value + 1 in "Face Sets". + + For codimension-1 submeshes (subdim == dim - 1), propagates the parent's + lower-dimensional labels ("Edge Sets" for 3D parents, "Vertex Sets" for 2D + parents) to "Face Sets" on the subdm. Exterior facets without a parent + label are marked with a fresh label value. """ cdef: PetscInt dim, subdim, pStart, pEnd, f, subfStart, subfEnd, subf, sub_ext_facet_size, next_label_val, i + PetscInt label_val PETSc.IS subpoint_is PETSc.IS sub_ext_facet_is const PetscInt *subpoint_indices = NULL @@ -3963,24 +3970,26 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): char *int_facet_label_name = "interior_facets" char *ext_facet_label_name = "exterior_facets" char *face_sets_label_name = "Face Sets" + char *edge_sets_label_name = "Edge Sets" + char *vertex_sets_label_name = "Vertex Sets" + char *parent_label_cname = NULL DMLabel ext_facet_label + DMLabel parent_facet_lbl PETSc.DMLabel sub_int_facet_label, sub_ext_facet_label PetscBool has_point dim = dm.getDimension() subdim = subdm.getDimension() - if subdim != dim: - # What labels the submesh should have by default is not trivial. - # Think harder and do something useful here if necessary. + if subdim != dim and subdim != dim - 1: return # Mark interior and exterior facets label_facets(subdm) sub_int_facet_label = subdm.getLabel("interior_facets") sub_ext_facet_label = subdm.getLabel("exterior_facets") - # Mark new exterior facets with current max label value + 1 in "Face Sets" subpoint_is = subdm.getSubpointIS() CHKERR(ISGetIndices(subpoint_is.iset, &subpoint_indices)) if subdim == dim: + # Mark new exterior facets with current max label value + 1 in "Face Sets" with dm.getLabelIdIS(FACE_SETS_LABEL) as label_value_indices: next_label_val = label_value_indices.max() + 1 if len(label_value_indices) > 0 else 0 next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) @@ -4006,7 +4015,56 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): if sub_ext_facet_is.iset: CHKERR(ISRestoreIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) else: - raise NotImplementedError("Currently, only implemented for cell submesh") + # subdim == dim - 1: codimension-1 submesh. + # Submesh facets are at depth (subdim - 1) in the parent mesh. + # Propagate parent labels at that depth into "Face Sets" on subdm. + if subdim - 1 == 1: + parent_label_cname = edge_sets_label_name + parent_label_pyname = "Edge Sets" + elif subdim - 1 == 0: + parent_label_cname = vertex_sets_label_name + parent_label_pyname = "Vertex Sets" + else: + parent_label_cname = NULL + parent_label_pyname = None + has_parent_facet_label = (parent_label_pyname is not None + and dm.hasLabel(parent_label_pyname)) + next_label_val = 0 + if has_parent_facet_label: + with dm.getLabelIdIS(parent_label_pyname) as label_value_indices: + next_label_val = label_value_indices.max() + 1 if len(label_value_indices) > 0 else 0 + next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) + # DMPlexFilter propagates parent "Face Sets" to the subdm, but those + # values sit on cells (codim-0), not facets. Remove them so they don't + # contaminate the exterior facet markers we are about to create. + if subdm.hasLabel(FACE_SETS_LABEL): + subdm.removeLabel(FACE_SETS_LABEL) + subdm.createLabel(FACE_SETS_LABEL) + sub_ext_facet_size = subdm.getStratumSize("exterior_facets", 1) + sub_ext_facet_is = subdm.getStratumIS("exterior_facets", 1) + if sub_ext_facet_is.iset: + CHKERR(ISGetIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) + subfStart, subfEnd = subdm.getHeightStratum(1) + if has_parent_facet_label: + CHKERR(DMGetLabel(dm.dm, parent_label_cname, &parent_facet_lbl)) + pStart, pEnd = dm.getChart() + CHKERR(DMLabelCreateIndex(parent_facet_lbl, pStart, pEnd)) + for i in range(sub_ext_facet_size): + subf = sub_ext_facet_indices[i] + if subf < subfStart or subf >= subfEnd: + continue + f = subpoint_indices[subf] + label_val = -1 + if has_parent_facet_label: + CHKERR(DMLabelGetValue(parent_facet_lbl, f, &label_val)) + if label_val >= 0: + CHKERR(DMSetLabelValue(subdm.dm, face_sets_label_name, subf, label_val)) + else: + CHKERR(DMSetLabelValue(subdm.dm, face_sets_label_name, subf, next_label_val)) + if has_parent_facet_label: + CHKERR(DMLabelDestroyIndex(parent_facet_lbl)) + if sub_ext_facet_is.iset: + CHKERR(ISRestoreIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) CHKERR(ISRestoreIndices(subpoint_is.iset, &subpoint_indices)) subdm.removeLabel("interior_facets") subdm.removeLabel("exterior_facets") From 19d97e8dbf722bd345e4169aa062f84ddd0d23d9 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Mon, 9 Mar 2026 15:55:33 +0100 Subject: [PATCH 02/15] add tests --- .../submesh/test_submesh_codim1_labels.py | 219 ++++++++++++++++++ 1 file changed, 219 insertions(+) create mode 100644 tests/firedrake/submesh/test_submesh_codim1_labels.py diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py new file mode 100644 index 0000000000..d9a9cc4794 --- /dev/null +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -0,0 +1,219 @@ +import pytest +import numpy as np +from firedrake import * + + +def _get_dm_vertex_coords(dm): + """Return ``{vertex_point: np.array([x, y, ...])}`` from the DM coordinate section.""" + coord_sec = dm.getCoordinateSection() + coord_arr = dm.getCoordinatesLocal().array_r + vStart, vEnd = dm.getDepthStratum(0) + result = {} + for v in range(vStart, vEnd): + off = coord_sec.getOffset(v) + dof = coord_sec.getDof(v) + result[v] = coord_arr[off:off + dof].copy() + return result + + +def _entity_vertex_coords(dm, point, vc): + """Return vertex coordinates for all vertices in the closure of *point*.""" + vStart, vEnd = dm.getDepthStratum(0) + closure = dm.getTransitiveClosure(point)[0] + verts = [c for c in closure if vStart <= c < vEnd] + if not verts: + return np.empty((0, dm.getCoordinateDim())) + return np.array([vc[v] for v in verts]) + + +def _all_at(cs, axis, value): + return np.allclose(cs[:, axis], value) + + +# --------------------------------------------------------------------------- +# 3D -> 2D helpers +# --------------------------------------------------------------------------- + +FACE_TAG_3D = 100 +EDGE_A = 1 +EDGE_B = 2 +EDGE_C = 3 + + +def _make_cube_with_edge_sets(extra_edge_c=False): + """UnitCubeMesh(2,2,2) with Face Sets on z=0 boundary faces and Edge Sets + on selected boundary edges of the z=0 face. + + If *extra_edge_c* is True an additional label EDGE_C is placed on edges + at z=1, x=0 (not part of the z=0 submesh). + """ + mesh = UnitCubeMesh(2, 2, 2) + dm = mesh.topology_dm + vc = _get_dm_vertex_coords(dm) + eStart, eEnd = dm.getDepthStratum(1) + fStart, fEnd = dm.getDepthStratum(2) + + for f in range(fStart, fEnd): + cs = _entity_vertex_coords(dm, f, vc) + if cs.shape[0] >= 3 and _all_at(cs, 2, 0.0) and dm.getSupportSize(f) == 1: + dm.setLabelValue("Face Sets", f, FACE_TAG_3D) + + dm.createLabel("Edge Sets") + for e in range(eStart, eEnd): + cs = _entity_vertex_coords(dm, e, vc) + if cs.shape[0] < 2: + continue + if _all_at(cs, 2, 0.0) and _all_at(cs, 0, 0.0): + dm.setLabelValue("Edge Sets", e, EDGE_A) + elif _all_at(cs, 2, 0.0) and _all_at(cs, 0, 1.0): + dm.setLabelValue("Edge Sets", e, EDGE_B) + elif extra_edge_c and _all_at(cs, 2, 1.0) and _all_at(cs, 0, 0.0): + dm.setLabelValue("Edge Sets", e, EDGE_C) + + return mesh + + +# --------------------------------------------------------------------------- +# 2D -> 1D helpers +# --------------------------------------------------------------------------- + +FACE_TAG_2D = 200 +VTX_A = 1 +VTX_B = 2 + + +def _make_square_with_vertex_sets(): + """UnitSquareMesh(2,2) with Face Sets on y=0 boundary edges and Vertex + Sets on the corner vertices of that edge. + """ + mesh = UnitSquareMesh(2, 2) + dm = mesh.topology_dm + vc = _get_dm_vertex_coords(dm) + vStart, vEnd = dm.getDepthStratum(0) + eStart, eEnd = dm.getDepthStratum(1) + + for e in range(eStart, eEnd): + cs = _entity_vertex_coords(dm, e, vc) + if cs.shape[0] >= 2 and _all_at(cs, 1, 0.0) and dm.getSupportSize(e) == 1: + dm.setLabelValue("Face Sets", e, FACE_TAG_2D) + + dm.createLabel("Vertex Sets") + for v in range(vStart, vEnd): + coords = vc[v] + if np.isclose(coords[1], 0.0) and np.isclose(coords[0], 0.0): + dm.setLabelValue("Vertex Sets", v, VTX_A) + elif np.isclose(coords[1], 0.0) and np.isclose(coords[0], 1.0): + dm.setLabelValue("Vertex Sets", v, VTX_B) + + return mesh + + +# --------------------------------------------------------------------------- +# Tests – 3D parent -> 2D submesh (Edge Sets → Face Sets) +# --------------------------------------------------------------------------- + +def test_submesh_codim1_edge_sets_propagated(): + """Parent 'Edge Sets' tags appear in submesh exterior_facets.unique_markers.""" + mesh = _make_cube_with_edge_sets() + submesh = Submesh(mesh, 2, FACE_TAG_3D) + markers = submesh.exterior_facets.unique_markers + assert markers is not None + marker_set = set(markers) + assert EDGE_A in marker_set, f"EDGE_A={EDGE_A} not in {marker_set}" + assert EDGE_B in marker_set, f"EDGE_B={EDGE_B} not in {marker_set}" + + +def test_submesh_codim1_edge_sets_excludes_nonsubmesh(): + """Edge Sets on parent edges outside the submesh are excluded.""" + mesh = _make_cube_with_edge_sets(extra_edge_c=True) + submesh = Submesh(mesh, 2, FACE_TAG_3D) + markers = submesh.exterior_facets.unique_markers + assert markers is not None + marker_set = set(markers) + assert EDGE_A in marker_set + assert EDGE_B in marker_set + assert EDGE_C not in marker_set, f"EDGE_C={EDGE_C} should not be in {marker_set}" + + +def test_submesh_codim1_unlabeled_get_default_value(): + """Exterior facets without a parent Edge Sets label get max(label_vals)+1.""" + mesh = _make_cube_with_edge_sets() + submesh = Submesh(mesh, 2, FACE_TAG_3D) + markers = submesh.exterior_facets.unique_markers + assert markers is not None + marker_set = set(markers) + next_val = max(EDGE_A, EDGE_B) + 1 + assert next_val in marker_set, f"Default {next_val} not in {marker_set}" + + +def test_submesh_codim1_ds_edge_sets(): + """ds(tag) with propagated Edge Sets integrates the correct boundary length.""" + mesh = _make_cube_with_edge_sets() + submesh = Submesh(mesh, 2, FACE_TAG_3D) + # EDGE_A labels edges at x=0, z=0: total length = 1 + assert abs(assemble(Constant(1.) * ds(EDGE_A, domain=submesh)) - 1.0) < 1e-12 + # EDGE_B labels edges at x=1, z=0: total length = 1 + assert abs(assemble(Constant(1.) * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 + + +# --------------------------------------------------------------------------- +# Tests – 2D parent -> 1D submesh (Vertex Sets → Face Sets) +# --------------------------------------------------------------------------- + +def test_submesh_codim1_vertex_sets_propagated(): + """Parent 'Vertex Sets' tags appear in 1D submesh exterior_facets.unique_markers.""" + mesh = _make_square_with_vertex_sets() + submesh = Submesh(mesh, 1, FACE_TAG_2D) + markers = submesh.exterior_facets.unique_markers + assert markers is not None + marker_set = set(markers) + assert VTX_A in marker_set, f"VTX_A={VTX_A} not in {marker_set}" + assert VTX_B in marker_set, f"VTX_B={VTX_B} not in {marker_set}" + + +def test_submesh_codim1_ds_vertex_sets(): + """ds(tag) with propagated Vertex Sets works on a 1D submesh.""" + mesh = _make_square_with_vertex_sets() + submesh = Submesh(mesh, 1, FACE_TAG_2D) + # Each tagged vertex is a single point (0-D measure = 1) + assert abs(assemble(Constant(1.) * ds(VTX_A, domain=submesh)) - 1.0) < 1e-12 + assert abs(assemble(Constant(1.) * ds(VTX_B, domain=submesh)) - 1.0) < 1e-12 + + +# --------------------------------------------------------------------------- +# Tests – no parent label +# --------------------------------------------------------------------------- + +def test_submesh_codim1_no_parent_edge_sets(): + """Codim-1 submesh with no Edge Sets on the parent still labels exterior facets.""" + mesh = UnitCubeMesh(2, 2, 2) + dm = mesh.topology_dm + vc = _get_dm_vertex_coords(dm) + fStart, fEnd = dm.getDepthStratum(2) + + for f in range(fStart, fEnd): + cs = _entity_vertex_coords(dm, f, vc) + if cs.shape[0] >= 3 and _all_at(cs, 2, 0.0) and dm.getSupportSize(f) == 1: + dm.setLabelValue("Face Sets", f, FACE_TAG_3D) + + submesh = Submesh(mesh, 2, FACE_TAG_3D) + markers = submesh.exterior_facets.unique_markers + assert markers is not None + # All exterior facets should receive the default label (0) + assert 0 in set(markers) + + +# --------------------------------------------------------------------------- +# Parallel +# --------------------------------------------------------------------------- + +@pytest.mark.parallel(nprocs=2) +def test_submesh_codim1_edge_sets_propagated_parallel(): + """Edge Sets propagation works in parallel (2 ranks).""" + mesh = _make_cube_with_edge_sets() + submesh = Submesh(mesh, 2, FACE_TAG_3D) + markers = submesh.exterior_facets.unique_markers + assert markers is not None + marker_set = set(markers) + assert EDGE_A in marker_set, f"EDGE_A={EDGE_A} not in {marker_set}" + assert EDGE_B in marker_set, f"EDGE_B={EDGE_B} not in {marker_set}" From 1ac8e61cd169d369f77af6c9267f3165b0830115 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 08:49:51 +0100 Subject: [PATCH 03/15] make an else condition explicit --- firedrake/cython/dmcommon.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index ef5fb1782c..b9826b1578 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4014,8 +4014,8 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): CHKERR(DMLabelDestroyIndex(ext_facet_label)) if sub_ext_facet_is.iset: CHKERR(ISRestoreIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) - else: - # subdim == dim - 1: codimension-1 submesh. + elif subdim == dim - 1: + # Codimension-1 submesh. # Submesh facets are at depth (subdim - 1) in the parent mesh. # Propagate parent labels at that depth into "Face Sets" on subdm. if subdim - 1 == 1: From 19b188f90fbe624d736d222b39edd71e575ad45c Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 09:19:53 +0100 Subject: [PATCH 04/15] add some tests --- firedrake/cython/dmcommon.pyx | 225 ++++++++++-------- .../submesh/test_submesh_codim1_labels.py | 138 ++++++++++- 2 files changed, 250 insertions(+), 113 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index b9826b1578..e38e8767fb 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -3937,134 +3937,155 @@ def submesh_correct_entity_classes(PETSc.DM dm, CHKERR(DMLabelDestroyIndex(lbl_ghost)) +cdef void _label_new_exterior_facets( + PETSc.DM dm, PETSc.DM subdm, + const PetscInt *subpoint_indices, + const PetscInt *sub_ext_facet_indices, + PetscInt sub_ext_facet_size, + PetscInt subfStart, PetscInt subfEnd, +): + """Same-dimension helper: tag exterior facets that were interior in the parent. + + Facets of the submesh that are *not* in the parent's ``exterior_facets`` + label are new boundary and receive ``max("Face Sets") + 1``. + """ + cdef: + PetscInt pStart, pEnd, next_label_val, subf, f, i + DMLabel parent_ext_label + PetscBool has_point + + with dm.getLabelIdIS(FACE_SETS_LABEL) as ids: + next_label_val = ids.max() + 1 if len(ids) > 0 else 0 + next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) + subdm.createLabel(FACE_SETS_LABEL) + + CHKERR(DMGetLabel(dm.dm, b"exterior_facets", &parent_ext_label)) + pStart, pEnd = dm.getChart() + CHKERR(DMLabelCreateIndex(parent_ext_label, pStart, pEnd)) + for i in range(sub_ext_facet_size): + subf = sub_ext_facet_indices[i] + if subfStart <= subf < subfEnd: + f = subpoint_indices[subf] + CHKERR(DMLabelHasPoint(parent_ext_label, f, &has_point)) + if not has_point: + CHKERR(DMSetLabelValue(subdm.dm, b"Face Sets", subf, next_label_val)) + CHKERR(DMLabelDestroyIndex(parent_ext_label)) + + +cdef void _propagate_parent_facet_labels( + PETSc.DM dm, PETSc.DM subdm, + PetscInt subdim, + const PetscInt *subpoint_indices, + const PetscInt *sub_ext_facet_indices, + PetscInt sub_ext_facet_size, + PetscInt subfStart, PetscInt subfEnd, +): + """Codimension-1 helper: map the parent's lower-dimensional label into "Face Sets". + + For a 3D→2D submesh the parent's "Edge Sets" are propagated; for 2D→1D + the parent's "Vertex Sets" are used. Exterior facets without a + corresponding parent label value receive a default tag. + """ + cdef: + PetscInt pStart, pEnd, next_label_val, label_val, subf, f, i + DMLabel parent_label + + if subdim == 2: + parent_label_name = b"Edge Sets" + elif subdim == 1: + parent_label_name = b"Vertex Sets" + else: + parent_label_name = None + + has_parent_label = (parent_label_name is not None + and dm.hasLabel(parent_label_name)) + + next_label_val = 0 + if has_parent_label: + with dm.getLabelIdIS(parent_label_name) as ids: + next_label_val = ids.max() + 1 if len(ids) > 0 else 0 + next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) + + # DMPlexFilter may propagate cell-level "Face Sets" to the subdm; + # remove them so they don't contaminate the facet markers. + if subdm.hasLabel(FACE_SETS_LABEL): + subdm.removeLabel(FACE_SETS_LABEL) + subdm.createLabel(FACE_SETS_LABEL) + + if has_parent_label: + CHKERR(DMGetLabel(dm.dm, parent_label_name, &parent_label)) + pStart, pEnd = dm.getChart() + CHKERR(DMLabelCreateIndex(parent_label, pStart, pEnd)) + for i in range(sub_ext_facet_size): + subf = sub_ext_facet_indices[i] + if subfStart <= subf < subfEnd: + f = subpoint_indices[subf] + label_val = -1 + if has_parent_label: + CHKERR(DMLabelGetValue(parent_label, f, &label_val)) + if label_val >= 0: + CHKERR(DMSetLabelValue(subdm.dm, b"Face Sets", subf, label_val)) + else: + CHKERR(DMSetLabelValue(subdm.dm, b"Face Sets", subf, next_label_val)) + if has_parent_label: + CHKERR(DMLabelDestroyIndex(parent_label)) + + @cython.boundscheck(False) @cython.wraparound(False) def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): - """Update facet labels of subdm taking the new exterior facet points into account. + """Update "Face Sets" on *subdm* for its exterior facets. Parameters ---------- dm : PETSc.DM - The parent dm. + The parent DM. subdm : PETSc.DM - The subdm. + The sub-DM whose facet labels are updated. Notes ----- - For same-dimension submeshes (subdim == dim), marks new exterior facets - with current max label value + 1 in "Face Sets". - - For codimension-1 submeshes (subdim == dim - 1), propagates the parent's - lower-dimensional labels ("Edge Sets" for 3D parents, "Vertex Sets" for 2D - parents) to "Face Sets" on the subdm. Exterior facets without a parent - label are marked with a fresh label value. + * **Same-dimension** (``subdim == dim``): new exterior facets (those that + were interior in the parent) are tagged with ``max("Face Sets") + 1``. + * **Codimension-1** (``subdim == dim - 1``): the parent's lower-dimensional + labels ("Edge Sets" for 3D, "Vertex Sets" for 2D) are mapped into + "Face Sets" on the subdm. Unlabeled facets get a default value. """ cdef: - PetscInt dim, subdim, pStart, pEnd, f, subfStart, subfEnd, subf, sub_ext_facet_size, next_label_val, i - PetscInt label_val - PETSc.IS subpoint_is - PETSc.IS sub_ext_facet_is + PetscInt dim, subdim, subfStart, subfEnd, sub_ext_facet_size + PETSc.IS subpoint_is, sub_ext_facet_is const PetscInt *subpoint_indices = NULL const PetscInt *sub_ext_facet_indices = NULL - char *int_facet_label_name = "interior_facets" - char *ext_facet_label_name = "exterior_facets" - char *face_sets_label_name = "Face Sets" - char *edge_sets_label_name = "Edge Sets" - char *vertex_sets_label_name = "Vertex Sets" - char *parent_label_cname = NULL - DMLabel ext_facet_label - DMLabel parent_facet_lbl - PETSc.DMLabel sub_int_facet_label, sub_ext_facet_label - PetscBool has_point dim = dm.getDimension() subdim = subdm.getDimension() if subdim != dim and subdim != dim - 1: return - # Mark interior and exterior facets + label_facets(subdm) - sub_int_facet_label = subdm.getLabel("interior_facets") - sub_ext_facet_label = subdm.getLabel("exterior_facets") + subpoint_is = subdm.getSubpointIS() CHKERR(ISGetIndices(subpoint_is.iset, &subpoint_indices)) + sub_ext_facet_size = subdm.getStratumSize("exterior_facets", 1) + sub_ext_facet_is = subdm.getStratumIS("exterior_facets", 1) + if sub_ext_facet_is.iset: + CHKERR(ISGetIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) + subfStart, subfEnd = subdm.getHeightStratum(1) + if subdim == dim: - # Mark new exterior facets with current max label value + 1 in "Face Sets" - with dm.getLabelIdIS(FACE_SETS_LABEL) as label_value_indices: - next_label_val = label_value_indices.max() + 1 if len(label_value_indices) > 0 else 0 - next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) - subdm.createLabel(FACE_SETS_LABEL) - sub_ext_facet_size = subdm.getStratumSize("exterior_facets", 1) - sub_ext_facet_is = subdm.getStratumIS("exterior_facets", 1) - if sub_ext_facet_is.iset: - CHKERR(ISGetIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) - CHKERR(DMGetLabel(dm.dm, ext_facet_label_name, &ext_facet_label)) - pStart, pEnd = dm.getChart() - CHKERR(DMLabelCreateIndex(ext_facet_label, pStart, pEnd)) - subfStart, subfEnd = subdm.getHeightStratum(1) - for i in range(sub_ext_facet_size): - subf = sub_ext_facet_indices[i] - if subf < subfStart or subf >= subfEnd: - continue - f = subpoint_indices[subf] - CHKERR(DMLabelHasPoint(ext_facet_label, f, &has_point)) - if not has_point: - # Found a new exterior facet - CHKERR(DMSetLabelValue(subdm.dm, face_sets_label_name, subf, next_label_val)) - CHKERR(DMLabelDestroyIndex(ext_facet_label)) - if sub_ext_facet_is.iset: - CHKERR(ISRestoreIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) + _label_new_exterior_facets( + dm, subdm, subpoint_indices, + sub_ext_facet_indices, sub_ext_facet_size, + subfStart, subfEnd) elif subdim == dim - 1: - # Codimension-1 submesh. - # Submesh facets are at depth (subdim - 1) in the parent mesh. - # Propagate parent labels at that depth into "Face Sets" on subdm. - if subdim - 1 == 1: - parent_label_cname = edge_sets_label_name - parent_label_pyname = "Edge Sets" - elif subdim - 1 == 0: - parent_label_cname = vertex_sets_label_name - parent_label_pyname = "Vertex Sets" - else: - parent_label_cname = NULL - parent_label_pyname = None - has_parent_facet_label = (parent_label_pyname is not None - and dm.hasLabel(parent_label_pyname)) - next_label_val = 0 - if has_parent_facet_label: - with dm.getLabelIdIS(parent_label_pyname) as label_value_indices: - next_label_val = label_value_indices.max() + 1 if len(label_value_indices) > 0 else 0 - next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) - # DMPlexFilter propagates parent "Face Sets" to the subdm, but those - # values sit on cells (codim-0), not facets. Remove them so they don't - # contaminate the exterior facet markers we are about to create. - if subdm.hasLabel(FACE_SETS_LABEL): - subdm.removeLabel(FACE_SETS_LABEL) - subdm.createLabel(FACE_SETS_LABEL) - sub_ext_facet_size = subdm.getStratumSize("exterior_facets", 1) - sub_ext_facet_is = subdm.getStratumIS("exterior_facets", 1) - if sub_ext_facet_is.iset: - CHKERR(ISGetIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) - subfStart, subfEnd = subdm.getHeightStratum(1) - if has_parent_facet_label: - CHKERR(DMGetLabel(dm.dm, parent_label_cname, &parent_facet_lbl)) - pStart, pEnd = dm.getChart() - CHKERR(DMLabelCreateIndex(parent_facet_lbl, pStart, pEnd)) - for i in range(sub_ext_facet_size): - subf = sub_ext_facet_indices[i] - if subf < subfStart or subf >= subfEnd: - continue - f = subpoint_indices[subf] - label_val = -1 - if has_parent_facet_label: - CHKERR(DMLabelGetValue(parent_facet_lbl, f, &label_val)) - if label_val >= 0: - CHKERR(DMSetLabelValue(subdm.dm, face_sets_label_name, subf, label_val)) - else: - CHKERR(DMSetLabelValue(subdm.dm, face_sets_label_name, subf, next_label_val)) - if has_parent_facet_label: - CHKERR(DMLabelDestroyIndex(parent_facet_lbl)) - if sub_ext_facet_is.iset: - CHKERR(ISRestoreIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) + _propagate_parent_facet_labels( + dm, subdm, subdim, subpoint_indices, + sub_ext_facet_indices, sub_ext_facet_size, + subfStart, subfEnd) + + if sub_ext_facet_is.iset: + CHKERR(ISRestoreIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) CHKERR(ISRestoreIndices(subpoint_is.iset, &subpoint_indices)) subdm.removeLabel("interior_facets") subdm.removeLabel("exterior_facets") diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index d9a9cc4794..26f905b732 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -40,14 +40,14 @@ def _all_at(cs, axis, value): EDGE_C = 3 -def _make_cube_with_edge_sets(extra_edge_c=False): - """UnitCubeMesh(2,2,2) with Face Sets on z=0 boundary faces and Edge Sets +def _make_cube_with_edge_sets(extra_edge_c=False, ncells=2): + """UnitCubeMesh with Face Sets on z=0 boundary faces and Edge Sets on selected boundary edges of the z=0 face. If *extra_edge_c* is True an additional label EDGE_C is placed on edges at z=1, x=0 (not part of the z=0 submesh). """ - mesh = UnitCubeMesh(2, 2, 2) + mesh = UnitCubeMesh(ncells, ncells, ncells) dm = mesh.topology_dm vc = _get_dm_vertex_coords(dm) eStart, eEnd = dm.getDepthStratum(1) @@ -82,11 +82,11 @@ def _make_cube_with_edge_sets(extra_edge_c=False): VTX_B = 2 -def _make_square_with_vertex_sets(): - """UnitSquareMesh(2,2) with Face Sets on y=0 boundary edges and Vertex +def _make_square_with_vertex_sets(ncells=2): + """UnitSquareMesh with Face Sets on y=0 boundary edges and Vertex Sets on the corner vertices of that edge. """ - mesh = UnitSquareMesh(2, 2) + mesh = UnitSquareMesh(ncells, ncells) dm = mesh.topology_dm vc = _get_dm_vertex_coords(dm) vStart, vEnd = dm.getDepthStratum(0) @@ -204,16 +204,132 @@ def test_submesh_codim1_no_parent_edge_sets(): # --------------------------------------------------------------------------- -# Parallel +# Parallel – marker propagation # --------------------------------------------------------------------------- -@pytest.mark.parallel(nprocs=2) -def test_submesh_codim1_edge_sets_propagated_parallel(): - """Edge Sets propagation works in parallel (2 ranks).""" - mesh = _make_cube_with_edge_sets() +def _check_edge_sets_propagated(ncells=2): + mesh = _make_cube_with_edge_sets(ncells=ncells) submesh = Submesh(mesh, 2, FACE_TAG_3D) markers = submesh.exterior_facets.unique_markers assert markers is not None marker_set = set(markers) assert EDGE_A in marker_set, f"EDGE_A={EDGE_A} not in {marker_set}" assert EDGE_B in marker_set, f"EDGE_B={EDGE_B} not in {marker_set}" + + +@pytest.mark.parallel(nprocs=2) +def test_submesh_codim1_edge_sets_propagated_2_procs(): + _check_edge_sets_propagated() + + +@pytest.mark.parallel(nprocs=3) +def test_submesh_codim1_edge_sets_propagated_3_procs(): + _check_edge_sets_propagated(ncells=4) + + +# --------------------------------------------------------------------------- +# Parallel – boundary length integration (ds) +# --------------------------------------------------------------------------- + +def _check_ds_edge_sets(ncells=2): + """Integrate 1 over each tagged boundary of the 2D submesh.""" + mesh = _make_cube_with_edge_sets(ncells=ncells) + submesh = Submesh(mesh, 2, FACE_TAG_3D) + assert abs(assemble(Constant(1.) * ds(EDGE_A, domain=submesh)) - 1.0) < 1e-12 + assert abs(assemble(Constant(1.) * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 + + +@pytest.mark.parallel(nprocs=2) +def test_submesh_codim1_ds_edge_sets_2_procs(): + _check_ds_edge_sets() + + +@pytest.mark.parallel(nprocs=3) +def test_submesh_codim1_ds_edge_sets_3_procs(): + _check_ds_edge_sets(ncells=4) + + +# --------------------------------------------------------------------------- +# Parallel – coordinate expression integration +# --------------------------------------------------------------------------- + +def _check_boundary_coordinate_integral(ncells=2): + """Integrate coordinate expressions over tagged boundaries. + + On the z=0 submesh: + EDGE_A is at x=0, y in [0,1] => int(y dy) = 1/2, int(x dy) = 0 + EDGE_B is at x=1, y in [0,1] => int(y dy) = 1/2, int(x dy) = 1 + """ + mesh = _make_cube_with_edge_sets(ncells=ncells) + submesh = Submesh(mesh, 2, FACE_TAG_3D) + x, y, z = SpatialCoordinate(submesh) + assert abs(assemble(y * ds(EDGE_A, domain=submesh)) - 0.5) < 1e-12 + assert abs(assemble(x * ds(EDGE_A, domain=submesh))) < 1e-12 + assert abs(assemble(y * ds(EDGE_B, domain=submesh)) - 0.5) < 1e-12 + assert abs(assemble(x * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 + + +@pytest.mark.parallel(nprocs=2) +def test_submesh_codim1_boundary_coord_integral_2_procs(): + _check_boundary_coordinate_integral() + + +@pytest.mark.parallel(nprocs=3) +def test_submesh_codim1_boundary_coord_integral_3_procs(): + _check_boundary_coordinate_integral(ncells=4) + + +# --------------------------------------------------------------------------- +# Parallel – FunctionSpace interpolation + boundary integration +# --------------------------------------------------------------------------- + +def _check_interpolate_and_integrate(ncells=2): + """Interpolate x^2+y^2 on the submesh and integrate over tagged boundaries. + + EDGE_A (x=0): int_0^1 y^2 dy = 1/3 + EDGE_B (x=1): int_0^1 (1+y^2) dy = 4/3 + """ + mesh = _make_cube_with_edge_sets(ncells=ncells) + submesh = Submesh(mesh, 2, FACE_TAG_3D) + V = FunctionSpace(submesh, "CG", 2) + x, y, z = SpatialCoordinate(submesh) + f = Function(V).interpolate(x**2 + y**2) + assert abs(assemble(f * ds(EDGE_A)) - 1. / 3) < 1e-12 + assert abs(assemble(f * ds(EDGE_B)) - 4. / 3) < 1e-12 + + +@pytest.mark.parallel(nprocs=2) +def test_submesh_codim1_interpolate_integrate_2_procs(): + _check_interpolate_and_integrate() + + +@pytest.mark.parallel(nprocs=3) +def test_submesh_codim1_interpolate_integrate_3_procs(): + _check_interpolate_and_integrate(ncells=4) + + +# --------------------------------------------------------------------------- +# Parallel – 2D -> 1D (Vertex Sets) +# --------------------------------------------------------------------------- + +def _check_vertex_sets_ds(ncells=2): + """Integrate over vertex-tagged boundaries of a 1D submesh.""" + mesh = _make_square_with_vertex_sets(ncells=ncells) + submesh = Submesh(mesh, 1, FACE_TAG_2D) + markers = submesh.exterior_facets.unique_markers + assert markers is not None + marker_set = set(markers) + assert VTX_A in marker_set + assert VTX_B in marker_set + assert abs(assemble(Constant(1.) * ds(VTX_A, domain=submesh)) - 1.0) < 1e-12 + assert abs(assemble(Constant(1.) * ds(VTX_B, domain=submesh)) - 1.0) < 1e-12 + + +@pytest.mark.parallel(nprocs=2) +def test_submesh_codim1_vertex_sets_ds_2_procs(): + _check_vertex_sets_ds() + + +@pytest.mark.parallel(nprocs=3) +def test_submesh_codim1_vertex_sets_ds_3_procs(): + _check_vertex_sets_ds(ncells=4) From fcbca356022c765cb9a4cd9923a8873c8c422ab0 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 09:57:31 +0100 Subject: [PATCH 05/15] fix tests --- firedrake/cython/dmcommon.pyx | 31 +++++++++++++------ .../submesh/test_submesh_codim1_labels.py | 25 ++++++++++++--- 2 files changed, 41 insertions(+), 15 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index e38e8767fb..d6db3fb19d 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -3983,12 +3983,13 @@ cdef void _propagate_parent_facet_labels( """Codimension-1 helper: map the parent's lower-dimensional label into "Face Sets". For a 3D→2D submesh the parent's "Edge Sets" are propagated; for 2D→1D - the parent's "Vertex Sets" are used. Exterior facets without a - corresponding parent label value receive a default tag. + the parent's "Vertex Sets" are used. Exterior facets that already carry + an inherited "Face Sets" value (from ``DMPlexFilter``) are left alone; + the remaining ones receive the parent label value or a default tag. """ cdef: - PetscInt pStart, pEnd, next_label_val, label_val, subf, f, i - DMLabel parent_label + PetscInt pStart, pEnd, next_label_val, label_val, existing_val, subf, f, i + DMLabel parent_label, face_sets_label if subdim == 2: parent_label_name = b"Edge Sets" @@ -4003,14 +4004,20 @@ cdef void _propagate_parent_facet_labels( next_label_val = 0 if has_parent_label: with dm.getLabelIdIS(parent_label_name) as ids: - next_label_val = ids.max() + 1 if len(ids) > 0 else 0 + if len(ids) > 0: + next_label_val = max(next_label_val, ids.max() + 1) + if subdm.hasLabel(FACE_SETS_LABEL): + with subdm.getLabelIdIS(FACE_SETS_LABEL) as ids: + if len(ids) > 0: + next_label_val = max(next_label_val, ids.max() + 1) next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) - # DMPlexFilter may propagate cell-level "Face Sets" to the subdm; - # remove them so they don't contaminate the facet markers. - if subdm.hasLabel(FACE_SETS_LABEL): - subdm.removeLabel(FACE_SETS_LABEL) - subdm.createLabel(FACE_SETS_LABEL) + # Ensure "Face Sets" exists but keep any values inherited from DMPlexFilter + # (e.g. mark_entities labels all closure entities, so boundary edges of a + # codim-1 submesh may already carry the correct tags). + if not subdm.hasLabel(FACE_SETS_LABEL): + subdm.createLabel(FACE_SETS_LABEL) + CHKERR(DMGetLabel(subdm.dm, b"Face Sets", &face_sets_label)) if has_parent_label: CHKERR(DMGetLabel(dm.dm, parent_label_name, &parent_label)) @@ -4019,6 +4026,10 @@ cdef void _propagate_parent_facet_labels( for i in range(sub_ext_facet_size): subf = sub_ext_facet_indices[i] if subfStart <= subf < subfEnd: + # Skip facets that already have an inherited "Face Sets" value. + CHKERR(DMLabelGetValue(face_sets_label, subf, &existing_val)) + if existing_val >= 0: + continue f = subpoint_indices[subf] label_val = -1 if has_parent_label: diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index 26f905b732..36a98dae32 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -49,6 +49,10 @@ def _make_cube_with_edge_sets(extra_edge_c=False, ncells=2): """ mesh = UnitCubeMesh(ncells, ncells, ncells) dm = mesh.topology_dm + # Remove default Face Sets so we have full control over labels. + if dm.hasLabel("Face Sets"): + dm.removeLabel("Face Sets") + dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) eStart, eEnd = dm.getDepthStratum(1) fStart, fEnd = dm.getDepthStratum(2) @@ -88,6 +92,10 @@ def _make_square_with_vertex_sets(ncells=2): """ mesh = UnitSquareMesh(ncells, ncells) dm = mesh.topology_dm + # Remove default Face Sets so we have full control over labels. + if dm.hasLabel("Face Sets"): + dm.removeLabel("Face Sets") + dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) vStart, vEnd = dm.getDepthStratum(0) eStart, eEnd = dm.getDepthStratum(1) @@ -136,14 +144,19 @@ def test_submesh_codim1_edge_sets_excludes_nonsubmesh(): def test_submesh_codim1_unlabeled_get_default_value(): - """Exterior facets without a parent Edge Sets label get max(label_vals)+1.""" + """Exterior facets without a parent Edge Sets label get a fresh default value.""" mesh = _make_cube_with_edge_sets() submesh = Submesh(mesh, 2, FACE_TAG_3D) markers = submesh.exterior_facets.unique_markers assert markers is not None marker_set = set(markers) - next_val = max(EDGE_A, EDGE_B) + 1 - assert next_val in marker_set, f"Default {next_val} not in {marker_set}" + # Default must be larger than all Edge Sets and inherited Face Sets values. + default_vals = marker_set - {EDGE_A, EDGE_B, FACE_TAG_3D} + assert len(default_vals) == 1, f"Expected one default value, got {default_vals}" + default_val = default_vals.pop() + assert default_val > max(EDGE_A, EDGE_B), ( + f"Default {default_val} should be > max(EDGE_A, EDGE_B)={max(EDGE_A, EDGE_B)}" + ) def test_submesh_codim1_ds_edge_sets(): @@ -188,6 +201,9 @@ def test_submesh_codim1_no_parent_edge_sets(): """Codim-1 submesh with no Edge Sets on the parent still labels exterior facets.""" mesh = UnitCubeMesh(2, 2, 2) dm = mesh.topology_dm + if dm.hasLabel("Face Sets"): + dm.removeLabel("Face Sets") + dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) fStart, fEnd = dm.getDepthStratum(2) @@ -199,8 +215,7 @@ def test_submesh_codim1_no_parent_edge_sets(): submesh = Submesh(mesh, 2, FACE_TAG_3D) markers = submesh.exterior_facets.unique_markers assert markers is not None - # All exterior facets should receive the default label (0) - assert 0 in set(markers) + assert len(markers) > 0 # --------------------------------------------------------------------------- From 88446cd8336d82c8fa434eae6b598df77d218e3a Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 10:09:29 +0100 Subject: [PATCH 06/15] small refactoring --- firedrake/cython/dmcommon.pyx | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index d6db3fb19d..4c5b9e304d 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -3937,6 +3937,14 @@ def submesh_correct_entity_classes(PETSc.DM dm, CHKERR(DMLabelDestroyIndex(lbl_ghost)) +cdef PetscInt _max_label_value(PETSc.DM dm, label_name): + """Return the maximum value in *label_name*, or -1 if absent/empty.""" + if not dm.hasLabel(label_name): + return -1 + with dm.getLabelIdIS(label_name) as ids: + return ids.max() if len(ids) > 0 else -1 + + cdef void _label_new_exterior_facets( PETSc.DM dm, PETSc.DM subdm, const PetscInt *subpoint_indices, @@ -3954,8 +3962,7 @@ cdef void _label_new_exterior_facets( DMLabel parent_ext_label PetscBool has_point - with dm.getLabelIdIS(FACE_SETS_LABEL) as ids: - next_label_val = ids.max() + 1 if len(ids) > 0 else 0 + next_label_val = _max_label_value(dm, FACE_SETS_LABEL) + 1 next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) subdm.createLabel(FACE_SETS_LABEL) @@ -4001,15 +4008,10 @@ cdef void _propagate_parent_facet_labels( has_parent_label = (parent_label_name is not None and dm.hasLabel(parent_label_name)) - next_label_val = 0 - if has_parent_label: - with dm.getLabelIdIS(parent_label_name) as ids: - if len(ids) > 0: - next_label_val = max(next_label_val, ids.max() + 1) - if subdm.hasLabel(FACE_SETS_LABEL): - with subdm.getLabelIdIS(FACE_SETS_LABEL) as ids: - if len(ids) > 0: - next_label_val = max(next_label_val, ids.max() + 1) + next_label_val = max( + _max_label_value(dm, parent_label_name) if has_parent_label else -1, + _max_label_value(subdm, FACE_SETS_LABEL), + ) + 1 next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) # Ensure "Face Sets" exists but keep any values inherited from DMPlexFilter From eec6bed899ad35954300d868179b63b7faf1ebad Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 13:21:11 +0100 Subject: [PATCH 07/15] Update tests/firedrake/submesh/test_submesh_codim1_labels.py Co-authored-by: Pablo Brubeck --- tests/firedrake/submesh/test_submesh_codim1_labels.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index 36a98dae32..422e533702 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -232,13 +232,8 @@ def _check_edge_sets_propagated(ncells=2): assert EDGE_B in marker_set, f"EDGE_B={EDGE_B} not in {marker_set}" -@pytest.mark.parallel(nprocs=2) -def test_submesh_codim1_edge_sets_propagated_2_procs(): - _check_edge_sets_propagated() - - -@pytest.mark.parallel(nprocs=3) -def test_submesh_codim1_edge_sets_propagated_3_procs(): +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_edge_sets_propagated(): _check_edge_sets_propagated(ncells=4) From 2cfffedcbf6c12e9600c86a6b1c01ac58a16c528 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 13:21:34 +0100 Subject: [PATCH 08/15] Update tests/firedrake/submesh/test_submesh_codim1_labels.py Co-authored-by: Pablo Brubeck --- tests/firedrake/submesh/test_submesh_codim1_labels.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index 422e533702..d96c2da3a6 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -249,13 +249,8 @@ def _check_ds_edge_sets(ncells=2): assert abs(assemble(Constant(1.) * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 -@pytest.mark.parallel(nprocs=2) -def test_submesh_codim1_ds_edge_sets_2_procs(): - _check_ds_edge_sets() - - -@pytest.mark.parallel(nprocs=3) -def test_submesh_codim1_ds_edge_sets_3_procs(): +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_ds_edge_sets(): _check_ds_edge_sets(ncells=4) From 7507425179a9ccf3705f62c06ac0469d8edc69ba Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 13:22:07 +0100 Subject: [PATCH 09/15] Apply suggestions from code review Co-authored-by: Pablo Brubeck --- .../submesh/test_submesh_codim1_labels.py | 25 ++++--------------- 1 file changed, 5 insertions(+), 20 deletions(-) diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index d96c2da3a6..1708eb6a4e 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -274,13 +274,8 @@ def _check_boundary_coordinate_integral(ncells=2): assert abs(assemble(x * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 -@pytest.mark.parallel(nprocs=2) -def test_submesh_codim1_boundary_coord_integral_2_procs(): - _check_boundary_coordinate_integral() - - -@pytest.mark.parallel(nprocs=3) -def test_submesh_codim1_boundary_coord_integral_3_procs(): +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_boundary_coord_integral(): _check_boundary_coordinate_integral(ncells=4) @@ -303,13 +298,8 @@ def _check_interpolate_and_integrate(ncells=2): assert abs(assemble(f * ds(EDGE_B)) - 4. / 3) < 1e-12 -@pytest.mark.parallel(nprocs=2) -def test_submesh_codim1_interpolate_integrate_2_procs(): - _check_interpolate_and_integrate() - - -@pytest.mark.parallel(nprocs=3) -def test_submesh_codim1_interpolate_integrate_3_procs(): +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_interpolate_integrate(): _check_interpolate_and_integrate(ncells=4) @@ -330,11 +320,6 @@ def _check_vertex_sets_ds(ncells=2): assert abs(assemble(Constant(1.) * ds(VTX_B, domain=submesh)) - 1.0) < 1e-12 -@pytest.mark.parallel(nprocs=2) -def test_submesh_codim1_vertex_sets_ds_2_procs(): - _check_vertex_sets_ds() - - -@pytest.mark.parallel(nprocs=3) +@pytest.mark.parallel([1, 2, 3]) def test_submesh_codim1_vertex_sets_ds_3_procs(): _check_vertex_sets_ds(ncells=4) From 307c68870271981912a0913066fde7d046513c8e Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 15:26:35 +0100 Subject: [PATCH 10/15] lint fixes --- .../submesh/test_submesh_codim1_labels.py | 37 ++++++++----------- 1 file changed, 15 insertions(+), 22 deletions(-) diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index 36a98dae32..566395efdd 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -120,9 +120,8 @@ def _make_square_with_vertex_sets(ncells=2): # Tests – 3D parent -> 2D submesh (Edge Sets → Face Sets) # --------------------------------------------------------------------------- -def test_submesh_codim1_edge_sets_propagated(): - """Parent 'Edge Sets' tags appear in submesh exterior_facets.unique_markers.""" - mesh = _make_cube_with_edge_sets() +def _check_edge_sets_propagated(ncells=2): + mesh = _make_cube_with_edge_sets(ncells=ncells) submesh = Submesh(mesh, 2, FACE_TAG_3D) markers = submesh.exterior_facets.unique_markers assert markers is not None @@ -131,6 +130,11 @@ def test_submesh_codim1_edge_sets_propagated(): assert EDGE_B in marker_set, f"EDGE_B={EDGE_B} not in {marker_set}" +def test_submesh_codim1_edge_sets_propagated(): + """Parent 'Edge Sets' tags appear in submesh exterior_facets.unique_markers.""" + _check_edge_sets_propagated() + + def test_submesh_codim1_edge_sets_excludes_nonsubmesh(): """Edge Sets on parent edges outside the submesh are excluded.""" mesh = _make_cube_with_edge_sets(extra_edge_c=True) @@ -159,9 +163,9 @@ def test_submesh_codim1_unlabeled_get_default_value(): ) -def test_submesh_codim1_ds_edge_sets(): - """ds(tag) with propagated Edge Sets integrates the correct boundary length.""" - mesh = _make_cube_with_edge_sets() +def _check_ds_edge_sets(ncells=2): + """Integrate 1 over each tagged boundary of the 2D submesh.""" + mesh = _make_cube_with_edge_sets(ncells=ncells) submesh = Submesh(mesh, 2, FACE_TAG_3D) # EDGE_A labels edges at x=0, z=0: total length = 1 assert abs(assemble(Constant(1.) * ds(EDGE_A, domain=submesh)) - 1.0) < 1e-12 @@ -169,6 +173,11 @@ def test_submesh_codim1_ds_edge_sets(): assert abs(assemble(Constant(1.) * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 +def test_submesh_codim1_ds_edge_sets(): + """ds(tag) with propagated Edge Sets integrates the correct boundary length.""" + _check_ds_edge_sets() + + # --------------------------------------------------------------------------- # Tests – 2D parent -> 1D submesh (Vertex Sets → Face Sets) # --------------------------------------------------------------------------- @@ -222,15 +231,6 @@ def test_submesh_codim1_no_parent_edge_sets(): # Parallel – marker propagation # --------------------------------------------------------------------------- -def _check_edge_sets_propagated(ncells=2): - mesh = _make_cube_with_edge_sets(ncells=ncells) - submesh = Submesh(mesh, 2, FACE_TAG_3D) - markers = submesh.exterior_facets.unique_markers - assert markers is not None - marker_set = set(markers) - assert EDGE_A in marker_set, f"EDGE_A={EDGE_A} not in {marker_set}" - assert EDGE_B in marker_set, f"EDGE_B={EDGE_B} not in {marker_set}" - @pytest.mark.parallel(nprocs=2) def test_submesh_codim1_edge_sets_propagated_2_procs(): @@ -246,13 +246,6 @@ def test_submesh_codim1_edge_sets_propagated_3_procs(): # Parallel – boundary length integration (ds) # --------------------------------------------------------------------------- -def _check_ds_edge_sets(ncells=2): - """Integrate 1 over each tagged boundary of the 2D submesh.""" - mesh = _make_cube_with_edge_sets(ncells=ncells) - submesh = Submesh(mesh, 2, FACE_TAG_3D) - assert abs(assemble(Constant(1.) * ds(EDGE_A, domain=submesh)) - 1.0) < 1e-12 - assert abs(assemble(Constant(1.) * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 - @pytest.mark.parallel(nprocs=2) def test_submesh_codim1_ds_edge_sets_2_procs(): From b2bdc653f5b080bc4d529eabc7413f02696287d6 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Tue, 10 Mar 2026 18:10:26 +0100 Subject: [PATCH 11/15] Update firedrake/cython/dmcommon.pyx Co-authored-by: Connor Ward --- 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 4c5b9e304d..5211e89764 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -4062,7 +4062,7 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): were interior in the parent) are tagged with ``max("Face Sets") + 1``. * **Codimension-1** (``subdim == dim - 1``): the parent's lower-dimensional labels ("Edge Sets" for 3D, "Vertex Sets" for 2D) are mapped into - "Face Sets" on the subdm. Unlabeled facets get a default value. + "Face Sets" on the subdm. Unlabeled facets get a default value of ``max(label) + 1``. """ cdef: From 562934616a4608ed2de6b09f3df88a89ff521347 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Wed, 11 Mar 2026 08:25:52 +0100 Subject: [PATCH 12/15] improve a docstring --- firedrake/cython/dmcommon.pyx | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 4c5b9e304d..b590583b9b 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -3990,9 +3990,17 @@ cdef void _propagate_parent_facet_labels( """Codimension-1 helper: map the parent's lower-dimensional label into "Face Sets". For a 3D→2D submesh the parent's "Edge Sets" are propagated; for 2D→1D - the parent's "Vertex Sets" are used. Exterior facets that already carry - an inherited "Face Sets" value (from ``DMPlexFilter``) are left alone; - the remaining ones receive the parent label value or a default tag. + the parent's "Vertex Sets" are used. + + ``DMPlexFilter`` copies every label from the parent into the subdm. + In particular, the inherited "Face Sets" already carries useful + facet-level boundary information (from ``complete_facet_labels`` / + ``DMPlexLabelComplete`` on the parent). We preserve those inherited + values because they are required by downstream code (e.g. a subsequent + ``Submesh`` call that selects by "Face Sets" value). Only exterior + facets that do *not* already have an inherited "Face Sets" value are + updated: they receive the parent's lower-dimensional label value, or a + fresh default tag if none exists. """ cdef: PetscInt pStart, pEnd, next_label_val, label_val, existing_val, subf, f, i @@ -4014,9 +4022,6 @@ cdef void _propagate_parent_facet_labels( ) + 1 next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) - # Ensure "Face Sets" exists but keep any values inherited from DMPlexFilter - # (e.g. mark_entities labels all closure entities, so boundary edges of a - # codim-1 submesh may already carry the correct tags). if not subdm.hasLabel(FACE_SETS_LABEL): subdm.createLabel(FACE_SETS_LABEL) CHKERR(DMGetLabel(subdm.dm, b"Face Sets", &face_sets_label)) @@ -4028,7 +4033,6 @@ cdef void _propagate_parent_facet_labels( for i in range(sub_ext_facet_size): subf = sub_ext_facet_indices[i] if subfStart <= subf < subfEnd: - # Skip facets that already have an inherited "Face Sets" value. CHKERR(DMLabelGetValue(face_sets_label, subf, &existing_val)) if existing_val >= 0: continue From e3f0b85e029a2dcddeb66cfbfff7f261974a0423 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Wed, 11 Mar 2026 08:30:59 +0100 Subject: [PATCH 13/15] refactor tests --- .../submesh/test_submesh_codim1_labels.py | 73 +++++-------------- 1 file changed, 17 insertions(+), 56 deletions(-) diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index aaf3c49343..f2640d9808 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -120,8 +120,10 @@ def _make_square_with_vertex_sets(ncells=2): # Tests – 3D parent -> 2D submesh (Edge Sets → Face Sets) # --------------------------------------------------------------------------- -def _check_edge_sets_propagated(ncells=2): - mesh = _make_cube_with_edge_sets(ncells=ncells) +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_edge_sets_propagated(): + """Parent 'Edge Sets' tags appear in submesh exterior_facets.unique_markers.""" + mesh = _make_cube_with_edge_sets(ncells=4) submesh = Submesh(mesh, 2, FACE_TAG_3D) markers = submesh.exterior_facets.unique_markers assert markers is not None @@ -130,11 +132,6 @@ def _check_edge_sets_propagated(ncells=2): assert EDGE_B in marker_set, f"EDGE_B={EDGE_B} not in {marker_set}" -def test_submesh_codim1_edge_sets_propagated(): - """Parent 'Edge Sets' tags appear in submesh exterior_facets.unique_markers.""" - _check_edge_sets_propagated() - - def test_submesh_codim1_edge_sets_excludes_nonsubmesh(): """Edge Sets on parent edges outside the submesh are excluded.""" mesh = _make_cube_with_edge_sets(extra_edge_c=True) @@ -163,9 +160,10 @@ def test_submesh_codim1_unlabeled_get_default_value(): ) -def _check_ds_edge_sets(ncells=2): - """Integrate 1 over each tagged boundary of the 2D submesh.""" - mesh = _make_cube_with_edge_sets(ncells=ncells) +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_ds_edge_sets(): + """ds(tag) with propagated Edge Sets integrates the correct boundary length.""" + mesh = _make_cube_with_edge_sets(ncells=4) submesh = Submesh(mesh, 2, FACE_TAG_3D) # EDGE_A labels edges at x=0, z=0: total length = 1 assert abs(assemble(Constant(1.) * ds(EDGE_A, domain=submesh)) - 1.0) < 1e-12 @@ -173,11 +171,6 @@ def _check_ds_edge_sets(ncells=2): assert abs(assemble(Constant(1.) * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 -def test_submesh_codim1_ds_edge_sets(): - """ds(tag) with propagated Edge Sets integrates the correct boundary length.""" - _check_ds_edge_sets() - - # --------------------------------------------------------------------------- # Tests – 2D parent -> 1D submesh (Vertex Sets → Face Sets) # --------------------------------------------------------------------------- @@ -227,38 +220,19 @@ def test_submesh_codim1_no_parent_edge_sets(): assert len(markers) > 0 -# --------------------------------------------------------------------------- -# Parallel – marker propagation -# --------------------------------------------------------------------------- - - -@pytest.mark.parallel([1, 2, 3]) -def test_submesh_codim1_edge_sets_propagated(): - _check_edge_sets_propagated(ncells=4) - - -# --------------------------------------------------------------------------- -# Parallel – boundary length integration (ds) -# --------------------------------------------------------------------------- - - -@pytest.mark.parallel([1, 2, 3]) -def test_submesh_codim1_ds_edge_sets(): - _check_ds_edge_sets(ncells=4) - - # --------------------------------------------------------------------------- # Parallel – coordinate expression integration # --------------------------------------------------------------------------- -def _check_boundary_coordinate_integral(ncells=2): +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_boundary_coord_integral(): """Integrate coordinate expressions over tagged boundaries. On the z=0 submesh: EDGE_A is at x=0, y in [0,1] => int(y dy) = 1/2, int(x dy) = 0 EDGE_B is at x=1, y in [0,1] => int(y dy) = 1/2, int(x dy) = 1 """ - mesh = _make_cube_with_edge_sets(ncells=ncells) + mesh = _make_cube_with_edge_sets(ncells=4) submesh = Submesh(mesh, 2, FACE_TAG_3D) x, y, z = SpatialCoordinate(submesh) assert abs(assemble(y * ds(EDGE_A, domain=submesh)) - 0.5) < 1e-12 @@ -267,22 +241,18 @@ def _check_boundary_coordinate_integral(ncells=2): assert abs(assemble(x * ds(EDGE_B, domain=submesh)) - 1.0) < 1e-12 -@pytest.mark.parallel([1, 2, 3]) -def test_submesh_codim1_boundary_coord_integral(): - _check_boundary_coordinate_integral(ncells=4) - - # --------------------------------------------------------------------------- # Parallel – FunctionSpace interpolation + boundary integration # --------------------------------------------------------------------------- -def _check_interpolate_and_integrate(ncells=2): +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_interpolate_integrate(): """Interpolate x^2+y^2 on the submesh and integrate over tagged boundaries. EDGE_A (x=0): int_0^1 y^2 dy = 1/3 EDGE_B (x=1): int_0^1 (1+y^2) dy = 4/3 """ - mesh = _make_cube_with_edge_sets(ncells=ncells) + mesh = _make_cube_with_edge_sets(ncells=4) submesh = Submesh(mesh, 2, FACE_TAG_3D) V = FunctionSpace(submesh, "CG", 2) x, y, z = SpatialCoordinate(submesh) @@ -291,18 +261,14 @@ def _check_interpolate_and_integrate(ncells=2): assert abs(assemble(f * ds(EDGE_B)) - 4. / 3) < 1e-12 -@pytest.mark.parallel([1, 2, 3]) -def test_submesh_codim1_interpolate_integrate(): - _check_interpolate_and_integrate(ncells=4) - - # --------------------------------------------------------------------------- # Parallel – 2D -> 1D (Vertex Sets) # --------------------------------------------------------------------------- -def _check_vertex_sets_ds(ncells=2): +@pytest.mark.parallel([1, 2, 3]) +def test_submesh_codim1_vertex_sets_ds(): """Integrate over vertex-tagged boundaries of a 1D submesh.""" - mesh = _make_square_with_vertex_sets(ncells=ncells) + mesh = _make_square_with_vertex_sets(ncells=4) submesh = Submesh(mesh, 1, FACE_TAG_2D) markers = submesh.exterior_facets.unique_markers assert markers is not None @@ -311,8 +277,3 @@ def _check_vertex_sets_ds(ncells=2): assert VTX_B in marker_set assert abs(assemble(Constant(1.) * ds(VTX_A, domain=submesh)) - 1.0) < 1e-12 assert abs(assemble(Constant(1.) * ds(VTX_B, domain=submesh)) - 1.0) < 1e-12 - - -@pytest.mark.parallel([1, 2, 3]) -def test_submesh_codim1_vertex_sets_ds_3_procs(): - _check_vertex_sets_ds(ncells=4) From 5f304028010aefa8a0fdb9d916d30bb1f5f47d09 Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Wed, 11 Mar 2026 14:26:24 +0100 Subject: [PATCH 14/15] connor's suggestion --- firedrake/cython/dmcommon.pyx | 149 ++++++++++++------ .../submesh/test_submesh_codim1_labels.py | 19 ++- 2 files changed, 115 insertions(+), 53 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index aa6a306c5f..ade1fb63ae 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -24,6 +24,8 @@ include "petschdr.pxi" FACE_SETS_LABEL = "Face Sets" CELL_SETS_LABEL = "Cell Sets" +EDGE_SETS_LABEL = "Edge Sets" +VERTEX_SETS_LABEL = "Vertex Sets" class DistributedMeshOverlapType(enum.Enum): @@ -1724,18 +1726,75 @@ def label_facets(PETSc.DM plex): CHKERR(DMLabelDestroyIndex(lbl_ext)) +cdef void _populate_lower_dim_labels(PETSc.DM dm, PetscInt tdim): + """Copy "Face Sets" values to "Edge Sets" and "Vertex Sets" for + lower-dimensional entities. + + After ``DMPlexLabelComplete`` has propagated "Face Sets" into the + closure, edges and vertices carry the boundary tag of their parent + face. This helper mirrors those values into the corresponding + lower-dimensional label so that a codimension-1 ``Submesh`` can + find the information it needs without relying on the inherited + "Face Sets" (which mixes cell- and facet-level values). + + Entities that already have a value in the target label are skipped + to avoid overwriting user-defined labels (e.g. from Gmsh). + """ + cdef: + PetscInt pStart, pEnd, p, face_val, existing_val + DMLabel face_sets_lbl, edge_sets_lbl, vertex_sets_lbl + + CHKERR(DMGetLabel(dm.dm, b"Face Sets", &face_sets_lbl)) + + if tdim >= 3: + if not dm.hasLabel(EDGE_SETS_LABEL): + dm.createLabel(EDGE_SETS_LABEL) + CHKERR(DMGetLabel(dm.dm, b"Edge Sets", &edge_sets_lbl)) + CHKERR(DMPlexGetDepthStratum(dm.dm, 1, &pStart, &pEnd)) + for p in range(pStart, pEnd): + CHKERR(DMLabelGetValue(face_sets_lbl, p, &face_val)) + if face_val >= 0: + CHKERR(DMLabelGetValue(edge_sets_lbl, p, &existing_val)) + if existing_val < 0: + CHKERR(DMLabelSetValue(edge_sets_lbl, p, face_val)) + + if tdim >= 2: + if not dm.hasLabel(VERTEX_SETS_LABEL): + dm.createLabel(VERTEX_SETS_LABEL) + CHKERR(DMGetLabel(dm.dm, b"Vertex Sets", &vertex_sets_lbl)) + CHKERR(DMPlexGetDepthStratum(dm.dm, 0, &pStart, &pEnd)) + for p in range(pStart, pEnd): + CHKERR(DMLabelGetValue(face_sets_lbl, p, &face_val)) + if face_val >= 0: + CHKERR(DMLabelGetValue(vertex_sets_lbl, p, &existing_val)) + if existing_val < 0: + CHKERR(DMLabelSetValue(vertex_sets_lbl, p, face_val)) + + def complete_facet_labels(PETSc.DM dm): """Transfer label values from the facet labels to everything in - the closure of the facets.""" - cdef PETSc.DMLabel label + the closure of the facets. - if get_topological_dimension(dm) == 0: + After completing "Face Sets" this also populates "Edge Sets" and + "Vertex Sets" so that boundary information is available at every + topological dimension. Entities that already carry a value in the + target label (e.g. user-defined labels from Gmsh) are not overwritten. + """ + cdef: + PETSc.DMLabel label + PetscInt tdim + + tdim = get_topological_dimension(dm) + if tdim == 0: return for name in [FACE_SETS_LABEL, "exterior_facets", "interior_facets"]: if dm.hasLabel(name): label = dm.getLabel(name) CHKERR( DMPlexLabelComplete(dm.dm, label.dmlabel) ) + if dm.hasLabel(FACE_SETS_LABEL): + _populate_lower_dim_labels(dm, tdim) + @cython.boundscheck(False) @cython.wraparound(False) @@ -3980,72 +4039,71 @@ cdef void _label_new_exterior_facets( cdef void _propagate_parent_facet_labels( - PETSc.DM dm, PETSc.DM subdm, + PETSc.DM subdm, PetscInt subdim, - const PetscInt *subpoint_indices, const PetscInt *sub_ext_facet_indices, PetscInt sub_ext_facet_size, PetscInt subfStart, PetscInt subfEnd, ): - """Codimension-1 helper: map the parent's lower-dimensional label into "Face Sets". - - For a 3D→2D submesh the parent's "Edge Sets" are propagated; for 2D→1D - the parent's "Vertex Sets" are used. - - ``DMPlexFilter`` copies every label from the parent into the subdm. - In particular, the inherited "Face Sets" already carries useful - facet-level boundary information (from ``complete_facet_labels`` / - ``DMPlexLabelComplete`` on the parent). We preserve those inherited - values because they are required by downstream code (e.g. a subsequent - ``Submesh`` call that selects by "Face Sets" value). Only exterior - facets that do *not* already have an inherited "Face Sets" value are - updated: they receive the parent's lower-dimensional label value, or a - fresh default tag if none exists. + """Codimension-1 helper: populate "Face Sets" from the subdm's own + lower-dimensional labels. + + ``DMPlexFilter`` copies every label from the parent into the subdm, + including "Edge Sets" and "Vertex Sets". Inherited "Face Sets" + values are preserved because removing and recreating the label + loses PETSc's internal parallel label-migration state, breaking + chained submesh operations (e.g. 3D → 2D → 1D). + + Only exterior facets that do *not* already have an inherited + "Face Sets" value are updated: they receive the subdm's own + "Edge Sets" (3D→2D) or "Vertex Sets" (2D→1D) value, or a fresh + default tag if none exists. """ cdef: - PetscInt pStart, pEnd, next_label_val, label_val, existing_val, subf, f, i - DMLabel parent_label, face_sets_label + PetscInt pStart, pEnd, next_label_val, label_val, existing_val, subf, i + DMLabel source_label, face_sets_label if subdim == 2: - parent_label_name = b"Edge Sets" + source_label_name = b"Edge Sets" elif subdim == 1: - parent_label_name = b"Vertex Sets" + source_label_name = b"Vertex Sets" else: - parent_label_name = None + source_label_name = None - has_parent_label = (parent_label_name is not None - and dm.hasLabel(parent_label_name)) + has_source = (source_label_name is not None + and subdm.hasLabel(source_label_name)) next_label_val = max( - _max_label_value(dm, parent_label_name) if has_parent_label else -1, + _max_label_value(subdm, source_label_name) if has_source else -1, _max_label_value(subdm, FACE_SETS_LABEL), ) + 1 - next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) + next_label_val = subdm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) if not subdm.hasLabel(FACE_SETS_LABEL): subdm.createLabel(FACE_SETS_LABEL) CHKERR(DMGetLabel(subdm.dm, b"Face Sets", &face_sets_label)) - if has_parent_label: - CHKERR(DMGetLabel(dm.dm, parent_label_name, &parent_label)) - pStart, pEnd = dm.getChart() - CHKERR(DMLabelCreateIndex(parent_label, pStart, pEnd)) + if has_source: + CHKERR(DMGetLabel(subdm.dm, source_label_name, &source_label)) + pStart, pEnd = subdm.getChart() + CHKERR(DMLabelCreateIndex(source_label, pStart, pEnd)) + for i in range(sub_ext_facet_size): subf = sub_ext_facet_indices[i] if subfStart <= subf < subfEnd: CHKERR(DMLabelGetValue(face_sets_label, subf, &existing_val)) if existing_val >= 0: continue - f = subpoint_indices[subf] label_val = -1 - if has_parent_label: - CHKERR(DMLabelGetValue(parent_label, f, &label_val)) + if has_source: + CHKERR(DMLabelGetValue(source_label, subf, &label_val)) if label_val >= 0: - CHKERR(DMSetLabelValue(subdm.dm, b"Face Sets", subf, label_val)) + CHKERR(DMLabelSetValue(face_sets_label, subf, label_val)) else: CHKERR(DMSetLabelValue(subdm.dm, b"Face Sets", subf, next_label_val)) - if has_parent_label: - CHKERR(DMLabelDestroyIndex(parent_label)) + + if has_source: + CHKERR(DMLabelDestroyIndex(source_label)) @cython.boundscheck(False) @@ -4064,9 +4122,10 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): ----- * **Same-dimension** (``subdim == dim``): new exterior facets (those that were interior in the parent) are tagged with ``max("Face Sets") + 1``. - * **Codimension-1** (``subdim == dim - 1``): the parent's lower-dimensional - labels ("Edge Sets" for 3D, "Vertex Sets" for 2D) are mapped into - "Face Sets" on the subdm. Unlabeled facets get a default value of ``max(label) + 1``. + * **Codimension-1** (``subdim == dim - 1``): inherited "Face Sets" values + are preserved and exterior facets without a value are filled from the + subdm's own "Edge Sets" (3D→2D) or "Vertex Sets" (2D→1D). + Unlabeled facets get a default value of ``max(label) + 1``. """ cdef: @@ -4082,8 +4141,6 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): label_facets(subdm) - subpoint_is = subdm.getSubpointIS() - CHKERR(ISGetIndices(subpoint_is.iset, &subpoint_indices)) sub_ext_facet_size = subdm.getStratumSize("exterior_facets", 1) sub_ext_facet_is = subdm.getStratumIS("exterior_facets", 1) if sub_ext_facet_is.iset: @@ -4091,19 +4148,21 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): subfStart, subfEnd = subdm.getHeightStratum(1) if subdim == dim: + subpoint_is = subdm.getSubpointIS() + CHKERR(ISGetIndices(subpoint_is.iset, &subpoint_indices)) _label_new_exterior_facets( dm, subdm, subpoint_indices, sub_ext_facet_indices, sub_ext_facet_size, subfStart, subfEnd) + CHKERR(ISRestoreIndices(subpoint_is.iset, &subpoint_indices)) elif subdim == dim - 1: _propagate_parent_facet_labels( - dm, subdm, subdim, subpoint_indices, + subdm, subdim, sub_ext_facet_indices, sub_ext_facet_size, subfStart, subfEnd) if sub_ext_facet_is.iset: CHKERR(ISRestoreIndices(sub_ext_facet_is.iset, &sub_ext_facet_indices)) - CHKERR(ISRestoreIndices(subpoint_is.iset, &subpoint_indices)) subdm.removeLabel("interior_facets") subdm.removeLabel("exterior_facets") diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index f2640d9808..41dced0ee2 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -49,9 +49,10 @@ def _make_cube_with_edge_sets(extra_edge_c=False, ncells=2): """ mesh = UnitCubeMesh(ncells, ncells, ncells) dm = mesh.topology_dm - # Remove default Face Sets so we have full control over labels. - if dm.hasLabel("Face Sets"): - dm.removeLabel("Face Sets") + # Remove default labels so we have full control. + for name in ("Face Sets", "Edge Sets", "Vertex Sets"): + if dm.hasLabel(name): + dm.removeLabel(name) dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) eStart, eEnd = dm.getDepthStratum(1) @@ -92,9 +93,10 @@ def _make_square_with_vertex_sets(ncells=2): """ mesh = UnitSquareMesh(ncells, ncells) dm = mesh.topology_dm - # Remove default Face Sets so we have full control over labels. - if dm.hasLabel("Face Sets"): - dm.removeLabel("Face Sets") + # Remove default labels so we have full control. + for name in ("Face Sets", "Vertex Sets"): + if dm.hasLabel(name): + dm.removeLabel(name) dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) vStart, vEnd = dm.getDepthStratum(0) @@ -203,8 +205,9 @@ def test_submesh_codim1_no_parent_edge_sets(): """Codim-1 submesh with no Edge Sets on the parent still labels exterior facets.""" mesh = UnitCubeMesh(2, 2, 2) dm = mesh.topology_dm - if dm.hasLabel("Face Sets"): - dm.removeLabel("Face Sets") + for name in ("Face Sets", "Edge Sets", "Vertex Sets"): + if dm.hasLabel(name): + dm.removeLabel(name) dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) fStart, fEnd = dm.getDepthStratum(2) From 055e5779c0cf44a5a776630e282474ba26047a8d Mon Sep 17 00:00:00 2001 From: lorenzoCorintis Date: Wed, 11 Mar 2026 18:18:58 +0100 Subject: [PATCH 15/15] Rebuild submesh "Face Sets" from lower-dimensional labels --- firedrake/cython/dmcommon.pyx | 262 +++++++++++------- firedrake/cython/petschdr.pxi | 2 + .../submesh/test_submesh_codim1_labels.py | 19 +- 3 files changed, 171 insertions(+), 112 deletions(-) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index ade1fb63ae..26f8666808 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -11,7 +11,7 @@ from firedrake.petsc import PETSc from mpi4py import MPI from firedrake.utils import IntType, ScalarType from libc.string cimport memset -from libc.stdlib cimport qsort +from libc.stdlib cimport qsort, malloc, free from finat.element_factory import as_fiat_cell cimport numpy as np @@ -1726,75 +1726,18 @@ def label_facets(PETSc.DM plex): CHKERR(DMLabelDestroyIndex(lbl_ext)) -cdef void _populate_lower_dim_labels(PETSc.DM dm, PetscInt tdim): - """Copy "Face Sets" values to "Edge Sets" and "Vertex Sets" for - lower-dimensional entities. - - After ``DMPlexLabelComplete`` has propagated "Face Sets" into the - closure, edges and vertices carry the boundary tag of their parent - face. This helper mirrors those values into the corresponding - lower-dimensional label so that a codimension-1 ``Submesh`` can - find the information it needs without relying on the inherited - "Face Sets" (which mixes cell- and facet-level values). - - Entities that already have a value in the target label are skipped - to avoid overwriting user-defined labels (e.g. from Gmsh). - """ - cdef: - PetscInt pStart, pEnd, p, face_val, existing_val - DMLabel face_sets_lbl, edge_sets_lbl, vertex_sets_lbl - - CHKERR(DMGetLabel(dm.dm, b"Face Sets", &face_sets_lbl)) - - if tdim >= 3: - if not dm.hasLabel(EDGE_SETS_LABEL): - dm.createLabel(EDGE_SETS_LABEL) - CHKERR(DMGetLabel(dm.dm, b"Edge Sets", &edge_sets_lbl)) - CHKERR(DMPlexGetDepthStratum(dm.dm, 1, &pStart, &pEnd)) - for p in range(pStart, pEnd): - CHKERR(DMLabelGetValue(face_sets_lbl, p, &face_val)) - if face_val >= 0: - CHKERR(DMLabelGetValue(edge_sets_lbl, p, &existing_val)) - if existing_val < 0: - CHKERR(DMLabelSetValue(edge_sets_lbl, p, face_val)) - - if tdim >= 2: - if not dm.hasLabel(VERTEX_SETS_LABEL): - dm.createLabel(VERTEX_SETS_LABEL) - CHKERR(DMGetLabel(dm.dm, b"Vertex Sets", &vertex_sets_lbl)) - CHKERR(DMPlexGetDepthStratum(dm.dm, 0, &pStart, &pEnd)) - for p in range(pStart, pEnd): - CHKERR(DMLabelGetValue(face_sets_lbl, p, &face_val)) - if face_val >= 0: - CHKERR(DMLabelGetValue(vertex_sets_lbl, p, &existing_val)) - if existing_val < 0: - CHKERR(DMLabelSetValue(vertex_sets_lbl, p, face_val)) - - def complete_facet_labels(PETSc.DM dm): """Transfer label values from the facet labels to everything in - the closure of the facets. + the closure of the facets.""" + cdef PETSc.DMLabel label - After completing "Face Sets" this also populates "Edge Sets" and - "Vertex Sets" so that boundary information is available at every - topological dimension. Entities that already carry a value in the - target label (e.g. user-defined labels from Gmsh) are not overwritten. - """ - cdef: - PETSc.DMLabel label - PetscInt tdim - - tdim = get_topological_dimension(dm) - if tdim == 0: + if get_topological_dimension(dm) == 0: return for name in [FACE_SETS_LABEL, "exterior_facets", "interior_facets"]: if dm.hasLabel(name): label = dm.getLabel(name) CHKERR( DMPlexLabelComplete(dm.dm, label.dmlabel) ) - if dm.hasLabel(FACE_SETS_LABEL): - _populate_lower_dim_labels(dm, tdim) - @cython.boundscheck(False) @cython.wraparound(False) @@ -3886,6 +3829,8 @@ def submesh_create(PETSc.DM dm, CHKERR(DMLabelSetValue(temp_label.dmlabel, p, label_value)) CHKERR(ISRestoreIndices(stratum_is, &stratum_indices)) CHKERR(ISDestroy(&stratum_is)) + if subdim == dm.getDimension() - 1 and dm.hasLabel(FACE_SETS_LABEL): + _populate_lower_dim_labels(dm, dm.getDimension()) # Make submesh using temp_label. subdm, ownership_transfer_sf = dm.filter(label=temp_label, value=label_value, @@ -3996,6 +3941,89 @@ def submesh_correct_entity_classes(PETSc.DM dm, CHKERR(DMLabelDestroyIndex(lbl_ghost)) +cdef void _populate_lower_dim_labels(PETSc.DM dm, PetscInt tdim): + """Derive "Edge Sets" / "Vertex Sets" from positive "Face Sets" strata. + + Iterates each "Face Sets" stratum value > 0 and copies the value to + the corresponding lower-dimensional label for depth-1 or depth-0 + entities. Entities that already carry a positive value in the + target label are skipped to preserve user-defined labels (e.g. from + Gmsh). Stale value-0 entries (mesh-generator artifacts) are + cleared before the positive value is set, so that + ``DMLabelGetValue`` returns the meaningful value rather than 0. + """ + cdef: + PetscInt pStart, pEnd, p, nvals, i, j, stratum_val, stratum_size + PetscInt existing_val + DMLabel face_sets_lbl, target_lbl + PETSc.PetscIS stratum_is = NULL + const PetscInt *stratum_pts = NULL + + CHKERR(DMGetLabel(dm.dm, b"Face Sets", &face_sets_lbl)) + CHKERR(DMLabelGetNumValues(face_sets_lbl, &nvals)) + if nvals == 0: + return + + cdef PetscInt *vals = malloc(nvals * sizeof(PetscInt)) + with dm.getLabelIdIS(FACE_SETS_LABEL) as ids: + for i in range(nvals): + vals[i] = ids[i] + + if tdim >= 3: + if not dm.hasLabel(EDGE_SETS_LABEL): + dm.createLabel(EDGE_SETS_LABEL) + CHKERR(DMGetLabel(dm.dm, b"Edge Sets", &target_lbl)) + CHKERR(DMPlexGetDepthStratum(dm.dm, 1, &pStart, &pEnd)) + for i in range(nvals): + stratum_val = vals[i] + if stratum_val <= 0: + continue + CHKERR(DMLabelGetStratumSize(face_sets_lbl, stratum_val, &stratum_size)) + if stratum_size == 0: + continue + CHKERR(DMLabelGetStratumIS(face_sets_lbl, stratum_val, &stratum_is)) + CHKERR(ISGetIndices(stratum_is, &stratum_pts)) + for j in range(stratum_size): + p = stratum_pts[j] + if pStart <= p < pEnd: + CHKERR(DMLabelGetValue(target_lbl, p, &existing_val)) + if existing_val > 0: + continue + if existing_val == 0: + CHKERR(DMLabelClearValue(target_lbl, p, 0)) + CHKERR(DMLabelSetValue(target_lbl, p, stratum_val)) + CHKERR(ISRestoreIndices(stratum_is, &stratum_pts)) + CHKERR(ISDestroy(&stratum_is)) + + if tdim >= 2: + if not dm.hasLabel(VERTEX_SETS_LABEL): + dm.createLabel(VERTEX_SETS_LABEL) + CHKERR(DMGetLabel(dm.dm, b"Vertex Sets", &target_lbl)) + CHKERR(DMPlexGetDepthStratum(dm.dm, 0, &pStart, &pEnd)) + for i in range(nvals): + stratum_val = vals[i] + if stratum_val <= 0: + continue + CHKERR(DMLabelGetStratumSize(face_sets_lbl, stratum_val, &stratum_size)) + if stratum_size == 0: + continue + CHKERR(DMLabelGetStratumIS(face_sets_lbl, stratum_val, &stratum_is)) + CHKERR(ISGetIndices(stratum_is, &stratum_pts)) + for j in range(stratum_size): + p = stratum_pts[j] + if pStart <= p < pEnd: + CHKERR(DMLabelGetValue(target_lbl, p, &existing_val)) + if existing_val > 0: + continue + if existing_val == 0: + CHKERR(DMLabelClearValue(target_lbl, p, 0)) + CHKERR(DMLabelSetValue(target_lbl, p, stratum_val)) + CHKERR(ISRestoreIndices(stratum_is, &stratum_pts)) + CHKERR(ISDestroy(&stratum_is)) + + free(vals) + + cdef PetscInt _max_label_value(PETSc.DM dm, label_name): """Return the maximum value in *label_name*, or -1 if absent/empty.""" if not dm.hasLabel(label_name): @@ -4039,29 +4067,34 @@ cdef void _label_new_exterior_facets( cdef void _propagate_parent_facet_labels( - PETSc.DM subdm, + PETSc.DM dm, PETSc.DM subdm, PetscInt subdim, const PetscInt *sub_ext_facet_indices, PetscInt sub_ext_facet_size, PetscInt subfStart, PetscInt subfEnd, ): - """Codimension-1 helper: populate "Face Sets" from the subdm's own - lower-dimensional labels. - - ``DMPlexFilter`` copies every label from the parent into the subdm, - including "Edge Sets" and "Vertex Sets". Inherited "Face Sets" - values are preserved because removing and recreating the label - loses PETSc's internal parallel label-migration state, breaking - chained submesh operations (e.g. 3D → 2D → 1D). - - Only exterior facets that do *not* already have an inherited - "Face Sets" value are updated: they receive the subdm's own - "Edge Sets" (3D→2D) or "Vertex Sets" (2D→1D) value, or a fresh - default tag if none exists. + """Codimension-1 helper: rebuild "Face Sets" purely from "Edge Sets" + (3D→2D) or "Vertex Sets" (2D→1D). + + ``DMPlexFilter`` copies every label from the parent into the subdm. + The inherited "Face Sets" mixes cell- and facet-level values from + the parent, so we clear every stratum with ``DMLabelClearStratum`` + (preserving the label object) and rebuild only from the subdm's own + lower-dimensional labels. ``DMPlexLabelComplete`` is called + afterwards to propagate values to closure entities and synchronise + ghost points via the point SF. + + When no lower-dimensional source label exists, inherited "Face Sets" + values are preserved and only unlabeled exterior facets receive a + fresh default tag. """ cdef: - PetscInt pStart, pEnd, next_label_val, label_val, existing_val, subf, i + PetscInt pStart, pEnd, next_label_val, label_val, existing_val + PetscInt subf, i, nvals, local_has DMLabel source_label, face_sets_label + PetscInt *stratum_vals_arr = NULL + PetscBool has_point + PETSc.DMLabel face_sets_py if subdim == 2: source_label_name = b"Edge Sets" @@ -4070,40 +4103,67 @@ cdef void _propagate_parent_facet_labels( else: source_label_name = None - has_source = (source_label_name is not None - and subdm.hasLabel(source_label_name)) + # Determine has_source consistently across all ranks. + if source_label_name is not None and subdm.hasLabel(source_label_name): + local_has = 1 if _max_label_value(subdm, source_label_name) > 0 else 0 + else: + local_has = 0 + has_source = bool(dm.comm.tompi4py().allreduce(local_has, op=MPI.MAX)) next_label_val = max( - _max_label_value(subdm, source_label_name) if has_source else -1, + _max_label_value(dm, source_label_name) + if (source_label_name is not None and dm.hasLabel(source_label_name)) + else -1, _max_label_value(subdm, FACE_SETS_LABEL), ) + 1 - next_label_val = subdm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) + next_label_val = dm.comm.tompi4py().allreduce(next_label_val, op=MPI.MAX) if not subdm.hasLabel(FACE_SETS_LABEL): subdm.createLabel(FACE_SETS_LABEL) CHKERR(DMGetLabel(subdm.dm, b"Face Sets", &face_sets_label)) if has_source: - CHKERR(DMGetLabel(subdm.dm, source_label_name, &source_label)) - pStart, pEnd = subdm.getChart() - CHKERR(DMLabelCreateIndex(source_label, pStart, pEnd)) - - for i in range(sub_ext_facet_size): - subf = sub_ext_facet_indices[i] - if subfStart <= subf < subfEnd: - CHKERR(DMLabelGetValue(face_sets_label, subf, &existing_val)) - if existing_val >= 0: - continue - label_val = -1 - if has_source: - CHKERR(DMLabelGetValue(source_label, subf, &label_val)) - if label_val >= 0: - CHKERR(DMLabelSetValue(face_sets_label, subf, label_val)) - else: - CHKERR(DMSetLabelValue(subdm.dm, b"Face Sets", subf, next_label_val)) - - if has_source: - CHKERR(DMLabelDestroyIndex(source_label)) + # --- Clear + rebuild path --- + CHKERR(DMLabelGetNumValues(face_sets_label, &nvals)) + if nvals > 0: + stratum_vals_arr = malloc(nvals * sizeof(PetscInt)) + with subdm.getLabelIdIS(FACE_SETS_LABEL) as ids: + for i in range(nvals): + stratum_vals_arr[i] = ids[i] + for i in range(nvals): + CHKERR(DMLabelClearStratum(face_sets_label, stratum_vals_arr[i])) + free(stratum_vals_arr) + + if subdm.hasLabel(source_label_name): + CHKERR(DMGetLabel(subdm.dm, source_label_name, + &source_label)) + pStart, pEnd = subdm.getChart() + CHKERR(DMLabelCreateIndex(source_label, pStart, pEnd)) + for subf in range(subfStart, subfEnd): + CHKERR(DMLabelHasPoint(source_label, subf, &has_point)) + if has_point: + CHKERR(DMLabelGetValue(source_label, subf, &label_val)) + if label_val > 0: + CHKERR(DMLabelSetValue(face_sets_label, subf, label_val)) + CHKERR(DMLabelDestroyIndex(source_label)) + + for i in range(sub_ext_facet_size): + subf = sub_ext_facet_indices[i] + if subfStart <= subf < subfEnd: + CHKERR(DMLabelGetValue(face_sets_label, subf, &label_val)) + if label_val < 0: + CHKERR(DMLabelSetValue(face_sets_label, subf, next_label_val)) + + face_sets_py = subdm.getLabel(FACE_SETS_LABEL) + CHKERR(DMPlexLabelComplete(subdm.dm, face_sets_py.dmlabel)) + else: + # --- Preserve inherited path --- + for i in range(sub_ext_facet_size): + subf = sub_ext_facet_indices[i] + if subfStart <= subf < subfEnd: + CHKERR(DMLabelGetValue(face_sets_label, subf, &existing_val)) + if existing_val < 0: + CHKERR(DMLabelSetValue(face_sets_label, subf, next_label_val)) @cython.boundscheck(False) @@ -4157,7 +4217,7 @@ def submesh_update_facet_labels(PETSc.DM dm, PETSc.DM subdm): CHKERR(ISRestoreIndices(subpoint_is.iset, &subpoint_indices)) elif subdim == dim - 1: _propagate_parent_facet_labels( - subdm, subdim, + dm, subdm, subdim, sub_ext_facet_indices, sub_ext_facet_size, subfStart, subfEnd) diff --git a/firedrake/cython/petschdr.pxi b/firedrake/cython/petschdr.pxi index d7b1478ca0..bca0dfee7a 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 DMLabelClearStratum(DMLabel, PetscInt) + PetscErrorCode DMLabelGetNumValues(DMLabel, PetscInt*) cdef extern from "petscdm.h" nogil: PetscErrorCode DMCreateLabel(PETSc.PetscDM,char[]) diff --git a/tests/firedrake/submesh/test_submesh_codim1_labels.py b/tests/firedrake/submesh/test_submesh_codim1_labels.py index 41dced0ee2..f2640d9808 100644 --- a/tests/firedrake/submesh/test_submesh_codim1_labels.py +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -49,10 +49,9 @@ def _make_cube_with_edge_sets(extra_edge_c=False, ncells=2): """ mesh = UnitCubeMesh(ncells, ncells, ncells) dm = mesh.topology_dm - # Remove default labels so we have full control. - for name in ("Face Sets", "Edge Sets", "Vertex Sets"): - if dm.hasLabel(name): - dm.removeLabel(name) + # Remove default Face Sets so we have full control over labels. + if dm.hasLabel("Face Sets"): + dm.removeLabel("Face Sets") dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) eStart, eEnd = dm.getDepthStratum(1) @@ -93,10 +92,9 @@ def _make_square_with_vertex_sets(ncells=2): """ mesh = UnitSquareMesh(ncells, ncells) dm = mesh.topology_dm - # Remove default labels so we have full control. - for name in ("Face Sets", "Vertex Sets"): - if dm.hasLabel(name): - dm.removeLabel(name) + # Remove default Face Sets so we have full control over labels. + if dm.hasLabel("Face Sets"): + dm.removeLabel("Face Sets") dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) vStart, vEnd = dm.getDepthStratum(0) @@ -205,9 +203,8 @@ def test_submesh_codim1_no_parent_edge_sets(): """Codim-1 submesh with no Edge Sets on the parent still labels exterior facets.""" mesh = UnitCubeMesh(2, 2, 2) dm = mesh.topology_dm - for name in ("Face Sets", "Edge Sets", "Vertex Sets"): - if dm.hasLabel(name): - dm.removeLabel(name) + if dm.hasLabel("Face Sets"): + dm.removeLabel("Face Sets") dm.createLabel("Face Sets") vc = _get_dm_vertex_coords(dm) fStart, fEnd = dm.getDepthStratum(2)