remove geometric boundary mask

test/operations/test_operations_2d-3d.py

@@ -17,16 +17,21 @@ def mesh(mesh2d):
     return utility.extrude_mesh_sigma(mesh2d, n_layers, bathymetry_2d)
 
 
-@pytest.fixture(params=["p1xp1", "P2xP2", "p1DGxp1DG", "P2DGxP2DG"])
+@pytest.fixture(params=[
+    "P1xP1", "P2xP2",
+    "P1DGxP1DG", "P2DGxP2DG",
+    "P1xP3", "P1DGxP3DG",
+    "P3xP1DG", "P3DGxP1",
+])
 def spaces(request):
-    if request.param == "p1xp1":
-        return (("CG", 1), ("CG", 1))
-    elif request.param == "P2xP2":
-        return (("CG", 2), ("CG", 2))
-    elif request.param == "p1DGxp1DG":
-        return (("DG", 1), ("DG", 1))
-    elif request.param == "P2DGxP2DG":
-        return (("DG", 2), ("DG", 2))
+    h_name, v_name = request.param.split('x')
+
+    def get_tuple(name):
+        degree = int(name[1:2])
+        family = 'DG' if 'DG' in name else 'CG'
+        return (family, degree)
+
+    return (get_tuple(h_name), get_tuple(v_name))
 
 
 @pytest.fixture

thetis/forcing.py

@@ -713,7 +713,7 @@ def __init__(self, elev_field, init_date, to_latlon, target_coordsys,
             # interpolate in the whole domain
             self.nodes = np.arange(self.elev_field.dat.data_with_halos.shape[0])
         else:
-            bc = DirichletBC(fs, 0., boundary_ids, method='geometric')
+            bc = DirichletBC(fs, 0., boundary_ids)
             self.nodes = bc.nodes
         self._empty_set = self.nodes.size == 0
 

thetis/limiter.py

@@ -146,8 +146,8 @@ def compute_bounds(self, field):
         if not self.is_2d:
             # Add nodal values from surface/bottom boundaries
             # NOTE calling firedrake par_loop with measure=ds_t raises an error
-            bottom_nodes = get_facet_mask(self.P1CG, 'geometric', 'bottom')
-            top_nodes = get_facet_mask(self.P1CG, 'geometric', 'top')
+            bottom_nodes = get_facet_mask(self.P1CG, 'bottom')
+            top_nodes = get_facet_mask(self.P1CG, 'top')
             bottom_idx = op2.Global(len(bottom_nodes), bottom_nodes, dtype=np.int32, name='node_idx')
             top_idx = op2.Global(len(top_nodes), top_nodes, dtype=np.int32, name='node_idx')
             code = """

thetis/utility.py

@@ -147,6 +147,23 @@ def get_functionspace(mesh, h_family, h_degree, v_family=None, v_degree=None,
     return constructor(mesh, elt, **kwargs)
 
 
+def get_extruded_base_element(ufl_element):
+    """
+    Return UFL TensorProductElement of an extruded UFL element.
+
+    In case of a non-extruded mesh, returns the element itself.
+    """
+    if isinstance(ufl_element, ufl.HDivElement):
+        ufl_element = ufl_element._element
+    if isinstance(ufl_element, ufl.MixedElement):
+        ufl_element = ufl_element.sub_elements()[0]
+    if isinstance(ufl_element, ufl.VectorElement):
+        ufl_element = ufl_element.sub_elements()[0]  # take the first component
+    if isinstance(ufl_element, ufl.EnrichedElement):
+        ufl_element = ufl_element._elements[0]
+    return ufl_element
+
+
 ElementContinuity = namedtuple("ElementContinuity", ["horizontal", "vertical"])
 """
 A named tuple describing the continuity of an element in the horizontal/vertical direction.
@@ -175,21 +192,16 @@ def element_continuity(ufl_element):
         'DQ': 'dg',
     }
 
-    if isinstance(elem, ufl.finiteelement.mixedelement.MixedElement):
-        elem = elem.sub_elements()[0]
-    if isinstance(elem, ufl.finiteelement.mixedelement.VectorElement):
-        elem = elem.sub_elements()[0]  # take the elem of first component
-    if isinstance(elem, ufl.finiteelement.enrichedelement.EnrichedElement):
-        elem = elem._elements[0]
-    if isinstance(elem, ufl.finiteelement.tensorproductelement.TensorProductElement):
-        a, b = elem.sub_elements()
-        horiz_type = elem_types[a.family()]
-        vert_type = elem_types[b.family()]
-    elif isinstance(elem, ufl.finiteelement.hdivcurl.HDivElement):
+    base_element = get_extruded_base_element(ufl_element)
+    if isinstance(elem, ufl.HDivElement):
         horiz_type = 'hdiv'
         vert_type = 'hdiv'
