Fix Submesh to propagate "Edge Sets" / "Vertex Sets" to codimension-1 submesh exterior facet markers

[lorenzoCorintis]

Description

• Fixes BUG: Submesh does not propagate "Edge Sets" to submesh exterior facet markers #4948 -- Submesh now propagates parent "Edge Sets" (3D→2D) and "Vertex Sets" (2D→1D) labels to "Face Sets" on codimension-1 submeshes, so that submesh.exterior_facets.unique_markers contains the correct boundary markers and ds(tag) integration works as expected.
• The previous implementation of submesh_update_facet_labels returned early for subdim != dim, leaving codimension-1 submeshes with no meaningful facet labels. The new else branch maps each submesh exterior facet back to its parent point via subpoint_is, looks up the parent label value, and writes it into "Face Sets" on the subdm. Exterior facets without a parent label receive a fresh default value (max(existing_values) + 1).
• Inherited "Face Sets" from DMPlexFilter (which sit on cells, not facets) are removed before relabeling to prevent contamination of exterior facet markers.

[connorjward]

I think that this is conceptually correct but the code as written is very very complicated. It would really benefit from a refactor to make it more readable.

[connorjward]

I would like to see more parallel tests, especially some that involve calculation, as that is where this functionality is really stressed.

[lorenzoCorintis]

I added parallel computation tests at 2 and 3 ranks: boundary length integration via ds(tag), coordinate expression integrals over tagged boundaries, CG2 interpolation + boundary integration, and 2D→1D vertex set propagation with ds(tag).
Let me know if there are more specific tests you would like to see

[lorenzoCorintis]

I think that this is conceptually correct but the code as written is very very complicated. It would really benefit from a refactor to make it more readable.

Thanks for the feedback! I tried to make it more readable, but please let me know if there are any other changes you would like me to do.

[connorjward]

In general I think this is really good.

[connorjward]

I don't understand this. Can you explain? I would expect that the given "Face Sets" is now just wrong and it seems a bit dodgy to be partially editing the "Face Sets" label.

[lorenzoCorintis]

DMPlexFilter copies every label from the parent into the subdm, including "Face Sets" values that DMPlexLabelComplete propagated to closure entities (edges, vertices). For a codim-1 submesh these inherited edge-level values carry useful boundary information.

I tried discarding the inherited "Face Sets" and rebuilding from scratch, but this breaks test_submesh_solve_2d_1d_poisson_hermite (which chains 3D→2D→1D submesh extractions): removing and recreating the label loses PETSc's internal parallel migration state, so the rebuilt label doesn't survive redistribution correctly.

I therefore preserve inherited values and only fill in exterior facets that don't already have one, using the parent's lower-dimensional label or a fresh default. The docstring now explains this.

[connorjward]

This seems quite counter intuitive and might cause some confusion down the line. I wonder if there's a simpler solution:

When we propagate "Face Sets" to the rest of the facet closure could we at the same time set "Edge Sets" and "Vertex Sets" for the extra entities.

Consider one of our 3D utility meshes with some marked boundary. "Face Sets" for this boundary would contain all facets, edges and vertices on the boundary, "Edge Sets" would contain all edges and vertices, and "Vertex Sets" just the vertices.

This way when we lose a dimension via submesh all the necessary boundary information is already encoded in "Edge Sets" and so "Face Sets" is safe to discard (or push to "Cell Sets"?).

[lorenzoCorintis]

I've implemented part 1 of your suggestion: complete_facet_labels now populates "Edge Sets" and "Vertex Sets" from the "Face Sets" closure (entities that already carry a value, e.g. from Gmsh, are not overwritten). So after mesh setup a 3D utility mesh with a marked boundary will have:

• "Face Sets": all facets, edges and vertices on the boundary
• "Edge Sets": all edges and vertices
• "Vertex Sets": just the vertices

For part 2 I also switched _propagate_parent_facet_labels to read from the subdm's own inherited "Edge Sets" / "Vertex Sets" instead of doing cross-DM queries against the parent via subpoint_indices. This simplifies the function signature and removes the parent-DM dependency entirely.

However, I was not able to fully discard inherited "Face Sets" on the subdm (the removeLabel / createLabel approach). Even when reading exclusively from subdm-local labels, removing and recreating "Face Sets" still breaks test_submesh_solve_2d_1d_poisson_hermite (which chains 3D→2D→1D). The issue is that removeLabel / createLabel loses PETSc's internal parallel label-migration state regardless of where the replacement values come from.

