diff --git a/firedrake/bcs.py b/firedrake/bcs.py index fad36359f1..74ff402b4b 100644 --- a/firedrake/bcs.py +++ b/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 diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index f440169ded..6e3107e8a9 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/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) diff --git a/firedrake/cython/extrusion_numbering.pyx b/firedrake/cython/extrusion_numbering.pyx index 5c08d14e56..e39aaf9367 100644 --- a/firedrake/cython/extrusion_numbering.pyx +++ b/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. diff --git a/firedrake/functionspacedata.py b/firedrake/functionspacedata.py index 21f74f034f..78209d4357 100644 --- a/firedrake/functionspacedata.py +++ b/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 diff --git a/tests/extrusion/test_variable_layers_numbering.py b/tests/extrusion/test_variable_layers_numbering.py index 45146644cd..47fc93c911 100644 --- a/tests/extrusion/test_variable_layers_numbering.py +++ b/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.