Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lorenzoCorintis wants to merge 17 commits into
base: main
Choose a base branch
from
Open

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

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
afde9fa
fix submesh labels
lorenzoCorintis Mar 9, 2026
19d97e8
add tests
lorenzoCorintis Mar 9, 2026
1ac8e61
make an else condition explicit
lorenzoCorintis Mar 10, 2026
19b188f
add some tests
lorenzoCorintis Mar 10, 2026
fcbca35
fix tests
lorenzoCorintis Mar 10, 2026
88446cd
small refactoring
lorenzoCorintis Mar 10, 2026
eec6bed
Update tests/firedrake/submesh/test_submesh_codim1_labels.py
lorenzoCorintis Mar 10, 2026
2cfffed
Update tests/firedrake/submesh/test_submesh_codim1_labels.py
lorenzoCorintis Mar 10, 2026
7507425
Apply suggestions from code review
lorenzoCorintis Mar 10, 2026
307c688
lint fixes
lorenzoCorintis Mar 10, 2026
6f9c0a8
Merge branch 'fix-submesh-labels' of github.com:lorenzoCorintis/fired…
lorenzoCorintis Mar 10, 2026
b2bdc65
Update firedrake/cython/dmcommon.pyx
lorenzoCorintis Mar 10, 2026
5629346
improve a docstring
lorenzoCorintis Mar 11, 2026
b4bd689
Merge branch 'fix-submesh-labels' of github.com:lorenzoCorintis/fired…
lorenzoCorintis Mar 11, 2026
e3f0b85
refactor tests
lorenzoCorintis Mar 11, 2026
5f30402
connor's suggestion
lorenzoCorintis Mar 11, 2026
055e577
Rebuild submesh "Face Sets" from lower-dimensional labels
lorenzoCorintis Mar 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
184 changes: 138 additions & 46 deletions firedrake/cython/dmcommon.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3937,76 +3937,168 @@ 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,
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

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)

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 that already carry
an inherited "Face Sets" value (from ``DMPlexFilter``) are left alone;
Copy link
Copy Markdown
Contributor

@connorjward connorjward Mar 10, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Author

@lorenzoCorintis lorenzoCorintis Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Contributor

@connorjward connorjward Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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"?).

Copy link
Copy Markdown
Author

@lorenzoCorintis lorenzoCorintis Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Contributor

@connorjward connorjward Mar 11, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Author

@lorenzoCorintis lorenzoCorintis Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

the remaining ones receive the parent label value or a default tag.
"""
cdef:
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"
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 = 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
# (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, <const char *>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:
# 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:
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
-----
This function marks the new exterior facets with current max label value + 1 in "Face Sets".
* **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.
Comment thread
lorenzoCorintis marked this conversation as resolved.
Outdated

"""
cdef:
PetscInt dim, subdim, pStart, pEnd, f, subfStart, subfEnd, subf, sub_ext_facet_size, next_label_val, i
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 = <char *>"interior_facets"
char *ext_facet_label_name = <char *>"exterior_facets"
char *face_sets_label_name = <char *>"Face Sets"
DMLabel ext_facet_label
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))
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:
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))
else:
raise NotImplementedError("Currently, only implemented for cell submesh")
_label_new_exterior_facets(
dm, subdm, subpoint_indices,
sub_ext_facet_indices, sub_ext_facet_size,
subfStart, subfEnd)
elif subdim == dim - 1:
_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")
Expand Down
Loading
Loading