So the current implementation preserves inherited "Face Sets" and only fills in exterior facets that don't already have a value. I think fully discarding inherited "Face Sets" would require PETSc-level changes (e.g. a way to reset label values without destroying the migration SF), which is out of scope here.

[connorjward]

The issue that I am having is that there is some quite complicated implicit union behaviour going on between the closure-propagated "Face Sets" and any new "Edge Sets" entries. For instance consider having a cube:

           +-------+   
          /       /|     
         /       / |     
      1 /       /  |  
       /       /   |     
      /       /    |     
     +-------+  1  +
     |       |    /   
     |       |   /    
   1 |       |  / 
     |       | /
     |       |/
     +-------+

where the left-most edges are labelled 1 in "Edge Sets" and the right-most face is labelled 1 in "Face Sets". If you take some sort of codim-1 submesh it is not obvious which edges (now faces) should be labelled 1.

I think my preferred approach is:

1. Only "Edge Sets" are used to build "Face Sets" for the submesh
2. Build "Edge Sets" and "Vertex Sets" only for our utility meshes, instead of always propagating from "Face Sets" - this is optional for this PR since my guess is you are mostly using gmsh.

Sorry for making this such a big deal, but I think we have to be careful about this.

However, I was not able to fully discard inherited "Face Sets" on the subdm (the removeLabel / createLabel approach).

What about something like DMLabelClearStratum?

Even when reading exclusively from subdm-local labels, removing and recreating "Face Sets" still breaks test_submesh_solve_2d_1d_poisson_hermite

This is weird. I wonder what is happening there.

[lorenzoCorintis]

@connorjward

Thanks for the suggestion — we've now implemented the approach you described:
1. Populate "Edge Sets" and "Vertex Sets" from "Face Sets" on the parent
Before DMPlexFilter, a new helper _populate_lower_dim_labels iterates each positive "Face Sets" stratum on the parent mesh and copies the value to "Edge Sets" (depth-1 entities) and "Vertex Sets" (depth-0 entities). Entities that already carry a positive value (e.g. user-defined labels from Gmsh) are preserved.
2. Clear inherited "Face Sets" on the subdm and rebuild from "Edge Sets"
_propagate_parent_facet_labels now uses DMLabelClearStratum to clear every stratum of the inherited "Face Sets" on the subdm, then rebuilds it purely from the subdm's own "Edge Sets" (3D→2D) or "Vertex Sets" (2D→1D). DMPlexLabelComplete is called afterwards to propagate to closure entities and synchronise ghost points via the point SF. When no source label exists, a fallback path preserves inherited values.
3. Handling mesh-generator artifacts (value 0)
The main difficulty we hit was that mesh files (and RelabeledMesh) can leave "Face Sets" value 0 on many boundary entities. Since DMLabelGetValue returns the minimum value for a point, edges at face intersections appeared to have value 0 rather than their meaningful positive value. The fix was to:

• Iterate per-stratum (via DMLabelGetStratumIS) rather than relying on DMLabelGetValue
• Clear stale value-0 entries with DMLabelClearValue before setting the correct positive value
• Skip stratum 0 entirely during derivation
  All existing submesh tests pass, including test_submesh_solve_2d_1d_poisson_hermite (nprocs=7) which chains 3D→2D→1D submesh extractions via RelabeledMesh.

[connorjward]

Just make this the test and stick the pytest.mark.parallel here.

[connorjward]

There are a couple more instances of this elsewhere.

[lorenzoCorintis]

Ah yes! Done

[connorjward]

@lorenzoCorintis would you say that this PR closes #4835?

[lorenzoCorintis]

@lorenzoCorintis would you say that this PR closes #4835?

@connorjward The issue in #4835 is different from what this PR fixes.

What #4835 reports: when creating a codim-1 submesh via RelabeledMesh + Submesh, the submesh's exterior_facets.unique_markers shows [1, 2, 3, 4] (the original mesh's markers) instead of just [1] (the relabeled marker).

Why this PR doesn't fix it: the root cause is that RelabeledMesh only clears the specific stratum being set (clearLabelStratum("Face Sets", 1)) but leaves the original strata 2, 3, 4 intact on the cloned plex. So the relabeled mesh still carries "Face Sets" = {1, 2, 3, 4}, DMPlexFilter copies all of those into the submesh, and they persist. That is a RelabeledMesh issue, not a submesh_update_facet_labels issue.

What this PR does fix: propagation of lower-dimensional labels ("Edge Sets" → "Face Sets" for 3D→2D, "Vertex Sets" → "Face Sets" for 2D→1D) on codimension-1 submeshes, which previously wasn't handled at all.

