Extend facet marking to 1D meshes (#2802)

* Extend facet marking to 1D meshes
firedrakeproject · Mar 9, 2023 · f13bbb5 · f13bbb5
1 parent 4779660
commit f13bbb5
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 8 deletions.
diff --git 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
@@ -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])