diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index ec7ab631af..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 @@ -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): @@ -3827,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, @@ -3937,77 +3941,288 @@ 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): + 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 *sub_ext_facet_indices, + PetscInt sub_ext_facet_size, + PetscInt subfStart, PetscInt subfEnd, +): + """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 + 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" + elif subdim == 1: + source_label_name = b"Vertex Sets" + else: + source_label_name = None + + # 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(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 = 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: + # --- 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) @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``): 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: - 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 = "interior_facets" - char *ext_facet_label_name = "exterior_facets" - char *face_sets_label_name = "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") - CHKERR(ISRestoreIndices(subpoint_is.iset, &subpoint_indices)) + 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, + 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)) subdm.removeLabel("interior_facets") subdm.removeLabel("exterior_facets") 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 new file mode 100644 index 0000000000..f2640d9808 --- /dev/null +++ b/tests/firedrake/submesh/test_submesh_codim1_labels.py @@ -0,0 +1,279 @@ +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, 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(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) + + 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(ncells=2): + """UnitSquareMesh with Face Sets on y=0 boundary edges and Vertex + Sets on the corner vertices of that edge. + """ + 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) + + 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) +# --------------------------------------------------------------------------- + +@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 + 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 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) + # 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)}" + ) + + +@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 + # 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 + if dm.hasLabel("Face Sets"): + dm.removeLabel("Face Sets") + dm.createLabel("Face Sets") + 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 + assert len(markers) > 0 + + +# --------------------------------------------------------------------------- +# Parallel – coordinate expression integration +# --------------------------------------------------------------------------- + +@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=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 + 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 + + +# --------------------------------------------------------------------------- +# Parallel – FunctionSpace interpolation + boundary integration +# --------------------------------------------------------------------------- + +@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=4) + 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 + + +# --------------------------------------------------------------------------- +# Parallel – 2D -> 1D (Vertex Sets) +# --------------------------------------------------------------------------- + +@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=4) + 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