diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index a4d4d92460..785f91c51b 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -13,6 +13,8 @@ from firedrake.utils import IntType, ScalarType from libc.string cimport memset from libc.stdlib cimport qsort from finat.element_factory import as_fiat_cell +from numbers import Integral +from collections.abc import Sequence cimport numpy as np cimport mpi4py.MPI as MPI @@ -2203,11 +2205,11 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): def _get_periodicity(dm: PETSc.DM) -> tuple[tuple[bool, bool], ...]: """Return mesh periodicity information. - + This function returns a 2-tuple of bools per dimension where the first entry indicates whether the mesh is periodic in that dimension, and the second indicates whether the mesh is single-cell periodic in that dimension. - + """ cdef: const PetscReal *maxCell, *L @@ -3971,7 +3973,7 @@ def create_halo_exchange_sf(PETSc.DM dm): def submesh_create(PETSc.DM dm, PetscInt subdim, label_name, - PetscInt label_value, + subdomain_id, PetscBool ignore_label_halo, comm=None): """Create submesh. @@ -3984,8 +3986,8 @@ def submesh_create(PETSc.DM dm, Topological dimension of the submesh label_name : str Name of the label - label_value : int - Value in the label + subdomain_id : int | Sequence + Values in the label ignore_label_halo : bool If labeled points in the halo are ignored. comm : PETSc.Comm | None @@ -3993,20 +3995,36 @@ def submesh_create(PETSc.DM dm, """ cdef: + PETSc.IS points, subpoints PETSc.DMLabel label, temp_label char *temp_label_name = "firedrake_submesh_temp_label" - PetscInt pStart, pEnd, p, i, stratum_size - PETSc.PetscIS stratum_is = NULL + PetscInt pStart, pEnd, p, i, stratum_size = 0, label_value = 1 const PetscInt *stratum_indices = NULL + # Cast subdomain_id into an iterable + if isinstance(subdomain_id, str) or not isinstance(subdomain_id, Sequence): + subdomain_id = (subdomain_id,) + # Take the union of the all the label values label = dm.getLabel(label_name) + points = PETSc.IS() + for sub in subdomain_id: + if isinstance(sub, Integral): + subpoints = label.getStratumIS(sub) + elif sub == "on_boundary": + subpoints = dm.getStratumIS("exterior_facets", 1) + else: + raise ValueError(f"Submesh construction got invalid subdomain_id {sub}.") + if points: + points = points.union(subpoints) + else: + points = subpoints # Create temp_label that contains no lower-dimensional points. dm.createLabel(temp_label_name) temp_label = dm.getLabel(temp_label_name) - CHKERR(DMLabelGetStratumSize(label.dmlabel, label_value, &stratum_size)) + if points: + CHKERR(ISGetSize(points.iset, &stratum_size)) if stratum_size > 0: - CHKERR(DMLabelGetStratumIS(label.dmlabel, label_value, &stratum_is)) - CHKERR(ISGetIndices(stratum_is, &stratum_indices)) + CHKERR(ISGetIndices(points.iset, &stratum_indices)) CHKERR(DMPlexGetDepthStratum(dm.dm, subdim, &pStart, &pEnd)) for i in range(stratum_size): p = stratum_indices[i] @@ -4014,8 +4032,7 @@ def submesh_create(PETSc.DM dm, # culling all lower-dimensional points. if pStart <= p < pEnd: CHKERR(DMLabelSetValue(temp_label.dmlabel, p, label_value)) - CHKERR(ISRestoreIndices(stratum_is, &stratum_indices)) - CHKERR(ISDestroy(&stratum_is)) + CHKERR(ISRestoreIndices(points.iset, &stratum_indices)) # Make submesh using temp_label. subdm, ownership_transfer_sf = dm.filter(label=temp_label, value=label_value, diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 94952565f8..942eb1c527 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -4884,10 +4884,11 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig subdim : int | None Topological dimension of the submesh. Defaults to ``mesh.topological_dimension``. - subdomain_id : int | None + subdomain_id : int | Sequence | None Subdomain ID representing the submesh. - If `None` the submesh will cover the entire domain. - This is useful to obtain a codim-1 submesh over all facets or + If multiple subdomain IDs are provided, their union is taken. + If `None` the submesh will cover the entire domain, + this is useful to obtain a codim-1 submesh over all facets or a submesh over a different communicator. label_name : str | None Name of the label to search ``subdomain_id`` in. @@ -4930,6 +4931,47 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig ridges to be contained in the quad mesh are shared by at most two facets to make the quad mesh orientation algorithm work. + Examples + -------- + >>> mesh = UnitSquareMesh(4, 4) + >>> x, y = SpatialCoordinate(mesh) + >>> DG = FunctionSpace(mesh, "DG", 0) + >>> DGT = FunctionSpace(mesh, "DGT", 0) + + Mark a cell subdomain and construct a codim-0 submesh from all cells in the subdomain + + >>> cell_marker = assemble(interpolate(conditional(lt(x, 0.5), 1, 0), DG)) + >>> mesh.mark_entities(cell_marker, 111) + >>> submesh = Submesh(mesh, subdomain_id=111) + + Mark a facet subdomain and construct a codim-1 submesh from all facets in the subdomain + + >>> facet_marker = assemble(interpolate(conditional(lt(abs(x-0.5), 1E-12), 1, 0), DGT)) + >>> mesh.mark_entities(facet_marker, 222) + >>> submesh = Submesh(mesh, subdim=mesh.topological_dimension-1, subdomain_id=222) + + Construct a codim-0 submesh of the union of multiple subdomains by passing a list + + >>> mesh.mark_entities(assemble(interpolate(conditional(lt(x, 0.5), 1, 0), DG)), 1) + >>> mesh.mark_entities(assemble(interpolate(conditional(lt(y, 0.5), 1, 0), DG)), 2) + >>> submesh = Submesh(mesh, subdomain_id=[1, 2]) + + Construct a codim-1 submesh of all the facets (the skeleton mesh) + + >>> submesh = Submesh(mesh, subdim=1) + + Construct a codim-1 submesh of the entire boundary + + >>> submesh = Submesh(mesh, subdomain_id="on_boundary") + + Construct a codim-1 submesh of the union of multiple boundaries + + >>> submesh = Submesh(mesh, subdim=mesh.topological_dimension-1, subdomain_id=[1, 2, 3]) + + Construct a codim-0 submesh of the part of the mesh owned by each MPI rank + + >>> submesh = Submesh(mesh, ignore_halo=True, comm=COMM_SELF) + """ if not isinstance(mesh, MeshGeometry): raise TypeError("Parent mesh must be a `MeshGeometry`") @@ -4937,6 +4979,18 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig raise NotImplementedError("Can not create a submesh of an ``ExtrudedMesh``") elif isinstance(mesh.topology, VertexOnlyMeshTopology): raise NotImplementedError("Can not create a submesh of a ``VertexOnlyMesh``") + + if subdomain_id == "on_boundary": + if subdim is None: + subdim = mesh.topological_dimension - 1 + elif subdim != mesh.topological_dimension - 1: + raise ValueError('subdomain_id="on_boundary" requires subdim=dim-1') + if label_name is None: + label_name = "exterior_facets" + elif label_name != "exterior_facets": + raise ValueError('subdomain_id="on_boundary" requires label_name="exterior_facets"') + subdomain_id = 1 + if subdim is None: subdim = mesh.topological_dimension plex = mesh.topology_dm @@ -4962,7 +5016,7 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig if subplex.getDimension() != subdim: raise RuntimeError(f"Found subplex dim ({subplex.getDimension()}) != expected ({subdim})") if reorder is None: - # Ideally we should set perm_is = mesh.dm_reordering[label_indices] + # Ideally we should set perm_is = mesh._dm_renumbering[label_indices] reorder = mesh._did_reordering submesh = Mesh( diff --git a/tests/firedrake/submesh/test_submesh_interface.py b/tests/firedrake/submesh/test_submesh_interface.py new file mode 100644 index 0000000000..394738fab8 --- /dev/null +++ b/tests/firedrake/submesh/test_submesh_interface.py @@ -0,0 +1,45 @@ +import pytest +import numpy as np +from firedrake import * + + +def test_submesh_subdomain_id_union(): + mesh = UnitSquareMesh(4, 4) + x, y = SpatialCoordinate(mesh) + M = FunctionSpace(mesh, "DG", 0) + m1 = Function(M).interpolate(conditional(lt(x, 0.5), 1, 0)) + m2 = Function(M).interpolate(conditional(lt(y, 0.5), 1, 0)) + mesh.mark_entities(m1, 111) + mesh.mark_entities(m2, 222) + + subdomain_id = [111, 222] + submesh1 = Submesh(mesh, mesh.topological_dimension, subdomain_id=subdomain_id) + + m3 = Function(M).interpolate(m1 + m2 - m1 * m2) + expected = assemble(m3*dx) + assert abs(assemble(1*dx(domain=submesh1)) - expected) < 1E-12 + + mesh.mark_entities(m3, 333) + submesh2 = Submesh(mesh, mesh.topological_dimension, 333) + assert submesh2.cell_set.size == submesh1.cell_set.size + assert np.allclose(submesh2.coordinates.dat.data, submesh1.coordinates.dat.data) + + +@pytest.mark.parametrize("subdomain_id", ["on_boundary", (1, 3, 6)]) +def test_submesh_facet_subdomain_id_union(subdomain_id): + mesh = UnitCubeMesh(2, 2, 2) + submesh1 = Submesh(mesh, mesh.topological_dimension - 1, subdomain_id=subdomain_id) + if subdomain_id == "on_boundary": + area = assemble(1*ds(domain=mesh)) + else: + area = assemble(1*ds(subdomain_id, domain=mesh)) + assert abs(assemble(1*dx(domain=submesh1)) - area) < 1E-12 + + DGT = FunctionSpace(mesh, "DGT", 0) + facet_function = Function(DGT) + DirichletBC(DGT, 1, subdomain_id).apply(facet_function) + facet_value = 999 + rmesh = RelabeledMesh(mesh, [facet_function], [facet_value]) + submesh2 = Submesh(rmesh, mesh.topological_dimension - 1, facet_value) + assert submesh2.cell_set.size == submesh1.cell_set.size + assert np.allclose(submesh2.coordinates.dat.data, submesh1.coordinates.dat.data)