From f13bbb50233a03fb757716dd087ae8b6ff4f546d Mon Sep 17 00:00:00 2001 From: Francesco Ballarin Date: Thu, 9 Mar 2023 14:24:47 +0100 Subject: [PATCH] Extend facet marking to 1D meshes (#2802) * Extend facet marking to 1D meshes --- firedrake/mesh.py | 19 ++++++++++++------- tests/regression/test_mark_entities.py | 16 +++++++++++++++- 2 files changed, 27 insertions(+), 8 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 99b170369b..b44d4ba298 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -822,7 +822,8 @@ def mark_entities(self, tf, label_name, label_value): :arg tf: The :class:`.CoordinatelessFunction` object that marks selected entities as 1. f.function_space().ufl_element() must be "DP" or "DQ" (degree 0) to mark cell entities and - "HDiv Trace " (degree 0) to mark facet entities. + "P" (degree 1) in 1D or "HDiv Trace" (degree 0) in 2D or 3D + to mark facet entities. :arg label_name: The name of the label to store entity selections. :arg lable_value: The value used in the label. @@ -1263,11 +1264,12 @@ def mark_entities(self, tf, label_name, label_value): if elem.family() in {"Discontinuous Lagrange", "DQ"} and elem.degree() == 0: # cells height = 0 - elif elem.family() == "HDiv Trace" and elem.degree() == 0: + elif (elem.family() == "HDiv Trace" and elem.degree() == 0 and self.cell_dimension() > 1) or \ + (elem.family() == "Lagrange" and elem.degree() == 1 and self.cell_dimension() == 1): # facets height = 1 else: - raise ValueError(f"indicator functions must be 'DP' or 'DQ' (degree 0) to mark cells and 'HDiv Trace' (degree 0) to mark facets: got (family, degree) = ({elem.family()}, {elem.degree()})") + raise ValueError(f"indicator functions must be 'DP' or 'DQ' (degree 0) to mark cells and 'P' (degree 1) in 1D or 'HDiv Trace' (degree 0) in 2D or 3D to mark facets: got (family, degree) = ({elem.family()}, {elem.degree()})") plex = self.topology_dm if not plex.hasLabel(label_name): plex.createLabel(label_name) @@ -2162,7 +2164,8 @@ def mark_entities(self, f, label_name, label_value): :arg f: The :class:`.Function` object that marks selected entities as 1. f.function_space().ufl_element() must be "DP" or "DQ" (degree 0) to mark cell entities and - "HDiv Trace " (degree 0) to mark facet entities. + "P" (degree 1) in 1D or "HDiv Trace" (degree 0) in 2D or 3D + to mark facet entities. :arg label_name: The name of the label to store entity selections. :arg lable_value: The value used in the label. @@ -2968,7 +2971,8 @@ def RelabeledMesh(mesh, indicator_functions, subdomain_ids, **kwargs): :arg indicator_functions: list of indicator functions that mark selected entities (cells or facets) as 1; must use "DP"/"DQ" (degree 0) functions to mark cell entities and - "HDiv Trace" (degree 0) functions to mark facet entities. + "P" (degree 1) functions in 1D or "HDiv Trace" (degree 0) functions + in 2D or 3D to mark facet entities. :arg subdomain_ids: list of subdomain ids associated with the indicator functions in indicator_functions; thus, must have the same length as indicator_functions. @@ -3020,12 +3024,13 @@ def RelabeledMesh(mesh, indicator_functions, subdomain_ids, **kwargs): # cells height = 0 dmlabel_name = dmcommon.CELL_SETS_LABEL - elif elem.family() == "HDiv Trace" and elem.degree() == 0: + elif (elem.family() == "HDiv Trace" and elem.degree() == 0 and mesh.topological_dimension() > 1) or \ + (elem.family() == "Lagrange" and elem.degree() == 1 and mesh.topological_dimension() == 1): # facets height = 1 dmlabel_name = dmcommon.FACE_SETS_LABEL else: - raise ValueError(f"indicator functions must be 'DP' or 'DQ' (degree 0) to mark cells and 'HDiv Trace' (degree 0) to mark facets: got (family, degree) = ({elem.family()}, {elem.degree()})") + raise ValueError(f"indicator functions must be 'DP' or 'DQ' (degree 0) to mark cells and 'P' (degree 1) in 1D or 'HDiv Trace' (degree 0) in 2D or 3D to mark facets: got (family, degree) = ({elem.family()}, {elem.degree()})") # Clear label stratum; this is a copy, so safe to change. plex1.clearLabelStratum(dmlabel_name, subid) dmlabel = plex1.getLabel(dmlabel_name) diff --git a/tests/regression/test_mark_entities.py b/tests/regression/test_mark_entities.py index ec027a3368..0c0cf5b40e 100644 --- a/tests/regression/test_mark_entities.py +++ b/tests/regression/test_mark_entities.py @@ -70,7 +70,7 @@ def test_mark_entities_overlapping_facet_subdomains(): assert abs(v - 1.0) < 1.e-10 -def test_mark_entities_mesh_mark_entities(): +def test_mark_entities_mesh_mark_entities_2d(): # +---+---+ # | \ | \ | mark \ # +---+---+ -----> +---+ @@ -90,3 +90,17 @@ def test_mark_entities_mesh_mark_entities(): label = plex.getLabel(label_name) assert label.getStratumIS(label_value).getSize() == 2 assert all(label.getStratumIS(label_value).getIndices() == [20, 30]) + + +def test_mark_entities_mesh_mark_entities_1d(): + label_name = "test_label" + label_value = 999 + mesh = UnitIntervalMesh(2) + x, = SpatialCoordinate(mesh) + V = FunctionSpace(mesh, "P", 1) + f = Function(V).interpolate(conditional(x < 0.25, 1., 0.)) + mesh.mark_entities(f, label_name, label_value) + plex = mesh.topology.topology_dm + label = plex.getLabel(label_name) + assert label.getStratumIS(label_value).getSize() == 1 + assert all(label.getStratumIS(label_value).getIndices() == [2])