From 4b5593fe5c1d6effec30461f88c109e7b7d8b283 Mon Sep 17 00:00:00 2001 From: ksagiyam Date: Fri, 28 Jul 2023 10:31:56 +0100 Subject: [PATCH] mesh: allow for marking facets with Q2 for hex for now --- firedrake/mesh.py | 12 ++++++++++-- tests/regression/test_mark_entities.py | 14 ++++++++++++++ 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 1bbbf32d91..366f2f1862 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -835,6 +835,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. @@ -1278,7 +1280,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 @@ -2363,6 +2366,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. @@ -3980,6 +3985,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. @@ -4032,7 +4039,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 diff --git a/tests/regression/test_mark_entities.py b/tests/regression/test_mark_entities.py index 3045ebb9fe..5741ac5cfa 100644 --- a/tests/regression/test_mark_entities.py +++ b/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 \