bcs: Handle interior markers and simplify construction of nodes

Now that we're only handling topological boundary conditions, we can
just inspect the section of the function space to determine which
nodes to mask out. This also enables us to easily support constraining
dofs that live on interior facets, so do that too.

Fixes #1661, #1135.

firedrake/bcs.py

@@ -234,6 +234,13 @@ def extract_forms(self, form_type):
 class DirichletBC(BCBase, DirichletBCMixin):
     r'''Implementation of a strong Dirichlet boundary condition.
 
+    .. note:
+
+       This uses facet markers in the domain, so may be used to
+       applied strong boundary conditions to interior facets (if they
+       have an appropriate mesh marker). The "on_boundary" string only
+       applies to the exterior boundaries of the domain.
+
     :arg V: the :class:`.FunctionSpace` on which the boundary condition
         should be applied.
     :arg g: the boundary condition values. This can be a :class:`.Function` on
@@ -251,6 +258,7 @@ class DirichletBC(BCBase, DirichletBCMixin):
     :arg method: the method for determining boundary nodes.
         DEPRECATED. The only way boundary nodes are identified is by
         topological association.
+
     '''
 
     @DirichletBCMixin._ad_annotate_init

firedrake/cython/dmcommon.pyx

@@ -900,89 +900,81 @@ def get_facet_nodes(mesh, np.ndarray[PetscInt, ndim=2, mode="c"] cell_nodes, lab
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def boundary_nodes(V, sub_domain):
-    """Extract boundary nodes from a function space..
+def facet_closure_nodes(V, sub_domain):
+    """Extract nodes in the closure of facets with a given marker.
+
+    This works fine for interior as well as exterior facets.
 
     :arg V: the function space
-    :arg sub_domain: a mesh marker selecting the part of the boundary
-        (may be "on_boundary" indicating the entire boundary).
-    :returns: a numpy array of unique nodes on the boundary of the
-        requested subdomain.
+    :arg sub_domain: a tuple of mesh markers selecting the facets, or
+        the magic string "on_boundary" indicating the entire boundary.
+    :returns: a numpy array of unique nodes in the closure of facets
+        with the given marker.
     """
     cdef:
-        np.ndarray[np.int32_t, ndim=2, mode="c"] local_nodes
-        np.ndarray[PetscInt, ndim=1, mode="c"] offsets
-        np.ndarray[np.uint32_t, ndim=1, mode="c"] local_facets
+        PETSc.Section sec = V.dm.getSection()
+        PETSc.DM dm = V.mesh().topology_dm
+        PetscInt nnodes, p, i, dof, offset, n, j, d
+        np.ndarray[PetscInt, ndim=1, mode="c"] points
         np.ndarray[PetscInt, ndim=1, mode="c"] nodes
-        np.ndarray[PetscInt, ndim=2, mode="c"] facet_node_list
-        np.ndarray[PetscInt, ndim=2, mode="c"] layer_extents
-        np.ndarray[PetscInt, ndim=1, mode="c"] facet_indices
-        int f, i, j, dof, facet, idx
-        int nfacet, nlocal, layers
-        PetscInt local_facet
-        PetscInt offset
-        bint all_facets
-        bint variable, extruded
-
-    mesh = V.mesh()
-    variable = mesh.variable_layers
-    extruded = mesh.cell_set._extruded
-
-    facet_dim = mesh.facet_dimension()
-    boundary_dofs = V.finat_element.entity_closure_dofs()[facet_dim]
-
-    local_nodes = np.empty((len(boundary_dofs),
-                            len(boundary_dofs[0])),
-                           dtype=np.int32)
-    for k, v in boundary_dofs.items():
-        local_nodes[k, :] = v
-
-    facets = V.mesh().exterior_facets
-    local_facets = facets.local_facet_dat.data_ro_with_halos
-    nlocal = local_nodes.shape[1]
-
     if sub_domain == "on_boundary":
-        subset = facets.set
-        all_facets = True
+        label = "exterior_facets"
+        sub_domain = (1, )
     else:
-        all_facets = False
-        subset = facets.subset(sub_domain)
-        facet_indices = subset.indices
-
-    nfacet = subset.total_size
-    offsets = V.offset
-    facet_node_list = V.exterior_facet_node_map().values_with_halo
-
-    if variable:
-        layer_extents = subset.layers_array
-        maxsize = local_nodes.shape[1] * np.sum((layer_extents[:, 1] - layer_extents[:, 0]) - 1)
-    elif extruded:
-        layers = subset.layers
-        maxsize = local_nodes.shape[1] * nfacet * (layers - 1)
-    else:
-        layers = 2
-        offset = 0
-        maxsize = local_nodes.shape[1] * nfacet
+        label = FACE_SETS_LABEL
+        if V.mesh().variable_layers:
+            # We can't use the generic code in this case because we
+            # need to manually take closure of facets on each external
+            # face (rather than using the label completion and section
+            # information).
+            # The reason for this is that (for example) a stack of
+            # vertices may extend higher than the faces
+            #
+            #            Y---.
+            #            |   |
+            #    .---x---x---.
+            #    |   | O |
+            #    .---x---x
+            #    |   |
+            #    .---Y
+            #
+            # BCs on the facet marked with 'O' should produce 4 values
+            # for a P1 field (marked X). But if we just include
+            # everything from the sections we'd get 6: additionally we
+            # see the points marked Y.
+            from firedrake.cython import extrusion_numbering as extnum
+            return extnum.facet_closure_nodes(V, sub_domain)
+
+    if not dm.hasLabel(label) or all(dm.getStratumSize(label, marker)
+                                     for marker in sub_domain) == 0:
+        return np.empty(0, dtype=IntType)
 
-    nodes = np.empty(maxsize, dtype=IntType)
-    idx = 0
-    for f in range(nfacet):
-        if all_facets:
-            facet = f
-        else:
-            facet = facet_indices[f]
-        local_facet = local_facets[facet]
-        if variable:
-            layers = layer_extents[f, 1] - layer_extents[f, 0]
-        for i in range(nlocal):
-            dof = local_nodes[local_facet, i]
-            for j in range(layers - 1):
-                if extruded:
-                    offset = j * offsets[dof]
-                nodes[idx] = facet_node_list[facet, dof] + offset
-                idx += 1
-
-    assert idx == nodes.shape[0]
+    nnodes = 0
+    for marker in sub_domain:
+        n = dm.getStratumSize(label, marker)
+        if n == 0:
+            continue
+        points = dm.getStratumIS(label, marker).indices
+        for i in range(n):
+            p = points[i]
+            CHKERR(PetscSectionGetDof(sec.sec, p, &dof))
+            nnodes += dof
+
+    nodes = np.empty(nnodes, dtype=IntType)
+    j = 0
+    for marker in sub_domain:
+        n = dm.getStratumSize(label, marker)
+        if n == 0:
+            continue
+        points = dm.getStratumIS(label, marker).indices
+        for i in range(n):
+            p = points[i]
+            CHKERR(PetscSectionGetDof(sec.sec, p, &dof))
+            CHKERR(PetscSectionGetOffset(sec.sec, p, &offset))
+            for d in range(dof):
+                nodes[j] = offset + d
+                j += 1
+    assert j == nnodes
     return np.unique(nodes)
 
 

firedrake/cython/extrusion_numbering.pyx

@@ -359,6 +359,114 @@ def node_classes(mesh, nodes_per_entity):
     return numpy.cumsum(node_classes)
 
 
+@cython.wraparound(False)
+def facet_closure_nodes(V, sub_domain):
+    """Extract nodes in the closure of facets with a given marker.
+
+    This works fine for interior as well as exterior facets.
+
+    .. note::
+       Don't call this function directly, but rather call
+       :func:`~.dmcommon.facet_closure_nodes`, which will dispatch
+       here if appropriate.
+
+    :arg V: the function space
+    :arg sub_domain: a mesh marker selecting the part of the boundary
+    :returns: a numpy array of unique nodes on the boundary of the
+        requested subdomain.
+    """
+    cdef:
+        numpy.ndarray[numpy.int32_t, ndim=2, mode="c"] local_nodes
+        numpy.ndarray[PetscInt, ndim=1, mode="c"] offsets
+        numpy.ndarray[numpy.uint32_t, ndim=1] local_facets
+        numpy.ndarray[PetscInt, ndim=1, mode="c"] nodes
+        numpy.ndarray[PetscInt, ndim=2] facet_node_list
+        numpy.ndarray[PetscInt, ndim=2, mode="c"] layer_extents
+        numpy.ndarray[PetscInt, ndim=1, mode="c"] facet_indices
+        int f, i, j, dof, facet, idx
+        int nfacet, nlocal, layers
+        PetscInt local_facet
+        PetscInt offset
+
+    # We don't have to handle the "on_boundary" case, because the
+    # caller handles it.
+    mesh = V.mesh()
+    facet_dim = mesh.facet_dimension()
+    boundary_dofs = V.finat_element.entity_closure_dofs()[facet_dim]
+
+    local_nodes = numpy.empty((len(boundary_dofs),
+                               len(boundary_dofs[0])),
+                              dtype=numpy.int32)
+    for k, v in boundary_dofs.items():
+        local_nodes[k, :] = v
+
+    all_nodes = []
+    nlocal = local_nodes.shape[1]
+    offsets = V.offset
+    # Walk over both facet types
+    for typ in ["exterior", "interior"]:
+        if typ == "exterior":
+            facets = V.mesh().exterior_facets
+            local_facets = facets.local_facet_dat.data_ro_with_halos
+            facet_node_list = V.exterior_facet_node_map().values_with_halo
+        elif typ == "interior":
+            facets = V.mesh().interior_facets
+            local_facets = facets.local_facet_dat.data_ro_with_halos[:, 0]
+            facet_node_list = V.interior_facet_node_map().values_with_halo
+            facet_node_list = facet_node_list[:, :V.finat_element.space_dimension()]
+
+        subset = facets.subset(sub_domain)
+        facet_indices = subset.indices
+
+        nfacet = subset.total_size
+
+        layer_extents = subset.layers_array
+        maxsize = local_nodes.shape[1] * numpy.sum((layer_extents[:, 1]
+                                                    - layer_extents[:, 0]) - 1)
+        nodes = numpy.empty(maxsize, dtype=IntType)
+        idx = 0
+        for f in range(nfacet):
+            # For each facet, pick up all dofs in the closure
+            facet = facet_indices[f]
+            local_facet = local_facets[facet]
+            layers = layer_extents[f, 1] - layer_extents[f, 0]
+            for i in range(nlocal):
+                dof = local_nodes[local_facet, i]
+                for j in range(layers - 1):
+                    nodes[idx] = facet_node_list[facet, dof] + j*offsets[dof]
+                    idx += 1
+
+        assert idx == nodes.shape[0]
+        all_nodes.append(nodes)
+    nodes = numpy.unique(numpy.concatenate(all_nodes))
+    # We need a halo exchange to determine all bc nodes.
+    # Consider
+    # +----+----+
+    # |\ 1 | 2 /
+    # | \  |  /
+    # |  \ | /
+    # | 0 \|/
+    # +----+
+    # With rank 0 owning cell 0 and rank 1 owning cells 1 and 2.
+    # Imagine now applying a DirichletBC on the right-most facet. That
+    # means that the bottom right node (in a CG1 function space) is
+    # killed. But rank 0 doesn't know that that dof is on a boundary
+    # (because it only sees cell 1 which does not have an external
+    # facet attached to that node).
+    # For all the other bcs, the topological completion of labels to
+    # all mesh points works. But for variable layer extruded meshes,
+    # we need to do this by hand.
+    # See github.com/firedrakeproject/firedrake/issues/1135 for even
+    # more details.
+    d = op2.Dat(V.dof_dset.set, dtype=numpy.int8)
+    d.data_with_halos[nodes] = 1
+    d.global_to_local_begin(op2.READ)
+    d.global_to_local_end(op2.READ)
+    indices, = numpy.where(d.data_ro_with_halos == 1)
+    # cast, because numpy.where returns an int64
+    return indices.astype(IntType)
+
+
 @cython.wraparound(False)
 def entity_layers(mesh, height, label=None):
     """Compute the layers for a given entity type.

firedrake/functionspacedata.py

@@ -297,18 +297,23 @@ def get_top_bottom_boundary_nodes(mesh, key, V):
 
 
 @cached
-def get_boundary_nodes(mesh, key, V):
+def get_facet_closure_nodes(mesh, key, V):
+    """Function space nodes in the closure of facets with a given
+    marker.
+    :arg mesh: Mesh to cache on
+    :arg key: (edofs, sub_domain) tuple
+    :arg V: function space.
+    :returns: numpy array of unique nodes in the closure of facets
+       with provided markers (both interior and exterior)."""
     _, sub_domain = key
-    indices = dmcommon.boundary_nodes(V, sub_domain)
-    # We need a halo exchange to determine all bc nodes.
-    # Should be improved by doing this on the DM topology once.
-    d = op2.Dat(V.dof_dset.set, dtype=IntType)
-    d.data_with_halos[indices] = 1
-    d.global_to_local_begin(op2.READ)
-    d.global_to_local_end(op2.READ)
-    indices, = numpy.where(d.data_ro_with_halos == 1)
-    # cast, because numpy where returns an int64
-    return indices.astype(IntType)
+    if sub_domain not in {"on_boundary", "top", "bottom"}:
+        valid = set(mesh.interior_facets.unique_markers)
+        valid |= set(mesh.exterior_facets.unique_markers)
+        invalid = set(sub_domain) - valid
+        if invalid:
+            raise LookupError(f"BC construction got invalid markers {invalid}. "
+                              f"Valid markers are '{valid}'")
+    return dmcommon.facet_closure_nodes(V, sub_domain)
 
 
 def get_max_work_functions(V):
@@ -433,7 +438,7 @@ def boundary_nodes(self, V, sub_domain):
             else:
                 sdkey = as_tuple(sub_domain)
             key = (entity_dofs_key(V.finat_element.entity_dofs()), sdkey)
-            return get_boundary_nodes(V.mesh(), key, V)
+            return get_facet_closure_nodes(V.mesh(), key, V)
 
     def get_map(self, V, entity_set, map_arity, name, offset):
         """Return a :class:`pyop2.Map` from some topological entity to

tests/extrusion/test_variable_layers_numbering.py

@@ -294,7 +294,9 @@ def test_numbering_two_d_bigger():
     for faces, val in [((11, 13), 1),
                        ((14, 20), 2),
                        ((16, ), 3),
-                       ((17, 18, 19), 4)]:
+                       ((17, 18, 19), 4),
+                       # This one is an interior face
+                       ((12, ), 5)]:
         for face in faces:
             dm.setLabelValue("Face Sets", face, val)
 
@@ -343,6 +345,10 @@ def test_numbering_two_d_bigger():
     assert numpy.equal(DirichletBC(V, 0, 4).nodes,
                        [12, 13, 14, 15, 16, 17, 18, 19, 20]).all()
 
+    # Interior face between base plex cells 0 and 1.
+    assert numpy.equal(DirichletBC(V, 0, 5).nodes,
+                       [4, 5, 8, 9]).all()
+
 
 def test_numbering_quad():
     # Number of cells in each column.