+    elif isinstance(base_element, ufl.TensorProductElement):
+        a, b = base_element.sub_elements()
+        horiz_type = elem_types[a.family()]
+        vert_type = elem_types[b.family()]
     else:
-        horiz_type = elem_types[elem.family()]
+        horiz_type = elem_types[base_element.family()]
         vert_type = horiz_type
     return ElementContinuity(horiz_type, vert_type)
 
@@ -210,25 +222,32 @@ def create_directory(path, comm=COMM_WORLD):
     return path
 
 
-def get_facet_mask(function_space, mode='geometric', facet='bottom'):
+def get_facet_mask(function_space, facet='bottom'):
     """
     Returns the top/bottom nodes of extruded 3D elements.
 
     :arg function_space: Firedrake :class:`FunctionSpace` object
-    :kwarg str mode: 'topological', to retrieve nodes that lie on the facet, or
-        'geometric' for nodes whose basis functions do not vanish on the facet.
     :kwarg str facet: 'top' or 'bottom'
 
     .. note::
         The definition of top/bottom depends on the direction of the extrusion.
         Here we assume that the mesh has been extruded upwards (along positive
         z axis).
     """
-    section, iset, facets = function_space.cell_boundary_masks[mode]
-    ifacet = -2 if facet == 'bottom' else -1
-    off = section.getOffset(facets[ifacet])
-    dof = section.getDof(facets[ifacet])
-    indices = iset[off:off+dof]
+    from tsfc.finatinterface import create_element as create_finat_element
+    # get base element
+    elem = get_extruded_base_element(function_space.ufl_element())
+    assert isinstance(elem, TensorProductElement), \
+        f'function space must be defined on an extruded 3D mesh: {elem}'
+    # figure out number of nodes in sub elements
+    h_elt, v_elt = elem.sub_elements()
+    nb_nodes_h = create_finat_element(h_elt).space_dimension()
+    nb_nodes_v = create_finat_element(v_elt).space_dimension()
+    # compute top/bottom facet indices
+    # extruded dimension is the inner loop in index
+    # on interval elements, the end points are the first two dofs
+    offset = 0 if facet == 'bottom' else 1
+    indices = np.arange(nb_nodes_h)*nb_nodes_v + offset
     return indices
 
 