[ksagiyam]

I only skimmed at the description, but it sounds more like a PETSc issue as PETSc basically reserves "{Cell, Face, Edge, Vertex} Sets" labels. I think the actual issue is simply that these label names are not correctly mapped when making codim-1 submeshes in PETSc. The following simple naive patch almost fixes "Face Sets" (of parent) -> "Cell Sets" (of submesh) case:

diff --git a/src/dm/impls/plex/plexsubmesh.c b/src/dm/impls/plex/plexsubmesh.c
index c21245e0b9a..85fbb49989f 100644
--- a/src/dm/impls/plex/plexsubmesh.c
+++ b/src/dm/impls/plex/plexsubmesh.c
@@ -3336,9 +3336,24 @@ static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPo
     PetscCall(PetscStrcmp(lname, "celltype", &isCelltype));
     PetscCall(PetscStrcmp(lname, "vtk", &isVTK));
     if (isDepth || isDim || isCelltype || isVTK) continue;
-    PetscCall(DMCreateLabel(subdm, lname));
-    PetscCall(DMGetLabel(dm, lname, &label));
-    PetscCall(DMGetLabel(subdm, lname, &newlabel));
+    {
+      PetscInt  dim, subdim;
+      PetscBool flag;
+      char     *sublname;
+
+      PetscCall(DMGetDimension(dm, &dim));
+      PetscCall(DMGetDimension(subdm, &subdim));
+      PetscCall(PetscStrcmp(lname, "Face Sets", &flag));
+      if (subdim == dim - 1 && flag) {
+        PetscCall(PetscStrallocpy("Cell Sets", &sublname));
+      } else {
+        PetscCall(PetscStrallocpy(lname, &sublname));
+      }
+      PetscCall(DMCreateLabel(subdm, sublname));
+      PetscCall(DMGetLabel(dm, lname, &label));
+      PetscCall(DMGetLabel(subdm, sublname, &newlabel));
+      PetscCall(PetscFree(sublname));
+    }
     PetscCall(DMLabelGetDefaultValue(label, &v));
     PetscCall(DMLabelSetDefaultValue(newlabel, v));
     PetscCall(DMLabelGetValueIS(label, &valueIS));

from firedrake import *


mesh = UnitSquareMesh(1, 1, quadrilateral=True)
subm = Submesh(mesh, mesh.topological_dimension - 1, 4)
subm.topology_dm.viewFromOptions("-dm_view")

DM Object: firedrake_default_submesh_topology 1 MPI process
  type: plex
firedrake_default_submesh_topology in 1 dimension:
Supports:
[0] Max support size: 1
[0]: 1 ----> 0
[0]: 2 ----> 0
Cones:
[0] Max cone size: 2
[0]: 0 <---- 1 (0)
[0]: 0 <---- 2 (0)
coordinates with 1 fields
  field 0 with 2 components
Process 0:
  (   1) dof  2 offset   0 0. 1.
  (   2) dof  2 offset   2 1. 1.
Labels:
Label 'celltype':
[0]: 1 (0)
[0]: 2 (0)
[0]: 0 (1)
Label 'Cell Sets':
[0]: 1 (1)
[0]: 0 (4)
[0]: 1 (4)
[0]: 2 (4)
[0]: 2 (2)
Label 'pyop2_core':
[0]: 0 (1)
[0]: 1 (1)
[0]: 2 (1)
Label 'pyop2_owned':
Label 'pyop2_ghost':
Label 'exterior_facets':
[0]: 1 (1)
[0]: 2 (1)
Label 'interior_facets':

You can remove redundant labels/points (label 1 and 2 in "Cell Sets" in the above) somewhere.

At least, it is worth speaking to Matt about this, I think

[lorenzoCorintis]

I only skimmed at the description, but it sounds more like a PETSc issue as PETSc basically reserves "{Cell, Face, Edge, Vertex} Sets" labels. I think the actual issue is simply that these label names are not correctly mapped when making codim-1 submeshes in PETSc. The following simple naive patch almost fixes "Face Sets" (of parent) -> "Cell Sets" (of submesh) case:

diff --git a/src/dm/impls/plex/plexsubmesh.c b/src/dm/impls/plex/plexsubmesh.c
index c21245e0b9a..85fbb49989f 100644
--- a/src/dm/impls/plex/plexsubmesh.c
+++ b/src/dm/impls/plex/plexsubmesh.c
@@ -3336,9 +3336,24 @@ static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPo
     PetscCall(PetscStrcmp(lname, "celltype", &isCelltype));
     PetscCall(PetscStrcmp(lname, "vtk", &isVTK));
     if (isDepth || isDim || isCelltype || isVTK) continue;
