Merge pull request #3049 from firedrakeproject/ksagiyam/relabeledmesh…

…_hex_workaround

mesh: allow for marking facets with Q2 for hex for now

firedrake/mesh.py

@@ -925,6 +925,8 @@ def mark_entities(self, tf, label_value, label_name=None):
             must be "DP" or "DQ" (degree 0) to mark cell entities and
             "P" (degree 1) in 1D or "HDiv Trace" (degree 0) in 2D or 3D
             to mark facet entities.
+            Can use "Q" (degree 2) functions for 3D hex meshes until
+            we support "HDiv Trace" elements on hex.
         :arg lable_value: The value used in the label.
         :arg label_name: The name of the label to store entity selections.
 
@@ -1315,7 +1317,8 @@ def mark_entities(self, tf, label_value, label_name=None):
             height = 0
             label_name = label_name or dmcommon.CELL_SETS_LABEL
         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):
+                (elem.family() == "Lagrange" and elem.degree() == 1 and self.cell_dimension() == 1) or \
+                (elem.family() == "Q" and elem.degree() == 2 and self.ufl_cell().cellname() == "hexahedron"):
             # facets
             height = 1
             label_name = label_name or dmcommon.FACE_SETS_LABEL
@@ -2427,6 +2430,8 @@ def mark_entities(self, f, label_value, label_name=None):
             must be "DP" or "DQ" (degree 0) to mark cell entities and
             "P" (degree 1) in 1D or "HDiv Trace" (degree 0) in 2D or 3D
             to mark facet entities.
+            Can use "Q" (degree 2) functions for 3D hex meshes until
+            we support "HDiv Trace" elements on hex.
         :arg lable_value: The value used in the label.
         :arg label_name: The name of the label to store entity selections.
 
@@ -4046,6 +4051,8 @@ def RelabeledMesh(mesh, indicator_functions, subdomain_ids, **kwargs):
         "DP"/"DQ" (degree 0) functions to mark cell entities and
         "P" (degree 1) functions in 1D or "HDiv Trace" (degree 0) functions
         in 2D or 3D to mark facet entities.
+        Can use "Q" (degree 2) functions for 3D hex meshes until
+        we support "HDiv Trace" elements on hex.
     :arg subdomain_ids: list of subdomain ids associated with
         the indicator functions in indicator_functions; thus,
         must have the same length as indicator_functions.
@@ -4098,7 +4105,8 @@ def RelabeledMesh(mesh, indicator_functions, subdomain_ids, **kwargs):
             height = 0
             dmlabel_name = dmcommon.CELL_SETS_LABEL
         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):
+                (elem.family() == "Lagrange" and elem.degree() == 1 and mesh.topological_dimension() == 1) or \
+                (elem.family() == "Q" and elem.degree() == 2 and mesh.topology.ufl_cell().cellname() == "hexahedron"):
             # facets
             height = 1
             dmlabel_name = dmcommon.FACE_SETS_LABEL

tests/regression/test_mark_entities.py

@@ -70,6 +70,20 @@ def test_mark_entities_overlapping_facet_subdomains():
     assert abs(v - 1.0) < 1.e-10
 
 
+def test_mark_entities_mesh_mark_entities_3d_hex():
+    label_name = "test_label"
+    label_value = 999
+    mesh = UnitCubeMesh(2, 2, 2, hexahedral=True)
+    x, y, z = SpatialCoordinate(mesh)
+    V = FunctionSpace(mesh, "Q", 2)
+    f = Function(V).interpolate(conditional(And(x > .6, And(y > .4, y < .6)), 1, 0))
+    mesh.mark_entities(f, label_value, label_name=label_name)
+    plex = mesh.topology.topology_dm
+    label = plex.getLabel(label_name)
+    assert label.getStratumIS(label_value).getSize() == 2
+    assert all(label.getStratumIS(label_value).getIndices() == [51, 57])
+
+
 def test_mark_entities_mesh_mark_entities_2d():
     # +---+---+
     # | \ | \ |  mark        \