@@ -266,7 +285,7 @@ def extrude_mesh_sigma(mesh2d, n_layers, bathymetry_2d, z_stretch_fact=1.0,
         for i, v in enumerate(min_depth):
             min_depth_arr[i] = v
 
-    nodes = get_facet_mask(fs_3d, 'geometric', 'bottom')
+    nodes = get_facet_mask(fs_3d, 'bottom')
     idx = op2.Global(len(nodes), nodes, dtype=np.int32, name='node_idx')
     min_depth_op2 = op2.Global(len(min_depth_arr), min_depth_arr, name='min_depth')
     kernel = op2.Kernel("""
@@ -758,7 +777,7 @@ def __init__(self, solver_obj):
         self.fs_3d = self.solver_obj.function_spaces.P1DG
         assert self.output.function_space() == self.fs_3d
 
-        nodes = get_facet_mask(self.fs_3d, 'geometric', 'bottom')
+        nodes = get_facet_mask(self.fs_3d, 'bottom')
         self.idx = op2.Global(len(nodes), nodes, dtype=np.int32, name='node_idx')
         self.kernel = op2.Kernel("""
             void my_kernel(double *output, double *z_field, int *idx) {
@@ -884,18 +903,11 @@ def __init__(self, input_2d, output_3d, elem_height=None):
         self.fs_3d = self.output_3d.function_space()
 
         family_2d = self.fs_2d.ufl_element().family()
-        ufl_elem = self.fs_3d.ufl_element()
-        if isinstance(ufl_elem, ufl.VectorElement):
-            # Unwind vector
-            ufl_elem = ufl_elem.sub_elements()[0]
-        if isinstance(ufl_elem, ufl.HDivElement):
-            # RT case
-            ufl_elem = ufl_elem._element
-        if ufl_elem.family() == 'TensorProductElement':
-            # a normal tensorproduct element
-            family_3dh = ufl_elem.sub_elements()[0].family()
-            if family_2d != family_3dh:
-                raise Exception('2D and 3D spaces do not match: {0:s} {1:s}'.format(family_2d, family_3dh))
+        base_element_3d = get_extruded_base_element(self.fs_3d.ufl_element())
+        assert isinstance(base_element_3d, ufl.TensorProductElement)
+        family_3dh = base_element_3d.sub_elements()[0].family()
+        if family_2d != family_3dh:
+            raise Exception('2D and 3D spaces do not match: {0:s} {1:s}'.format(family_2d, family_3dh))
         self.do_hdiv_scaling = family_2d in ['Raviart-Thomas', 'RTCF', 'Brezzi-Douglas-Marini', 'BDMCF']
         if self.do_hdiv_scaling and elem_height is None:
             raise Exception('elem_height must be provided for HDiv spaces')
@@ -905,7 +917,7 @@ def __init__(self, input_2d, output_3d, elem_height=None):
         # number of nodes in vertical direction
         n_vert_nodes = self.fs_3d.finat_element.space_dimension() / self.fs_2d.finat_element.space_dimension()
 
-        nodes = get_facet_mask(self.fs_3d, 'geometric', 'bottom')
+        nodes = get_facet_mask(self.fs_3d, 'bottom')
         self.idx = op2.Global(len(nodes), nodes, dtype=np.int32, name='node_idx')
         self.kernel = op2.Kernel("""
             void my_kernel(double *func, double *func2d, int *idx) {
@@ -1012,26 +1024,21 @@ def __init__(self, input_3d, output_2d,
             elem_facet = boundary
 
         family_2d = self.fs_2d.ufl_element().family()
-        elem = self.fs_3d.ufl_element()
-        if isinstance(elem, ufl.VectorElement):
-            elem = elem.sub_elements()[0]
-        if isinstance(elem, ufl.HDivElement):
-            elem = elem._element
-        if isinstance(elem, ufl.TensorProductElement):
-            # a normal tensorproduct element
-            family_3dh = elem.sub_elements()[0].family()
-            if family_2d != family_3dh:
-                raise Exception('2D and 3D spaces do not match: {0:s} {1:s}'.format(family_2d, family_3dh))
+        base_element_3d = get_extruded_base_element(self.fs_3d.ufl_element())
+        assert isinstance(base_element_3d, ufl.TensorProductElement)
+        family_3dh = base_element_3d.sub_elements()[0].family()
+        if family_2d != family_3dh:
+            raise Exception('2D and 3D spaces do not match: {0:s} {1:s}'.format(family_2d, family_3dh))
         self.do_hdiv_scaling = family_2d in ['Raviart-Thomas', 'RTCF', 'Brezzi-Douglas-Marini', 'BDMCF']
         if self.do_hdiv_scaling and elem_height is None:
             raise Exception('elem_height must be provided for HDiv spaces')
 
         assert elem_facet in ['top', 'bottom', 'average'], 'Unsupported elem_facet: {:}'.format(elem_facet)
         if elem_facet == 'average':
-            nodes = np.hstack((get_facet_mask(self.fs_3d, 'geometric', 'bottom'),
-                               get_facet_mask(self.fs_3d, 'geometric', 'top')))
+            nodes = np.hstack((get_facet_mask(self.fs_3d, 'bottom'),
+                               get_facet_mask(self.fs_3d, 'top')))
         else:
-            nodes = get_facet_mask(self.fs_3d, 'geometric', elem_facet)
+            nodes = get_facet_mask(self.fs_3d, elem_facet)
         if boundary == 'top':
             self.iter_domain = op2.ON_TOP
         elif boundary == 'bottom':
@@ -1398,23 +1405,16 @@ def __init__(self, solver):
         self.fs_2d = self.fields.elev_cg_2d.function_space()
 
         family_2d = self.fs_2d.ufl_element().family()
-        ufl_elem = self.fs_3d.ufl_element()
-        if isinstance(ufl_elem, ufl.VectorElement):
-            # Unwind vector
-            ufl_elem = ufl_elem.sub_elements()[0]
-        if isinstance(ufl_elem, ufl.HDivElement):
-            # RT case
-            ufl_elem = ufl_elem._element
-        if ufl_elem.family() == 'TensorProductElement':
-            # a normal tensorproduct element
-            family_3dh = ufl_elem.sub_elements()[0].family()
-            if family_2d != family_3dh:
-                raise Exception('2D and 3D spaces do not match: "{0:s}" != "{1:s}"'.format(family_2d, family_3dh))
+        base_element_3d = get_extruded_base_element(self.fs_3d.ufl_element())
+        assert isinstance(base_element_3d, ufl.TensorProductElement)
+        family_3dh = base_element_3d.sub_elements()[0].family()
+        if family_2d != family_3dh:
+            raise Exception('2D and 3D spaces do not match: "{0:s}" != "{1:s}"'.format(family_2d, family_3dh))
 
         # number of nodes in vertical direction
         n_vert_nodes = self.fs_3d.finat_element.space_dimension() / self.fs_2d.finat_element.space_dimension()
 
-        nodes = get_facet_mask(self.fs_3d, 'geometric', 'bottom')
+        nodes = get_facet_mask(self.fs_3d, 'bottom')
         self.idx = op2.Global(len(nodes), nodes, dtype=np.int32, name='node_idx')
         self.kernel_z_coord = op2.Kernel("""
             void my_kernel(double *z_coord_3d, double *z_ref_3d, double *elev_2d, double *bath_2d, int *idx) {