-    PetscCall(DMCreateLabel(subdm, lname));
-    PetscCall(DMGetLabel(dm, lname, &label));
-    PetscCall(DMGetLabel(subdm, lname, &newlabel));
+    {
+      PetscInt  dim, subdim;
+      PetscBool flag;
+      char     *sublname;
+
+      PetscCall(DMGetDimension(dm, &dim));
+      PetscCall(DMGetDimension(subdm, &subdim));
+      PetscCall(PetscStrcmp(lname, "Face Sets", &flag));
+      if (subdim == dim - 1 && flag) {
+        PetscCall(PetscStrallocpy("Cell Sets", &sublname));
+      } else {
+        PetscCall(PetscStrallocpy(lname, &sublname));
+      }
+      PetscCall(DMCreateLabel(subdm, sublname));
+      PetscCall(DMGetLabel(dm, lname, &label));
+      PetscCall(DMGetLabel(subdm, sublname, &newlabel));
+      PetscCall(PetscFree(sublname));
+    }
     PetscCall(DMLabelGetDefaultValue(label, &v));
     PetscCall(DMLabelSetDefaultValue(newlabel, v));
     PetscCall(DMLabelGetValueIS(label, &valueIS));

from firedrake import *


mesh = UnitSquareMesh(1, 1, quadrilateral=True)
subm = Submesh(mesh, mesh.topological_dimension - 1, 4)
subm.topology_dm.viewFromOptions("-dm_view")

DM Object: firedrake_default_submesh_topology 1 MPI process
  type: plex
firedrake_default_submesh_topology in 1 dimension:
Supports:
[0] Max support size: 1
[0]: 1 ----> 0
[0]: 2 ----> 0
Cones:
[0] Max cone size: 2
[0]: 0 <---- 1 (0)
[0]: 0 <---- 2 (0)
coordinates with 1 fields
  field 0 with 2 components
Process 0:
  (   1) dof  2 offset   0 0. 1.
  (   2) dof  2 offset   2 1. 1.
Labels:
Label 'celltype':
[0]: 1 (0)
[0]: 2 (0)
[0]: 0 (1)
Label 'Cell Sets':
[0]: 1 (1)
[0]: 0 (4)
[0]: 1 (4)
[0]: 2 (4)
[0]: 2 (2)
Label 'pyop2_core':
[0]: 0 (1)
[0]: 1 (1)
[0]: 2 (1)
Label 'pyop2_owned':
Label 'pyop2_ghost':
Label 'exterior_facets':
[0]: 1 (1)
[0]: 2 (1)
Label 'interior_facets':

You can remove redundant labels/points (label 1 and 2 in "Cell Sets" in the above) somewhere.

At least, it is worth speaking to Matt about this, I think

That's a great point — a PETSc-level label remapping in DMPlexFilterLabels_Internal would be the cleanest solution.

Our approach here is effectively doing the same thing from the Firedrake side: we derive "Edge Sets" / "Vertex Sets" from "Face Sets" on the parent (via _populate_lower_dim_labels), then after DMPlexFilter we rebuild the subdm's "Face Sets" purely from the inherited "Edge Sets" (3D→2D) or "Vertex Sets" (2D→1D).

A PETSc patch that shifts the label names during DMPlexFilter (Face Sets → Cell Sets, Edge Sets → Face Sets, Vertex Sets → Edge Sets) would make our Firedrake-side code unnecessary. It would probably also need to handle:

• The full chain of shifts (not just Face → Cell)
• Stripping closure entities that no longer belong to the target depth after the dimension shift
• The complete_facet_labels call that Firedrake does after mesh creation

Agreed that it's worth discussing it further. In the meantime, this PR provides a working Firedrake-level solution — if PETSc gains native support we can simplify later.

[connorjward]

In my feedback I suggested that this only be done for our utility meshes (i.e. not gmsh). This still has the confusing intersection behaviour.

[connorjward]

Maybe this code isn't even necessary because if "Edge Sets" isn't defined then we can use bits from "Face Sets". The thing is we have to avoid blending the two which means if "Edge Sets" is defined we shouldn't be using the original "Face Sets" to build the "Face Sets" for the sub dm.

Are there cases where this doesn't work? If this is critical functionality we could consider relabelling "Edge Sets" with some offset to preserve uniqueness. WDYT?

[connorjward]

We prefer PetscMalloc and PetscFree