diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 5b6d1f1fb1..4e957affc9 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst @@ -4,7 +4,8 @@ Oceanic Basin Modes: Quasi-Geostrophic approach .. rst-class:: emphasis This tutorial was contributed by Christine Kaufhold and `Francis - Poulin `__. + Poulin `__. The tutorial was later updated by + Emma Rothwell to add in the restrict flag in the LinearEigenproblem. As a continuation of the Quasi-Geostrophic (QG) model described in the other tutorial, we will now see how we can use Firedrake to compute the spatial @@ -138,12 +139,14 @@ We define the Test Function :math:`\phi` and the Trial Function To build the weak formulation of our equation we need to build two PETSc matrices in the form of a generalized eigenvalue problem, -:math:`A\psi = \lambda M\psi` :: +:math:`A\psi = \lambda M\psi`. This eigenproblem takes `restrict=True` to help +users to avoid convergence failures by removing eigenvalues on the +boundary, while preserving the original function space for the eigenmodes. :: eigenproblem = LinearEigenproblem( A=beta*phi*psi.dx(0)*dx, M=-inner(grad(psi), grad(phi))*dx - F*psi*phi*dx, - bcs=bc) + bcs=bc, restrict=True) Next we program our eigenvalue solver through the PETSc options system. The first is specifying that we have an generalized eigenvalue problem that is diff --git a/docs/source/conf.py b/docs/source/conf.py index fab77b03e3..19dd7dd7b5 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -264,6 +264,7 @@ \sphinxDUC{22C5}{$\cdot$} \sphinxDUC{25A3}{$\boxdot$} \sphinxDUC{03BB}{$\lambda$} +\sphinxDUC{0393}{$\Gamma$} % Sphinx equivalent of % \DeclareUnicodeCharacter{}{} diff --git a/docs/source/firedrake_24.rst b/docs/source/firedrake_24.rst index 8ab5b8e368..33265ed714 100644 --- a/docs/source/firedrake_24.rst +++ b/docs/source/firedrake_24.rst @@ -21,9 +21,76 @@ in the numerical solution of partial differential equations. The event will provide Firedrake users with the opportunity to interact directly with developers and with other users. +The conference will begin around lunchtime on 16 September and end on the +evening of 18 September. -More details about registration, submissions, accommodation and the -conference dinner will be posted in due course. +Conference venue +---------------- + +The conference will take place in the Andrew Wiles Building. Its main entrance +is located on Woodstock Road (`see here for the precise location +`__). Take some care in navigation, as +Google Maps sometimes directs people to the old mathematics building on St. +Giles. + +Accommodation +------------- + +We have arranged a number of single en-suite rooms at `Oriel college, Oxford +`__, for the nights of Monday 16 September and Tuesday +17 September. The fee includes breakfast on the following days. Oriel College, +founded in 1326, is one of the constituent colleges of the university of Oxford. +The accommodation is approximately fifteen minutes walk from the Andrew Wiles +Building. + + +Conference dinner +----------------- + +Registration includes a drinks reception and three-course dinner at `Corpus +Christi college, Oxford `__, founded in 1517. The dinner will take place on the +evening of 17 September. Corpus Christi college is right next to Oriel, about +fifteen minutes walk from the Andrew Wiles building. + +It is possible during registration to purchase additional tickets for the +dinner, to bring non-registered companions. + +.. image:: images/corpusdininghall.jpg + :width: 100% + :alt: Corpus Christi dining hall + :align: center + + +Registration +------------ + +Registration is **now open**! `To register, please follow this link +`__. + +The registration fees are as follows: + +.. list-table:: + :widths: 25 25 50 + :header-rows: 0 + + * - + - Without accommodation + - With accommodation + * - Student + - £25 + - £200 + * - Non-student + - £115 + - £350 + + +Support +------- + +The conference has been kindly supported by `EPSRC grant EP/W026163/1 +`__ and +`EPSRC grant EP/W026066/1 +`__. .. Programme diff --git a/docs/source/images/corpusdininghall.jpg b/docs/source/images/corpusdininghall.jpg new file mode 100644 index 0000000000..8d4970298d Binary files /dev/null and b/docs/source/images/corpusdininghall.jpg differ diff --git a/firedrake/adjoint_utils/variational_solver.py b/firedrake/adjoint_utils/variational_solver.py index ff16ce6b25..71538e3a22 100644 --- a/firedrake/adjoint_utils/variational_solver.py +++ b/firedrake/adjoint_utils/variational_solver.py @@ -14,15 +14,15 @@ def wrapper(self, *args, **kwargs): from firedrake import derivative, adjoint, TrialFunction init(self, *args, **kwargs) self._ad_F = self.F - self._ad_u = self.u + self._ad_u = self.u_restrict self._ad_bcs = self.bcs self._ad_J = self.J try: # Some forms (e.g. SLATE tensors) are not currently # differentiable. dFdu = derivative(self.F, - self.u, - TrialFunction(self.u.function_space())) + self.u_restrict, + TrialFunction(self.u_restrict.function_space())) self._ad_adj_F = adjoint(dFdu) except (TypeError, NotImplementedError): self._ad_adj_F = None @@ -130,7 +130,7 @@ def _ad_problem_clone(self, problem, dependencies): _ad_count_map[J_replace_map[coeff]] = coeff.count() nlvp = NonlinearVariationalProblem(replace(problem.F, F_replace_map), - F_replace_map[problem.u], + F_replace_map[problem.u_restrict], bcs=problem.bcs, J=replace(problem.J, J_replace_map)) nlvp._ad_count_map_update(_ad_count_map) diff --git a/firedrake/assemble.py b/firedrake/assemble.py index 9090377dae..5d0a4f233c 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -1631,7 +1631,7 @@ def _as_global_kernel_arg(self, tsfc_arg): # TODO Make singledispatchmethod with Python 3.8 return _as_global_kernel_arg(tsfc_arg, self) - def _get_map_arg(self, finat_element): + def _get_map_arg(self, finat_element, boundary_set): """Get the appropriate map argument for the given FInAT element. :arg finat_element: A FInAT element. @@ -1639,7 +1639,7 @@ def _get_map_arg(self, finat_element): the given FInAT element. This function uses a cache to ensure that PyOP2 knows when it can reuse maps. """ - key = self._get_map_id(finat_element) + key = self._get_map_id(finat_element), boundary_set try: return self._map_arg_cache[key] @@ -1679,20 +1679,20 @@ def _get_dim(self, finat_element): else: return (1,) - def _make_dat_global_kernel_arg(self, finat_element, index=None): + def _make_dat_global_kernel_arg(self, finat_element, boundary_set, index=None): if isinstance(finat_element, finat.EnrichedElement) and finat_element.is_mixed: assert index is None - subargs = tuple(self._make_dat_global_kernel_arg(subelem.element) + subargs = tuple(self._make_dat_global_kernel_arg(subelem.element, boundary_set) for subelem in finat_element.elements) return op2.MixedDatKernelArg(subargs) else: dim = self._get_dim(finat_element) - map_arg = self._get_map_arg(finat_element) + map_arg = self._get_map_arg(finat_element, boundary_set) return op2.DatKernelArg(dim, map_arg, index) - def _make_mat_global_kernel_arg(self, relem, celem): + def _make_mat_global_kernel_arg(self, relem, celem, rbset, cbset): if any(isinstance(e, finat.EnrichedElement) and e.is_mixed for e in {relem, celem}): - subargs = tuple(self._make_mat_global_kernel_arg(rel.element, cel.element) + subargs = tuple(self._make_mat_global_kernel_arg(rel.element, cel.element, rbset, cbset) for rel, cel in product(relem.elements, celem.elements)) shape = len(relem.elements), len(celem.elements) return op2.MixedMatKernelArg(subargs, shape) @@ -1700,7 +1700,7 @@ def _make_mat_global_kernel_arg(self, relem, celem): # PyOP2 matrix objects have scalar dims so we flatten them here rdim = numpy.prod(self._get_dim(relem), dtype=int) cdim = numpy.prod(self._get_dim(celem), dtype=int) - map_args = self._get_map_arg(relem), self._get_map_arg(celem) + map_args = self._get_map_arg(relem, rbset), self._get_map_arg(celem, cbset) return op2.MatKernelArg((((rdim, cdim),),), map_args, unroll=self._unroll) @staticmethod @@ -1737,17 +1737,21 @@ def _as_global_kernel_arg_output(_, self): if V.ufl_element().family() == "Real": return op2.GlobalKernelArg((1,)) else: - return self._make_dat_global_kernel_arg(create_element(V.ufl_element())) + return self._make_dat_global_kernel_arg(create_element(V.ufl_element()), V.boundary_set) elif rank == 2: if all(V.ufl_element().family() == "Real" for V in Vs): return op2.GlobalKernelArg((1,)) elif any(V.ufl_element().family() == "Real" for V in Vs): - el, = (create_element(V.ufl_element()) for V in Vs - if V.ufl_element().family() != "Real") - return self._make_dat_global_kernel_arg(el) + for V in Vs: + if V.ufl_element().family() != "Real": + el = create_element(V.ufl_element()) + boundary_set = V.boundary_set + break + return self._make_dat_global_kernel_arg(el, boundary_set) else: rel, cel = (create_element(V.ufl_element()) for V in Vs) - return self._make_mat_global_kernel_arg(rel, cel) + rbset, cbset = (V.boundary_set for V in Vs) + return self._make_mat_global_kernel_arg(rel, cel, rbset, cbset) else: raise AssertionError @@ -1755,7 +1759,7 @@ def _as_global_kernel_arg_output(_, self): @_as_global_kernel_arg.register(kernel_args.CoordinatesKernelArg) def _as_global_kernel_arg_coordinates(_, self): finat_element = create_element(self._mesh.ufl_coordinate_element()) - return self._make_dat_global_kernel_arg(finat_element) + return self._make_dat_global_kernel_arg(finat_element, self._mesh.coordinates.function_space().boundary_set) @_as_global_kernel_arg.register(kernel_args.CoefficientKernelArg) @@ -1773,7 +1777,7 @@ def _as_global_kernel_arg_coefficient(_, self): return op2.GlobalKernelArg((ufl_element.value_size,)) else: finat_element = create_element(ufl_element) - return self._make_dat_global_kernel_arg(finat_element, index) + return self._make_dat_global_kernel_arg(finat_element, V.boundary_set, index) @_as_global_kernel_arg.register(kernel_args.ConstantKernelArg) @@ -1788,7 +1792,7 @@ def _as_global_kernel_arg_cell_sizes(_, self): # this mirrors tsfc.kernel_interface.firedrake_loopy.KernelBuilder.set_cell_sizes ufl_element = finat.ufl.FiniteElement("P", self._mesh.ufl_cell(), 1) finat_element = create_element(ufl_element) - return self._make_dat_global_kernel_arg(finat_element) + return self._make_dat_global_kernel_arg(finat_element, self._mesh.coordinates.function_space().boundary_set) @_as_global_kernel_arg.register(kernel_args.ExteriorFacetKernelArg) @@ -1815,7 +1819,7 @@ def _as_global_kernel_arg_cell_orientations(_, self): # this mirrors firedrake.mesh.MeshGeometry.init_cell_orientations ufl_element = finat.ufl.FiniteElement("DG", cell=self._mesh.ufl_cell(), degree=0) finat_element = create_element(ufl_element) - return self._make_dat_global_kernel_arg(finat_element) + return self._make_dat_global_kernel_arg(finat_element, self._mesh.coordinates.function_space().boundary_set) @_as_global_kernel_arg.register(LayerCountKernelArg) diff --git a/firedrake/bcs.py b/firedrake/bcs.py index d3000fb8ee..61369a4083 100644 --- a/firedrake/bcs.py +++ b/firedrake/bcs.py @@ -279,6 +279,8 @@ def __init__(self, V, g, sub_domain, method=None): warnings.simplefilter('always', DeprecationWarning) warnings.warn("Selecting a bcs method is deprecated. Only topological association is supported", DeprecationWarning) + if len(V.boundary_set) and sub_domain not in V.boundary_set: + raise ValueError(f"Sub-domain {sub_domain} not in the boundary set of the restricted space.") super().__init__(V, sub_domain) if len(V) > 1: raise ValueError("Cannot apply boundary conditions on mixed spaces directly.\n" diff --git a/firedrake/checkpointing.py b/firedrake/checkpointing.py index e207b84a14..85f3f206e4 100644 --- a/firedrake/checkpointing.py +++ b/firedrake/checkpointing.py @@ -1437,7 +1437,7 @@ def _get_dm_for_checkpointing(self, tV): sd_key = self._get_shared_data_key_for_checkpointing(tV.mesh(), tV.ufl_element()) if isinstance(tV.ufl_element(), (finat.ufl.VectorElement, finat.ufl.TensorElement)): nodes_per_entity, real_tensorproduct, block_size = sd_key - global_numbering = tV.mesh().create_section(nodes_per_entity, real_tensorproduct, block_size=block_size) + global_numbering, _ = tV.mesh().create_section(nodes_per_entity, real_tensorproduct, block_size=block_size) topology_dm = tV.mesh().topology_dm dm = PETSc.DMShell().create(tV.mesh()._comm) dm.setPointSF(topology_dm.getPointSF()) diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index d3c68b5cc2..02db0f8a3b 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -1179,7 +1179,7 @@ def entity_orientations(mesh, @cython.boundscheck(False) @cython.wraparound(False) -def create_section(mesh, nodes_per_entity, on_base=False, block_size=1): +def create_section(mesh, nodes_per_entity, on_base=False, block_size=1, boundary_set=None): """Create the section describing a global numbering. :arg mesh: The mesh. @@ -1188,11 +1188,15 @@ def create_section(mesh, nodes_per_entity, on_base=False, block_size=1): extruded, the number of nodes on, and on top of, each topological entity in the base mesh. :arg on_base: If True, assume extruded space is actually Foo x Real. + :arg boundary_set: A set of boundary markers, indicating the sub-domains + a boundary condition is specified on. :arg block_size: The integer by which nodes_per_entity is uniformly multiplied to get the true data layout. :returns: A PETSc Section providing the number of dofs, and offset of each dof, on each mesh point. + :returns: An integer providing the total number of constrained nodes in the + Section. """ # We don't use DMPlexCreateSection because we only ever put one # field in each section. @@ -1200,10 +1204,13 @@ def create_section(mesh, nodes_per_entity, on_base=False, block_size=1): PETSc.DM dm PETSc.Section section PETSc.IS renumbering - PetscInt i, p, layers, pStart, pEnd + PetscInt i, p, layers, pStart, pEnd, dof, j PetscInt dimension, ndof + PetscInt *dof_array = NULL + PetscInt *entity_point_map np.ndarray[PetscInt, ndim=2, mode="c"] nodes np.ndarray[PetscInt, ndim=2, mode="c"] layer_extents + np.ndarray[PetscInt, ndim=1, mode="c"] points bint variable, extruded, on_base_ dm = mesh.topology_dm @@ -1227,12 +1234,20 @@ def create_section(mesh, nodes_per_entity, on_base=False, block_size=1): section = PETSc.Section().create(comm=mesh._comm) get_chart(dm.dm, &pStart, &pEnd) section.setChart(pStart, pEnd) - renumbering = mesh._dm_renumbering + + if boundary_set: + renumbering, (constrainedStart, constrainedEnd) = plex_renumbering(dm, + mesh._entity_classes, reordering=mesh._default_reordering, boundary_set=boundary_set) + else: + renumbering = mesh._dm_renumbering + constrainedStart = -1 + constrainedEnd = -1 + CHKERR(PetscSectionSetPermutation(section.sec, renumbering.iset)) dimension = get_topological_dimension(dm) nodes = nodes_per_entity.reshape(dimension + 1, -1) for i in range(dimension + 1): - get_depth_stratum(dm.dm, i, &pStart, &pEnd) + get_depth_stratum(dm.dm, i, &pStart, &pEnd) # gets all points at dim i if not variable: ndof = nodes[i, 0] for p in range(pStart, pEnd): @@ -1243,8 +1258,55 @@ def create_section(mesh, nodes_per_entity, on_base=False, block_size=1): layers = layer_extents[p, 1] - layer_extents[p, 0] ndof = layers*nodes[i, 0] + (layers - 1)*nodes[i, 1] CHKERR(PetscSectionSetDof(section.sec, p, block_size * ndof)) + + if boundary_set: + for marker in boundary_set: + if marker == "on_boundary": + label = "exterior_facets" + marker = 1 + else: + label = FACE_SETS_LABEL + 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(section.sec, p, &dof)) + CHKERR(PetscSectionSetConstraintDof(section.sec, p, dof)) section.setUp() - return section + + if boundary_set: + # have to loop again as we need to call section.setUp() first + for marker in boundary_set: + if marker == "on_boundary": + label = "exterior_facets" + marker = 1 + else: + label = FACE_SETS_LABEL + n = dm.getStratumSize(label, marker) + if n == 0: + continue + points = dm.getStratumIS(label, marker).indices + CHKERR(PetscSectionGetMaxDof(section.sec, &dof)) + CHKERR(PetscMalloc1(dof, &dof_array)) + for i in range(n): + p = points[i] + CHKERR(PetscSectionGetDof(section.sec, p, &dof)) + for j in range(dof): + dof_array[j] = j + CHKERR(PetscSectionSetConstraintIndices(section.sec, p, dof_array)) + CHKERR(PetscFree(dof_array)) + + constrained_nodes = 0 + + CHKERR(ISGetIndices(renumbering.iset, &entity_point_map)) + for entity in range(constrainedStart, constrainedEnd): + CHKERR(PetscSectionGetDof(section.sec, entity_point_map[entity], &dof)) + constrained_nodes += dof + CHKERR(ISRestoreIndices(renumbering.iset, &entity_point_map)) + + return section, constrained_nodes @cython.boundscheck(False) @@ -2262,7 +2324,8 @@ def validate_mesh(PETSc.DM dm): @cython.wraparound(False) def plex_renumbering(PETSc.DM plex, np.ndarray entity_classes, - np.ndarray[PetscInt, ndim=1, mode="c"] reordering=None): + np.ndarray[PetscInt, ndim=1, mode="c"] reordering=None, + boundary_set=None): """ Build a global node renumbering as a permutation of Plex points. @@ -2275,10 +2338,16 @@ def plex_renumbering(PETSc.DM plex, DMPlexGetOrdering). Optional, if not provided (or ``None``), no reordering is applied and the plex is traversed in original order. + :arg boundary_set: A set of boundary subdomains, where a DirichletBC is to + be applied. This is None if working with a FunctionSpace object, and + non-empty if working with a RestrictedFunctionSpace object. If one is + provided, this will move all core or owned boundary points to the end + of the core+owned block. The node permutation is derived from a depth-first traversal of the Plex graph over each entity class in turn. The returned IS - is the Plex -> PyOP2 permutation. + is the Plex -> PyOP2 permutation. A tuple indicating the start and end of + the core/owned constrained block is returned, for use in create_section. """ cdef: PetscInt dim, cStart, cEnd, nfacets, nclosure, c, ci, l, p, f @@ -2290,6 +2359,7 @@ def plex_renumbering(PETSc.DM plex, PETSc.IS facet_is = None PETSc.IS perm_is = None PetscBT seen = NULL + PetscBT seen_boundary = NULL PetscBool has_point DMLabel labels[3] bint reorder = reordering is not None @@ -2299,40 +2369,75 @@ def plex_renumbering(PETSc.DM plex, get_height_stratum(plex.dm, 0, &cStart, &cEnd) CHKERR(PetscMalloc1(pEnd - pStart, &perm)) CHKERR(PetscBTCreate(pEnd - pStart, &seen)) + if boundary_set: + CHKERR(PetscBTCreate(pEnd - pStart, &seen_boundary)) ncells = np.zeros(3, dtype=IntType) # Get label pointers and label-specific array indices CHKERR(DMGetLabel(plex.dm, b"pyop2_core", &labels[0])) CHKERR(DMGetLabel(plex.dm, b"pyop2_owned", &labels[1])) CHKERR(DMGetLabel(plex.dm, b"pyop2_ghost", &labels[2])) - for l in range(3): - CHKERR(DMLabelCreateIndex(labels[l], pStart, pEnd)) + for idx in range(3): + CHKERR(DMLabelCreateIndex(labels[idx], pStart, pEnd)) entity_classes = entity_classes.astype(IntType) - lidx = np.zeros(3, dtype=IntType) - lidx[1] = sum(entity_classes[:, 0]) - lidx[2] = sum(entity_classes[:, 1]) + # Get boundary points (if the boundary_set exists) and count each type + constrained_core = 0 + constrained_owned = 0 + if boundary_set: + for marker in boundary_set: + if marker == "on_boundary": + label = "exterior_facets" + marker = 1 + else: + label = FACE_SETS_LABEL + n = plex.getStratumSize(label, marker) + if n == 0: + continue + points = plex.getStratumIS(label, marker).indices + for i in range(n): + p = points[i] + if not PetscBTLookup(seen_boundary, p): + for idx in range(3): + CHKERR(DMLabelHasPoint(labels[idx], p, &has_point)) + if has_point: + PetscBTSet(seen_boundary, p) + if idx == 1: + constrained_owned += 1 + elif idx == 0: + constrained_core += 1 + break + + # assign lists + lidx = np.zeros(4, dtype=IntType) + lidx[1] = sum(entity_classes[:, 0]) - constrained_core + lidx[2] = sum(entity_classes[:, 1]) + lidx[3] = lidx[2] - constrained_core - constrained_owned + for c in range(pStart, pEnd): if reorder: cell = reordering[c] else: cell = c - # We always re-order cell-wise so that we inherit any cache # coherency from the reordering provided by the Plex if cStart <= cell < cEnd: - # Get cell closure get_transitive_closure(plex.dm, cell, PETSC_TRUE, &nclosure, &closure) for ci in range(nclosure): p = closure[2*ci] if not PetscBTLookup(seen, p): - for l in range(3): - CHKERR(DMLabelHasPoint(labels[l], p, &has_point)) + for idx in range(3): + CHKERR(DMLabelHasPoint(labels[idx], p, &has_point)) if has_point: PetscBTSet(seen, p) - perm[lidx[l]] = p - lidx[l] += 1 + if boundary_set and PetscBTLookup(seen_boundary, p) and idx <= 1: + # push boundary point to end of constrained owned dofs + perm[lidx[3]] = p + lidx[3] += 1 + else: + perm[lidx[idx]] = p + lidx[idx] += 1 break if closure != NULL: @@ -2341,11 +2446,15 @@ def plex_renumbering(PETSc.DM plex, CHKERR(DMLabelDestroyIndex(labels[c])) CHKERR(PetscBTDestroy(&seen)) + + if boundary_set: + CHKERR(PetscBTDestroy(&seen_boundary)) + perm_is = PETSc.IS().create(comm=plex.comm) perm_is.setType("general") CHKERR(ISGeneralSetIndices(perm_is.iset, pEnd - pStart, perm, PETSC_OWN_POINTER)) - return perm_is + return perm_is, (lidx[1], lidx[3]) @cython.boundscheck(False) @cython.wraparound(False) @@ -3194,7 +3303,7 @@ def make_global_numbering(PETSc.Section lsec, PETSc.Section gsec): :arg lsec: Section describing local dof layout and numbers. :arg gsec: Section describing global dof layout and numbers.""" cdef: - PetscInt c, p, pStart, pEnd, dof, loff, goff + PetscInt c, p, pStart, pEnd, dof, cdof, loff, goff np.ndarray[PetscInt, ndim=1, mode="c"] val val = np.empty(lsec.getStorageSize(), dtype=IntType) @@ -3202,12 +3311,17 @@ def make_global_numbering(PETSc.Section lsec, PETSc.Section gsec): for p in range(pStart, pEnd): CHKERR(PetscSectionGetDof(lsec.sec, p, &dof)) + CHKERR(PetscSectionGetConstraintDof(lsec.sec, p, &cdof)) if dof > 0: CHKERR(PetscSectionGetOffset(lsec.sec, p, &loff)) CHKERR(PetscSectionGetOffset(gsec.sec, p, &goff)) - goff = cabs(goff) - for c in range(dof): - val[loff + c] = goff + c + if cdof > 0: + for c in range(dof): + val[loff + c] = -1 + else: + goff = cabs(goff) + for c in range(dof): + val[loff + c] = goff + c return val diff --git a/firedrake/cython/petschdr.pxi b/firedrake/cython/petschdr.pxi index b64a5620f5..65d68572bd 100644 --- a/firedrake/cython/petschdr.pxi +++ b/firedrake/cython/petschdr.pxi @@ -90,6 +90,9 @@ cdef extern from "petscis.h" nogil: int PetscSectionSetDof(PETSc.PetscSection,PetscInt,PetscInt) int PetscSectionSetFieldDof(PETSc.PetscSection,PetscInt,PetscInt,PetscInt) int PetscSectionGetFieldDof(PETSc.PetscSection,PetscInt,PetscInt,PetscInt*) + int PetscSectionGetConstraintDof(PETSc.PetscSection,PetscInt,PetscInt*) + int PetscSectionSetConstraintDof(PETSc.PetscSection,PetscInt,PetscInt) + int PetscSectionSetConstraintIndices(PETSc.PetscSection,PetscInt, PetscInt[]) int PetscSectionGetMaxDof(PETSc.PetscSection,PetscInt*) int PetscSectionSetPermutation(PETSc.PetscSection,PETSc.PetscIS) int ISGetIndices(PETSc.PetscIS,PetscInt*[]) diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py index c883bf96be..ae8a9f020b 100644 --- a/firedrake/eigensolver.py +++ b/firedrake/eigensolver.py @@ -1,9 +1,13 @@ """Specify and solve finite element eigenproblems.""" from firedrake.assemble import assemble +from firedrake.bcs import DirichletBC from firedrake.function import Function +from firedrake.functionspace import RestrictedFunctionSpace +from firedrake.ufl_expr import TrialFunction, TestFunction from firedrake import utils from firedrake.petsc import OptionsManager, flatten_parameters from firedrake.exceptions import ConvergenceError +from ufl import replace, inner, dx try: from slepc4py import SLEPc except ImportError: @@ -28,35 +32,56 @@ class LinearEigenproblem(): bcs : DirichletBC or list of DirichletBC The boundary conditions. bc_shift: float - The value to shift the boundary condition eigenvalues by. + The value to shift the boundary condition eigenvalues by. This value + will be ignored if restrict==True. + restrict: bool + If True, replace the function spaces of u and v with their restricted + version. The output space remains unchanged. Notes ----- - - If Dirichlet boundary conditions are supplied then these will result in the - eigenproblem having a nullspace spanned by the basis functions with support - on the boundary. To facilitate solution, this is shifted by the specified - amount. It is the user's responsibility to ensure that the shift is not - close to an actual eigenvalue of the system. + If restrict==True, the arguments of A and M will be replaced, such that + their function space is replaced by the equivalent RestrictedFunctionSpace + class. This avoids the computation of eigenvalues associated with the + Dirichlet boundary conditions. This in turn prevents convergence failures, + and allows only the non-boundary eigenvalues to be returned. The + eigenvectors will be in the original, non-restricted space. + + If restrict==False and Dirichlet boundary conditions are supplied, then + these conditions will result in the eigenproblem having a nullspace spanned + by the basis functions with support on the boundary. To facilitate + solution, this is shifted by the specified amount. It is the user's + responsibility to ensure that the shift is not close to an actual + eigenvalue of the system. """ - def __init__(self, A, M=None, bcs=None, bc_shift=0.0): + def __init__(self, A, M=None, bcs=None, bc_shift=0.0, restrict=True): if not SLEPc: raise ImportError( "Unable to import SLEPc, eigenvalue computation not possible " "(try firedrake-update --slepc)" ) - self.A = A # LHS args = A.arguments() v, u = args - if M: - self.M = M - else: - from ufl import inner, dx - self.M = inner(u, v) * dx self.output_space = u.function_space() - self.bcs = bcs self.bc_shift = bc_shift + self.restrict = restrict + + if not M: + M = inner(u, v) * dx + + if restrict and bcs: # assumed u and v are in the same space here + V_res = RestrictedFunctionSpace(self.output_space, boundary_set=set([bc.sub_domain for bc in bcs])) + u_res = TrialFunction(V_res) + v_res = TestFunction(V_res) + self.M = replace(M, {u: u_res, v: v_res}) + self.A = replace(A, {u: u_res, v: v_res}) + self.bcs = [DirichletBC(V_res, bc.function_arg, bc.sub_domain) for bc in bcs] + self.restricted_space = V_res + else: + self.A = A # LHS + self.M = M + self.bcs = bcs def dirichlet_bcs(self): """Return an iterator over the Dirichlet boundary conditions.""" @@ -66,7 +91,10 @@ def dirichlet_bcs(self): @utils.cached_property def dm(self): r"""Return the dm associated with the output space.""" - return self.output_space.dm + if self.restrict: + return self.restricted_space.dm + else: + return self.output_space.dm class LinearEigensolver(OptionsManager): @@ -186,9 +214,16 @@ def eigenfunction(self, i): (Function, Function) The real and imaginary parts of the eigenfunction. """ - eigenmodes_real = Function(self._problem.output_space) # fn of V - eigenmodes_imag = Function(self._problem.output_space) + if self._problem.restrict: + eigenmodes_real = Function(self._problem.restricted_space) + eigenmodes_imag = Function(self._problem.restricted_space) + else: + eigenmodes_real = Function(self._problem.output_space) # fn of V + eigenmodes_imag = Function(self._problem.output_space) with eigenmodes_real.dat.vec_wo as vr: with eigenmodes_imag.dat.vec_wo as vi: self.es.getEigenvector(i, vr, vi) # gets the i-th eigenvector + if self._problem.restrict: + eigenmodes_real = Function(self._problem.output_space).interpolate(eigenmodes_real) + eigenmodes_imag = Function(self._problem.output_space).interpolate(eigenmodes_imag) return eigenmodes_real, eigenmodes_imag # returns Firedrake fns diff --git a/firedrake/functionspace.py b/firedrake/functionspace.py index 8345a0eb7e..2d6726b2a1 100644 --- a/firedrake/functionspace.py +++ b/firedrake/functionspace.py @@ -16,7 +16,7 @@ __all__ = ("MixedFunctionSpace", "FunctionSpace", - "VectorFunctionSpace", "TensorFunctionSpace") + "VectorFunctionSpace", "TensorFunctionSpace", "RestrictedFunctionSpace") @PETSc.Log.EventDecorator() @@ -292,9 +292,9 @@ def rec(eles): spaces = tuple(s.topological for s in flatten(spaces)) # Error checking for space in spaces: - if type(space) in (impl.FunctionSpace, impl.RealFunctionSpace): + if type(space) in (impl.FunctionSpace, impl.RealFunctionSpace, impl.RestrictedFunctionSpace): continue - elif type(space) is impl.ProxyFunctionSpace: + elif type(space) in (impl.ProxyFunctionSpace, impl.ProxyRestrictedFunctionSpace): if space.component is not None: raise ValueError("Can't make mixed space with %s" % space) continue @@ -305,3 +305,23 @@ def rec(eles): if mesh is not mesh.topology: new = cls.create(new, mesh) return new + + +@PETSc.Log.EventDecorator("CreateFunctionSpace") +def RestrictedFunctionSpace(function_space, name=None, boundary_set=[]): + """Create a :class:`.RestrictedFunctionSpace`. + + Parameters + ---------- + function_space : + FunctionSpace object to restrict + name : + An optional name for the function space. + boundary_set : + A set of subdomains of the mesh in which Dirichlet boundary conditions + will be applied. + + """ + return impl.WithGeometry.create(impl.RestrictedFunctionSpace(function_space, name=name, + boundary_set=boundary_set), + function_space.mesh()) diff --git a/firedrake/functionspacedata.py b/firedrake/functionspacedata.py index e9bf789e41..be00e3ffff 100644 --- a/firedrake/functionspacedata.py +++ b/firedrake/functionspacedata.py @@ -75,16 +75,17 @@ def get_global_numbering(mesh, key, global_numbering=None): entities. :arg mesh: The mesh to use. - :arg key: a (nodes_per_entity, real_tensorproduct) tuple where + :arg key: a (nodes_per_entity, real_tensorproduct, boundary_set) tuple where nodes_per_entity is a tuple of the number of nodes per topological entity; real_tensorproduct is True if the function space is a - degenerate fs x Real tensorproduct. + degenerate fs x Real tensorproduct; boundary_set is a set of boundary + markers, indicating sub-domains a boundary condition is specified on. :returns: A new PETSc Section. """ if global_numbering: return global_numbering - nodes_per_entity, real_tensorproduct = key - return mesh.create_section(nodes_per_entity, real_tensorproduct) + nodes_per_entity, real_tensorproduct, boundary_set = key + return mesh.create_section(nodes_per_entity, real_tensorproduct, boundary_set=boundary_set) @cached @@ -92,17 +93,18 @@ def get_node_set(mesh, key): """Get the :class:`node set `. :arg mesh: The mesh to use. - :arg key: a (nodes_per_entity, real_tensorproduct) tuple where - nodes_per_entity is a tuple of the number of nodes per topological - entity; real_tensorproduct is True if the function space is a - degenerate fs x Real tensorproduct. + :arg key: a (nodes_per_entity, real_tensorproduct, boundary_set) tuple + where nodes_per_entity is a tuple of the number of nodes per + topological entity; real_tensorproduct is True if the function space is + a degenerate fs x Real tensorproduct; boundary_set is a set of boundary + markers, indicating sub-domains a boundary condition is specified on. :returns: A :class:`pyop2.Set` for the function space nodes. """ - nodes_per_entity, real_tensorproduct = key - global_numbering = get_global_numbering(mesh, (nodes_per_entity, real_tensorproduct)) + nodes_per_entity, real_tensorproduct, _ = key + global_numbering, constrained_size = get_global_numbering(mesh, key) node_classes = mesh.node_classes(nodes_per_entity, real_tensorproduct=real_tensorproduct) halo = halo_mod.Halo(mesh.topology_dm, global_numbering, comm=mesh.comm) - node_set = op2.Set(node_classes, halo=halo, comm=mesh.comm) + node_set = op2.Set(node_classes, halo=halo, comm=mesh.comm, constrained_size=constrained_size) extruded = mesh.cell_set._extruded assert global_numbering.getStorageSize() == node_set.total_size @@ -150,7 +152,8 @@ def get_entity_node_lists(mesh, key, entity_dofs, entity_permutations, global_nu """Get the map from mesh entity sets to function space nodes. :arg mesh: The mesh to use. - :arg key: a (entity_dofs_key, real_tensorproduct, entity_permutations_key) tuple. + :arg key: a (entity_dofs_key, real_tensorproduct, entity_permutations_key, + boundary_set) tuple. :arg entity_dofs: FInAT entity dofs. :arg entity_permutations: FInAT entity permutations. :arg global_numbering: The PETSc Section describing node layout @@ -183,10 +186,12 @@ def get_map_cache(mesh, key): """Get the map cache for this mesh. :arg mesh: The mesh to use. - :arg key: a (entity_dofs_key, real_tensorproduct, entity_permutations_key) tuple where - entity_dofs is Canonicalised entity_dofs (see :func:`entity_dofs_key`); - real_tensorproduct is True if the function space is a degenerate - fs x Real tensorproduct. + :arg key: a (entity_dofs_key, real_tensorproduct, entity_permutations_key, + boundary_set) tuple where entity_dofs is Canonicalised entity_dofs + (see :func:`entity_dofs_key`); real_tensorproduct is True if the + function space is a degenerate fs x Real tensorproduct; boundary_set is + the set of subdomains a restricted function space is applied to, or + None if using a regular function space. """ if type(mesh.topology) is mesh_mod.VertexOnlyMeshTopology: return {mesh.cell_set: None} @@ -268,13 +273,13 @@ def get_top_bottom_boundary_nodes(mesh, key, V): """Get top or bottom boundary nodes of an extruded function space. :arg mesh: The mesh to cache on. - :arg key: The key a 2-tuple of ``(entity_dofs_key, sub_domain)``. + :arg key: A 3-tuple of ``(entity_dofs_key, sub_domain, boundary_set)`` key. Where sub_domain indicates top or bottom. :arg V: The FunctionSpace to select from. :arg entity_dofs: The flattened entity dofs. :returnsL: A numpy array of the (unique) boundary nodes. """ - _, sub_domain = key + _, sub_domain, boundary_set = key cell_node_list = V.cell_node_list offset = V.offset if mesh.variable_layers: @@ -302,11 +307,11 @@ 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 key: (edofs, sub_domain, boundary_set) 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 + _, sub_domain, boundary_set = key if sub_domain not in {"on_boundary", "top", "bottom"}: valid = set(mesh.interior_facets.unique_markers) valid |= set(mesh.exterior_facets.unique_markers) @@ -396,17 +401,22 @@ class FunctionSpaceData(object): :arg mesh: The mesh to share the data on. :arg ufl_element: The UFL element. + :arg boundary_set: The set of subdomains that a Dirichlet boundary condition + will act on. This is None if the function space is not a + :class:`.RestrictedFunctionSpace`. """ __slots__ = ("real_tensorproduct", "map_cache", "entity_node_lists", "node_set", "cell_boundary_masks", "interior_facet_boundary_masks", "offset", "offset_quotient", - "extruded", "mesh", "global_numbering") + "extruded", "mesh", "global_numbering", "boundary_set") @PETSc.Log.EventDecorator() - def __init__(self, mesh, ufl_element): + def __init__(self, mesh, ufl_element, boundary_set=None): if type(ufl_element) is finat.ufl.MixedElement: raise ValueError("Can't create FunctionSpace for MixedElement") + self.boundary_set = boundary_set + finat_element = create_element(ufl_element) real_tensorproduct = eutils.is_real_tensor_product_element(finat_element) entity_dofs = finat_element.entity_dofs() @@ -418,9 +428,9 @@ def __init__(self, mesh, ufl_element): # Create the PetscSection mapping topological entities to functionspace nodes # For non-scalar valued function spaces, there are multiple dofs per node. - key = (nodes_per_entity, real_tensorproduct) + key = (nodes_per_entity, real_tensorproduct, boundary_set) # These are keyed only on nodes per topological entity. - global_numbering = get_global_numbering(mesh, key) + global_numbering, constrained_size = get_global_numbering(mesh, key) node_set = get_node_set(mesh, key) edofs_key = entity_dofs_key(entity_dofs) @@ -432,7 +442,7 @@ def __init__(self, mesh, ufl_element): # implementation because of the need to support boundary # conditions. # Map caches are specific to a cell_node_list, which is keyed by entity_dof - self.map_cache = get_map_cache(mesh, (edofs_key, real_tensorproduct, eperm_key)) + self.map_cache = get_map_cache(mesh, (edofs_key, real_tensorproduct, eperm_key, boundary_set)) if isinstance(mesh, mesh_mod.ExtrudedMeshTopology): self.offset = eutils.calculate_dof_offset(finat_element) @@ -443,7 +453,7 @@ def __init__(self, mesh, ufl_element): else: self.offset_quotient = None - self.entity_node_lists = get_entity_node_lists(mesh, (edofs_key, real_tensorproduct, eperm_key), entity_dofs, entity_permutations, global_numbering, self.offset) + self.entity_node_lists = get_entity_node_lists(mesh, (edofs_key, real_tensorproduct, eperm_key, boundary_set), entity_dofs, entity_permutations, global_numbering, self.offset) self.node_set = node_set self.cell_boundary_masks = get_boundary_masks(mesh, (edofs_key, "cell"), finat_element) self.interior_facet_boundary_masks = get_boundary_masks(mesh, (edofs_key, "interior_facet"), finat_element) @@ -473,14 +483,14 @@ def boundary_nodes(self, V, sub_domain): raise ValueError("Invalid subdomain '%s' for non-extruded mesh", sub_domain) entity_dofs = eutils.flat_entity_dofs(V.finat_element.entity_dofs()) - key = (entity_dofs_key(entity_dofs), sub_domain) + key = (entity_dofs_key(entity_dofs), sub_domain, V.boundary_set) return get_top_bottom_boundary_nodes(V.mesh(), key, V) else: if sub_domain == "on_boundary": sdkey = sub_domain else: sdkey = as_tuple(sub_domain) - key = (entity_dofs_key(V.finat_element.entity_dofs()), sdkey) + key = (entity_dofs_key(V.finat_element.entity_dofs()), sdkey, V.boundary_set) return get_facet_closure_nodes(V.mesh(), key, V) @PETSc.Log.EventDecorator() @@ -511,12 +521,14 @@ def get_map(self, V, entity_set, map_arity, name, offset, offset_quotient): @PETSc.Log.EventDecorator() -def get_shared_data(mesh, ufl_element): +def get_shared_data(mesh, ufl_element, boundary_set=None): """Return the ``FunctionSpaceData`` for the given element. :arg mesh: The mesh to build the function space data on. :arg ufl_element: A UFL element. + :arg boundary_set: A set of boundary markers, indicating the subdomains a + boundary condition is specified on. :raises ValueError: if mesh or ufl_element are invalid. :returns: a ``FunctionSpaceData`` object with the shared data. @@ -526,4 +538,4 @@ def get_shared_data(mesh, ufl_element): if not isinstance(ufl_element, finat.ufl.finiteelement.FiniteElementBase): raise ValueError("Can't create function space data from a %s" % type(ufl_element)) - return FunctionSpaceData(mesh, ufl_element) + return FunctionSpaceData(mesh, ufl_element, boundary_set) diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py index a2ece362d8..5308cf4cdf 100644 --- a/firedrake/functionspaceimpl.py +++ b/firedrake/functionspaceimpl.py @@ -95,7 +95,7 @@ def __init__(self, mesh, element, component=None, cargo=None): assert component is None or isinstance(component, int) assert cargo is None or isinstance(cargo, FunctionSpaceCargo) - super().__init__(mesh, element) + super().__init__(mesh, element, label=cargo.topological._label or "") self.component = component self.cargo = cargo self.comm = mesh.comm @@ -466,6 +466,9 @@ class FunctionSpace(object): which provides extra error checking and argument sanitising. """ + + boundary_set = frozenset() + @PETSc.Log.EventDecorator() def __init__(self, mesh, element, name=None): super(FunctionSpace, self).__init__() @@ -490,7 +493,8 @@ def __init__(self, mesh, element, name=None): self.shape = rvs[:len(rvs) - len(sub)] else: self.shape = () - self._ufl_function_space = ufl.FunctionSpace(mesh.ufl_mesh(), element) + self._label = "" + self._ufl_function_space = ufl.FunctionSpace(mesh.ufl_mesh(), element, label=self._label) self._mesh = mesh self.rank = len(self.shape) @@ -615,6 +619,9 @@ def ufl_function_space(self): r"""The :class:`~ufl.classes.FunctionSpace` associated with this space.""" return self._ufl_function_space + def label(self): + return self._label + def __len__(self): return 1 @@ -627,9 +634,7 @@ def __repr__(self): self.name) def __str__(self): - return "FunctionSpace(%s, %s, name=%s)" % (self.mesh(), - self.ufl_element(), - self.name) + return self.__repr__() @utils.cached_property def subfunctions(self): @@ -670,7 +675,10 @@ def node_count(self): this process. If the :class:`FunctionSpace` has :attr:`FunctionSpace.rank` 0, this is equal to the :attr:`FunctionSpace.dof_count`, otherwise the :attr:`FunctionSpace.dof_count` is :attr:`dim` times the :attr:`node_count`.""" - return self.node_set.total_size + constrained_node_set = set() + for sub_domain in self.boundary_set: + constrained_node_set.update(self._shared_data.boundary_nodes(self, sub_domain)) + return self.node_set.total_size - len(constrained_node_set) @utils.cached_property def dof_count(self): @@ -792,6 +800,90 @@ def collapse(self): return FunctionSpace(self.mesh(), self.ufl_element()) +class RestrictedFunctionSpace(FunctionSpace): + r"""A representation of a function space, with additional information + about where boundary conditions are to be applied. + + If a :class:`FunctionSpace` is represented as V, we can decompose V into + V = V0 + VΓ, where V0 contains functions in the basis of V that vanish on + the boundary where a boundary condition is applied, and VΓ contains all + other basis functions. The :class:`RestrictedFunctionSpace` + corresponding to V takes functions only from V0 when solving problems, or + when creating a TestFunction and TrialFunction. The values on the boundary + set will remain constant when solving, but are present in the + output of the solver. + + :arg function_space: The :class:`FunctionSpace` to restrict. + :kwarg name: An optional name for this :class:`RestrictedFunctionSpace`, + useful for later identification. + :kwarg boundary_set: A set of subdomains on which a DirichletBC will be + applied. + + Notes + ----- + If using this class to solve or similar, a list of DirichletBCs will still + need to be specified on this space and passed into the function. + """ + def __init__(self, function_space, name=None, boundary_set=frozenset()): + label = "" + for boundary_domain in boundary_set: + label += str(boundary_domain) + label += "_" + self.boundary_set = frozenset(boundary_set) + super().__init__(function_space._mesh.topology, + function_space.ufl_element(), function_space.name) + self._label = label + self._ufl_function_space = ufl.FunctionSpace(function_space._mesh.ufl_mesh(), + function_space.ufl_element(), + label=self._label) + self.function_space = function_space + self.name = name or (function_space.name or "Restricted" + "_" + + "_".join(sorted( + [str(i) for i in self.boundary_set]))) + + def set_shared_data(self): + sdata = get_shared_data(self._mesh, self.ufl_element(), self.boundary_set) + self._shared_data = sdata + self.node_set = sdata.node_set + r"""A :class:`pyop2.types.set.Set` representing the function space nodes.""" + self.dof_dset = op2.DataSet(self.node_set, self.shape or 1, + name="%s_nodes_dset" % self.name) + r"""A :class:`pyop2.types.dataset.DataSet` representing the function space + degrees of freedom.""" + + # check not all degrees of freedom are constrained + unconstrained_dofs = self.dof_dset.size - self.dof_dset.constrained_size + if self.comm.allreduce(unconstrained_dofs) == 0: + raise ValueError("All degrees of freedom are constrained.") + self.finat_element = create_element(self.ufl_element()) + # Used for reconstruction of mixed/component spaces. + # sdata carries real_tensorproduct. + self.real_tensorproduct = sdata.real_tensorproduct + self.extruded = sdata.extruded + self.offset = sdata.offset + self.offset_quotient = sdata.offset_quotient + self.cell_boundary_masks = sdata.cell_boundary_masks + self.interior_facet_boundary_masks = sdata.interior_facet_boundary_masks + self.global_numbering = sdata.global_numbering + + def __eq__(self, other): + if not isinstance(other, RestrictedFunctionSpace): + return False + return self.function_space == other.function_space and \ + self.boundary_set == other.boundary_set + + def __repr__(self): + return self.__class__.__name__ + "(%r, name=%r, boundary_set=%r)" % ( + str(self.function_space), self.name, self.boundary_set) + + def __hash__(self): + return hash((self.mesh(), self.dof_dset, self.ufl_element(), + self.boundary_set)) + + def local_to_global_map(self, bcs, lgmap=None): + return lgmap or self.dof_dset.lgmap + + class MixedFunctionSpace(object): r"""A function space on a mixed finite element. @@ -815,6 +907,11 @@ def __init__(self, spaces, name=None): self._ufl_function_space = ufl.FunctionSpace(mesh.ufl_mesh(), finat.ufl.MixedElement(*[s.ufl_element() for s in spaces])) self.name = name or "_".join(str(s.name) for s in spaces) + label = "" + for s in spaces: + label += "(" + s._label + ")_" + self._label = label + self.boundary_set = frozenset() self._subspaces = {} self._mesh = mesh self.comm = mesh.comm @@ -1048,6 +1145,55 @@ def make_dat(self, *args, **kwargs): return super(ProxyFunctionSpace, self).make_dat(*args, **kwargs) +class ProxyRestrictedFunctionSpace(RestrictedFunctionSpace): + r"""A :class:`RestrictedFunctionSpace` that one can attach extra properties to. + + :arg function_space: The function space to be restricted. + :kwarg name: The name of the restricted function space. + :kwarg boundary_set: The boundary domains on which boundary conditions will + be specified + + .. warning:: + + Users should not build a :class:`ProxyRestrictedFunctionSpace` directly, + it is mostly used as an internal implementation detail. + """ + def __new__(cls, function_space, name=None, boundary_set=frozenset()): + topology = function_space._mesh.topology + self = super(ProxyRestrictedFunctionSpace, cls).__new__(cls) + if function_space._mesh is not topology: + return WithGeometry.create(self, function_space._mesh) + else: + return self + + def __repr__(self): + return "%sProxyRestrictedFunctionSpace(%r, name=%r, boundary_set=%r, index=%r, component=%r)" % \ + (str(self.identifier).capitalize(), + str(self.function_space), + self.name, + self.boundary_set, + self.index, + self.component) + + def __str__(self): + return self.__repr__() + + identifier = None + r"""An optional identifier, for debugging purposes.""" + + no_dats = False + r"""Can this proxy make :class:`pyop2.types.dat.Dat` objects""" + + def make_dat(self, *args, **kwargs): + r"""Create a :class:`pyop2.types.dat.Dat`. + + :raises ValueError: if :attr:`no_dats` is ``True``. + """ + if self.no_dats: + raise ValueError("Can't build Function on %s function space" % self.identifier) + return super(ProxyRestrictedFunctionSpace, self).make_dat(*args, **kwargs) + + def IndexedFunctionSpace(index, space, parent): r"""Build a new FunctionSpace that remembers it is a particular subspace of a :class:`MixedFunctionSpace`. @@ -1060,6 +1206,8 @@ def IndexedFunctionSpace(index, space, parent): """ if space.ufl_element().family() == "Real": new = RealFunctionSpace(space.mesh(), space.ufl_element(), name=space.name) + elif len(space.boundary_set) > 0: + new = ProxyRestrictedFunctionSpace(space.function_space, name=space.name, boundary_set=space.boundary_set) else: new = ProxyFunctionSpace(space.mesh(), space.ufl_element(), name=space.name) new.index = index diff --git a/firedrake/mesh.py b/firedrake/mesh.py index a09092fe62..1b6c0e2175 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -568,16 +568,16 @@ def callback(self): tdim = dmcommon.get_topological_dimension(self.topology_dm) entity_dofs = np.zeros(tdim+1, dtype=IntType) entity_dofs[-1] = 1 - self._cell_numbering = self.create_section(entity_dofs) + self._cell_numbering, _ = self.create_section(entity_dofs) if tdim == 0: self._vertex_numbering = self._cell_numbering else: entity_dofs[:] = 0 entity_dofs[0] = 1 - self._vertex_numbering = self.create_section(entity_dofs) + self._vertex_numbering, _ = self.create_section(entity_dofs) entity_dofs[:] = 0 entity_dofs[-2] = 1 - facet_numbering = self.create_section(entity_dofs) + facet_numbering, _ = self.create_section(entity_dofs) self._facet_ordering = dmcommon.get_facet_ordering(self.topology_dm, facet_numbering) self._callback = callback self.name = name @@ -737,16 +737,18 @@ def cell_to_facets(self): """ pass - def create_section(self, nodes_per_entity, real_tensorproduct=False, block_size=1): + def create_section(self, nodes_per_entity, real_tensorproduct=False, block_size=1, boundary_set=None): """Create a PETSc Section describing a function space. :arg nodes_per_entity: number of function space nodes per topological entity. :arg real_tensorproduct: If True, assume extruded space is actually Foo x Real. :arg block_size: The integer by which nodes_per_entity is uniformly multiplied to get the true data layout. + :arg boundary_set: A set of boundary markers, indicating the subdomains + a boundary condition is specified on. :returns: a new PETSc Section. """ - return dmcommon.create_section(self, nodes_per_entity, on_base=real_tensorproduct, block_size=block_size) + return dmcommon.create_section(self, nodes_per_entity, on_base=real_tensorproduct, block_size=block_size, boundary_set=boundary_set) def node_classes(self, nodes_per_entity, real_tensorproduct=False): """Compute node classes given nodes per entity. @@ -1063,7 +1065,7 @@ def _renumber_entities(self, reorder): else: # No reordering reordering = None - return dmcommon.plex_renumbering(self.topology_dm, self._entity_classes, reordering) + return dmcommon.plex_renumbering(self.topology_dm, self._entity_classes, reordering)[0] @utils.cached_property def cell_closure(self): @@ -1617,7 +1619,7 @@ def _renumber_entities(self, reorder): perm_is.setIndices(perm) return perm_is else: - return dmcommon.plex_renumbering(self.topology_dm, self._entity_classes, None) + return dmcommon.plex_renumbering(self.topology_dm, self._entity_classes, None)[0] @utils.cached_property # TODO: Recalculate if mesh moves def cell_closure(self): diff --git a/firedrake/solving.py b/firedrake/solving.py index b5ece29988..a5447681c1 100644 --- a/firedrake/solving.py +++ b/firedrake/solving.py @@ -120,6 +120,10 @@ def solve(*args, **kwargs): In the same fashion you can add the near nullspace using the ``near_nullspace`` keyword argument. + + To exclude Dirichlet boundary condition nodes through the use of a + :class`.RestrictedFunctionSpace`, set the ``restrict`` keyword + argument to be True. """ assert len(args) > 0 @@ -138,7 +142,7 @@ def _solve_varproblem(*args, **kwargs): eq, u, bcs, J, Jp, M, form_compiler_parameters, \ solver_parameters, nullspace, nullspace_T, \ near_nullspace, \ - options_prefix = _extract_args(*args, **kwargs) + options_prefix, restrict = _extract_args(*args, **kwargs) # Check whether solution is valid if not isinstance(u, (function.Function, vector.Vector)): @@ -153,7 +157,8 @@ def _solve_varproblem(*args, **kwargs): if isinstance(eq.lhs, ufl.Form) and isinstance(eq.rhs, ufl.BaseForm): # Create problem problem = vs.LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs, Jp, - form_compiler_parameters=form_compiler_parameters) + form_compiler_parameters=form_compiler_parameters, + restrict=restrict) # Create solver and call solve solver = vs.LinearVariationalSolver(problem, solver_parameters=solver_parameters, nullspace=nullspace, @@ -169,7 +174,8 @@ def _solve_varproblem(*args, **kwargs): raise TypeError("Only '0' support on RHS of nonlinear Equation, not %r" % eq.rhs) # Create problem problem = vs.NonlinearVariationalProblem(eq.lhs, u, bcs, J, Jp, - form_compiler_parameters=form_compiler_parameters) + form_compiler_parameters=form_compiler_parameters, + restrict=restrict) # Create solver and call solve solver = vs.NonlinearVariationalSolver(problem, solver_parameters=solver_parameters, nullspace=nullspace, @@ -278,7 +284,7 @@ def _extract_args(*args, **kwargs): valid_kwargs = ["bcs", "J", "Jp", "M", "form_compiler_parameters", "solver_parameters", "nullspace", "transpose_nullspace", "near_nullspace", - "options_prefix", "appctx"] + "options_prefix", "appctx", "restrict"] for kwarg in kwargs.keys(): if kwarg not in valid_kwargs: raise RuntimeError("Illegal keyword argument '%s'; valid keywords \ @@ -321,9 +327,11 @@ def _extract_args(*args, **kwargs): form_compiler_parameters = kwargs.get("form_compiler_parameters", {}) solver_parameters = kwargs.get("solver_parameters", {}) options_prefix = kwargs.get("options_prefix", None) + restrict = kwargs.get("restrict", False) return eq, u, bcs, J, Jp, M, form_compiler_parameters, \ - solver_parameters, nullspace, nullspace_T, near_nullspace, options_prefix + solver_parameters, nullspace, nullspace_T, near_nullspace, \ + options_prefix, restrict def _extract_bcs(bcs): diff --git a/firedrake/solving_utils.py b/firedrake/solving_utils.py index 0d62c9be13..01d1bf4a43 100644 --- a/firedrake/solving_utils.py +++ b/firedrake/solving_utils.py @@ -197,7 +197,7 @@ def __init__(self, problem, mat_type, pmat_type, appctx=None, self.fcp = problem.form_compiler_parameters # Function to hold current guess - self._x = problem.u + self._x = problem.u_restrict if appctx is None: appctx = {} @@ -327,7 +327,7 @@ def split(self, fields): for field in fields: F = splitter.split(problem.F, argument_indices=(field, )) J = splitter.split(problem.J, argument_indices=(field, field)) - us = problem.u.subfunctions + us = problem.u_restrict.subfunctions V = F.arguments()[0].function_space() # Exposition: # We are going to make a new solution Function on the sub @@ -371,11 +371,11 @@ def split(self, fields): # solving for, and some spaces that have just become # coefficients in the new form. u = as_vector(vec) - F = replace(F, {problem.u: u}) - J = replace(J, {problem.u: u}) + F = replace(F, {problem.u_restrict: u}) + J = replace(J, {problem.u_restrict: u}) if problem.Jp is not None: Jp = splitter.split(problem.Jp, argument_indices=(field, field)) - Jp = replace(Jp, {problem.u: u}) + Jp = replace(Jp, {problem.u_restrict: u}) else: Jp = None bcs = [] diff --git a/firedrake/variational_solver.py b/firedrake/variational_solver.py index 9f5fe14c01..ef420f5970 100644 --- a/firedrake/variational_solver.py +++ b/firedrake/variational_solver.py @@ -5,8 +5,12 @@ from firedrake import dmhooks, slate, solving, solving_utils, ufl_expr, utils from firedrake import function from firedrake.petsc import PETSc, OptionsManager, flatten_parameters -from firedrake.bcs import DirichletBC +from firedrake.function import Function +from firedrake.functionspace import RestrictedFunctionSpace +from firedrake.ufl_expr import TrialFunction, TestFunction +from firedrake.bcs import DirichletBC, EquationBC from firedrake.adjoint_utils import NonlinearVariationalProblemMixin, NonlinearVariationalSolverMixin +from ufl import replace __all__ = ["LinearVariationalProblem", "LinearVariationalSolver", @@ -44,7 +48,7 @@ class NonlinearVariationalProblem(NonlinearVariationalProblemMixin): def __init__(self, F, u, bcs=None, J=None, Jp=None, form_compiler_parameters=None, - is_linear=False): + is_linear=False, restrict=False): r""" :param F: the nonlinear form :param u: the :class:`.Function` to solve for @@ -57,22 +61,51 @@ def __init__(self, F, u, bcs=None, J=None, compiler (optional) :is_linear: internally used to check if all domain/bc forms are given either in 'A == b' style or in 'F == 0' style. + :param restrict: (optional) If `True`, use restricted function spaces, + that exclude Dirichlet boundary condition nodes, internally for + the test and trial spaces. """ - - self.bcs = solving._extract_bcs(bcs) - # Check form style consistency - self.is_linear = is_linear - is_form_consistent(self.is_linear, self.bcs) - self.Jp_eq_J = Jp is None - + V = u.function_space() + self.output_space = V self.u = u - self.F = F - self.Jp = Jp + if not isinstance(self.u, function.Function): raise TypeError("Provided solution is a '%s', not a Function" % type(self.u).__name__) + # Use the user-provided Jacobian. If none is provided, derive # the Jacobian from the residual. self.J = J or ufl_expr.derivative(F, u) + self.F = F + self.Jp = Jp + if bcs: + for bc in bcs: + if isinstance(bc, EquationBC): + restrict = False + self.restrict = restrict + + if restrict and bcs: + V_res = RestrictedFunctionSpace(V, boundary_set=set([bc.sub_domain for bc in bcs])) + bcs = [DirichletBC(V_res, bc.function_arg, bc.sub_domain) for bc in bcs] + self.u_restrict = Function(V_res).interpolate(u) + v_res, u_res = TestFunction(V_res), TrialFunction(V_res) + F_arg, = F.arguments() + replace_dict = {F_arg: v_res} + replace_dict[self.u] = self.u_restrict + self.F = replace(F, replace_dict) + v_arg, u_arg = self.J.arguments() + self.J = replace(self.J, {v_arg: v_res, u_arg: u_res}) + if self.Jp: + v_arg, u_arg = self.Jp.arguments() + self.Jp = replace(self.Jp, {v_arg: v_res, u_arg: u_res}) + self.restricted_space = V_res + else: + self.u_restrict = u + + self.bcs = solving._extract_bcs(bcs) + # Check form style consistency + self.is_linear = is_linear + is_form_consistent(self.is_linear, self.bcs) + self.Jp_eq_J = Jp is None # Argument checking check_pde_args(self.F, self.J, self.Jp) @@ -87,7 +120,7 @@ def dirichlet_bcs(self): @utils.cached_property def dm(self): - return self.u.function_space().dm + return self.u_restrict.function_space().dm class NonlinearVariationalSolver(OptionsManager, NonlinearVariationalSolverMixin): @@ -198,7 +231,7 @@ def update_diffusivity(current_solution): self._problem = problem self._ctx = ctx - self._work = problem.u.dof_dset.layout_vec.duplicate() + self._work = problem.u_restrict.dof_dset.layout_vec.duplicate() self.snes.setDM(problem.dm) ctx.set_function(self.snes) @@ -263,7 +296,7 @@ def solve(self, bounds=None): problem_dms.append(solution_dm) for dbc in problem.dirichlet_bcs(): - dbc.apply(problem.u) + dbc.apply(problem.u_restrict) if bounds is not None: lower, upper = bounds @@ -271,7 +304,7 @@ def solve(self, bounds=None): self.snes.setVariableBounds(lb, ub) work = self._work - with problem.u.dat.vec as u: + with problem.u_restrict.dat.vec as u: u.copy(work) with ExitStack() as stack: # Ensure options database has full set of options (so monitors @@ -283,10 +316,12 @@ def solve(self, bounds=None): self.snes.solve(None, work) work.copy(u) self._setup = True + if problem.restrict: + problem.u.interpolate(problem.u_restrict) solving_utils.check_snes_convergence(self.snes) # Grab the comm associated with the `_problem` and call PETSc's garbage cleanup routine - comm = self._problem.u.function_space().mesh()._comm + comm = self._problem.u_restrict.function_space().mesh()._comm PETSc.garbage_cleanup(comm=comm) @@ -296,7 +331,7 @@ class LinearVariationalProblem(NonlinearVariationalProblem): @PETSc.Log.EventDecorator() def __init__(self, a, L, u, bcs=None, aP=None, form_compiler_parameters=None, - constant_jacobian=False): + constant_jacobian=False, restrict=False): r""" :param a: the bilinear form :param L: the linear form @@ -311,6 +346,9 @@ def __init__(self, a, L, u, bcs=None, aP=None, Jacobian is constant (i.e. does not depend on varying fields). If your Jacobian does not change, set this flag to ``True``. + :param restrict: (optional) If `True`, use restricted function spaces, + that exclude Dirichlet boundary condition nodes, internally for + the test and trial spaces. """ # In the linear case, the Jacobian is the equation LHS. J = a @@ -326,7 +364,7 @@ def __init__(self, a, L, u, bcs=None, aP=None, super(LinearVariationalProblem, self).__init__(F, u, bcs, J, aP, form_compiler_parameters=form_compiler_parameters, - is_linear=True) + is_linear=True, restrict=restrict) self._constant_jacobian = constant_jacobian diff --git a/scripts/firedrake-install b/scripts/firedrake-install index d17b3df2ff..a25018e33a 100755 --- a/scripts/firedrake-install +++ b/scripts/firedrake-install @@ -1727,6 +1727,7 @@ if mode == "install": # numpy requires setuptools==59.2.0 and wheel==0.37.0 until it moves to meson build run_pip(["install", "-U", "setuptools==59.2.0"]) run_pip(["install", "-U", "hatch"]) +run_pip(["install", "-U", "editables"]) run_pip(["install", "-U", "pip"]) run_pip(["install", "-U", "wheel==0.37.0"]) diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py index f58339238f..00150e49f1 100644 --- a/tests/regression/test_eigensolver.py +++ b/tests/regression/test_eigensolver.py @@ -3,7 +3,7 @@ from firedrake import * -def evals(n, degree=1, mesh=None): +def evals(n, degree=1, mesh=None, restrict=False): '''We base this test on the 1D Poisson problem with Dirichlet boundary conditions, outlined in part 1 of Daniele Boffi's 'Finite element approximation of eigenvalue problems' Acta Numerica 2010''' @@ -18,7 +18,7 @@ def evals(n, degree=1, mesh=None): # Create eigenproblem with boundary conditions bc = DirichletBC(V, 0.0, "on_boundary") - eigenprob = LinearEigenproblem(a, bcs=bc, bc_shift=-6666.) + eigenprob = LinearEigenproblem(a, bcs=bc, bc_shift=-6666., restrict=restrict) # Create corresponding eigensolver, looking for n eigenvalues eigensolver = LinearEigensolver( @@ -38,19 +38,20 @@ def evals(n, degree=1, mesh=None): estimates[k] = eigensolver.eigenvalue(k).real - true_values[-1] = eigenprob.bc_shift - + if not restrict: + true_values[-1] = eigenprob.bc_shift # sort in case order of numerical and analytic values differs. return sorted(true_values), sorted(estimates) +@pytest.mark.parametrize("restrict", [True, False]) @pytest.mark.parametrize(('n', 'degree', 'tolerance'), [(5, 1, 1e-13), (10, 1, 1e-13), (20, 1, 1e-13), (30, 1, 1e-13)]) -def test_evals_1d(n, degree, tolerance): - true_values, estimates = evals(n, degree=degree) +def test_evals_1d(n, degree, tolerance, restrict): + true_values, estimates = evals(n, degree=degree, restrict=restrict) assert np.allclose(true_values, estimates, rtol=tolerance) diff --git a/tests/regression/test_mesh_from_plex.py b/tests/regression/test_mesh_from_plex.py index 2e43847211..f639187ebc 100644 --- a/tests/regression/test_mesh_from_plex.py +++ b/tests/regression/test_mesh_from_plex.py @@ -17,7 +17,7 @@ def get_plex_with_update_coordinates(mesh): gdim = mesh.geometric_dimension() entity_dofs = np.zeros(tdim + 1, dtype=np.int32) entity_dofs[0] = gdim - coord_section = mesh.create_section(entity_dofs) + coord_section, _ = mesh.create_section(entity_dofs) plex = mesh.topology_dm.clone() remove_pyop2_label(plex) diff --git a/tests/regression/test_restricted_function_space.py b/tests/regression/test_restricted_function_space.py new file mode 100644 index 0000000000..4a26e7834c --- /dev/null +++ b/tests/regression/test_restricted_function_space.py @@ -0,0 +1,207 @@ +import pytest +import numpy as np +from firedrake import * + + +def compare_function_space_assembly(function_space, restricted_function_space, + bcs, res_bcs=[]): + u = TrialFunction(function_space) + v = TestFunction(function_space) + original_form = inner(u, v) * dx + + u2 = TrialFunction(restricted_function_space) + v2 = TestFunction(restricted_function_space) + restricted_form = inner(u2, v2) * dx + + normal_fs_matrix = assemble(original_form, bcs=bcs) + restricted_fs_matrix = assemble(restricted_form, bcs=res_bcs) + + identity = np.identity(normal_fs_matrix.M.nrows) + delete_rows = [] + for i in range(normal_fs_matrix.M.nrows): + row = normal_fs_matrix.M.values[i, :] + comparison = row == identity[i, :] + if comparison.all(): + delete_rows.append(i) + + # Remove all rows/columns associated with boundaries (i.e. identity) + normal_fs_matrix_reduced = np.delete(normal_fs_matrix.M.values, delete_rows, + axis=0) + normal_fs_matrix_reduced = np.delete(normal_fs_matrix_reduced, delete_rows, + axis=1) + + assert (restricted_fs_matrix.M.nrows == np.shape(normal_fs_matrix_reduced)[0]) + assert (restricted_fs_matrix.M.ncols == np.shape(normal_fs_matrix_reduced)[1]) + assert (np.array_equal(normal_fs_matrix_reduced, restricted_fs_matrix.M.values)) + + +@pytest.mark.parametrize("j", [1, 2, 5]) +def test_restricted_function_space_1_1_square(j): + mesh = UnitSquareMesh(1, 1) + V = FunctionSpace(mesh, "CG", j) + bc = DirichletBC(V, 0, 2) + V_res = RestrictedFunctionSpace(V, name="Restricted", boundary_set=[2]) + res_bc = DirichletBC(V_res, 0, 2) + compare_function_space_assembly(V, V_res, [bc], [res_bc]) + + +@pytest.mark.parametrize("j", [1, 2, 5]) +def test_restricted_function_space_j_j_square(j): + mesh = UnitSquareMesh(j, j) + V = FunctionSpace(mesh, "CG", 1) + bc = DirichletBC(V, 0, 2) + V_res = RestrictedFunctionSpace(V, name="Restricted", boundary_set=[2]) + + compare_function_space_assembly(V, V_res, [bc]) + + +def test_poisson_homogeneous_bcs(): + mesh = UnitSquareMesh(1, 1) + V = FunctionSpace(mesh, "CG", 2) + f = Function(V) + x, y = SpatialCoordinate(mesh) + + u = TrialFunction(V) + v = TestFunction(V) + + f.interpolate(16 * pi**2 * (y-1)**2 * y**2 - 2 * (y-1)**2 - 8 * (y-1)*y + - 2*y**2)*sin(4*pi*x) + + original_form = inner(grad(u), grad(v)) * dx + + bc = DirichletBC(V, 0, 1) + V_res = RestrictedFunctionSpace(V, name="Restricted", boundary_set=[1]) + f2 = Function(V_res) + + f2.interpolate(16 * pi**2 * (y-1)**2 * y**2 - 2 * (y-1)**2 - 8 * (y-1)*y + - 2*y**2)*sin(4*pi*x) + + u2 = TrialFunction(V_res) + v2 = TestFunction(V_res) + restricted_form = inner(grad(u2), grad(v2)) * dx + + L = inner(f, v) * dx + L_res = inner(f2, v2) * dx + + u = Function(V) + bc2 = DirichletBC(V_res, Constant(0), 1) + u2 = Function(V_res) + + solve(original_form == L, u, bcs=[bc]) + solve(restricted_form == L_res, u2, bcs=[bc2]) + + assert errornorm(u, u2) < 1.e-12 + + +def test_poisson_inhomogeneous_bcs(): + mesh = UnitSquareMesh(1, 2) + V = FunctionSpace(mesh, "CG", 2) + V_res = RestrictedFunctionSpace(V, boundary_set=[1, 2]) + bc = DirichletBC(V_res, 0, 1) + bc2 = DirichletBC(V_res, 1, 2) + u2 = TrialFunction(V_res) + v2 = TestFunction(V_res) + restricted_form = inner(grad(u2), grad(v2)) * dx + u = Function(V_res) + rhs = inner(Constant(0), v2) * dx + + solve(restricted_form == rhs, u, bcs=[bc, bc2]) + + assert errornorm(SpatialCoordinate(mesh)[0], u) < 1.e-12 + + +@pytest.mark.parametrize("j", ["2", "sin(x) * y", "x**2"]) +def test_poisson_inhomogeneous_bcs_2(j): + mesh = UnitSquareMesh(1, 1) + V = FunctionSpace(mesh, "CG", 2) + + x, y = SpatialCoordinate(mesh) + + bc3 = DirichletBC(V, 0, 1) + bc4 = DirichletBC(V, Function(V).interpolate(eval(j)), 2) + u = TrialFunction(V) + v = TestFunction(V) + + V_res = RestrictedFunctionSpace(V, name="Restricted", boundary_set=[1, 2]) + + bc = DirichletBC(V_res, 0, 1) + bc2 = DirichletBC(V_res, Function(V_res).interpolate(eval(j)), 2) + + u2 = TrialFunction(V_res) + v2 = TestFunction(V_res) + + original_form = inner(grad(u), grad(v)) * dx + restricted_form = inner(grad(u2), grad(v2)) * dx + + u = Function(V) + u2 = Function(V_res) + + rhs = inner(Constant(0), v) * dx + rhs2 = inner(Constant(0), v2) * dx + + solve(original_form == rhs, u, bcs=[bc3, bc4]) + solve(restricted_form == rhs2, u2, bcs=[bc, bc2]) + + assert errornorm(u, u2) < 1.e-12 + + +@pytest.mark.parallel(nprocs=3) +def test_poisson_inhomogeneous_bcs_high_level_interface(): + mesh = UnitSquareMesh(8, 8) + V = FunctionSpace(mesh, "CG", 2) + bc1 = DirichletBC(V, 0., 1) + bc2 = DirichletBC(V, 1., 2) + u = TrialFunction(V) + v = TestFunction(V) + a = inner(grad(u), grad(v)) * dx + u = Function(V) + L = inner(Constant(0), v) * dx + solve(a == L, u, bcs=[bc1, bc2], restrict=True) + assert errornorm(SpatialCoordinate(mesh)[0], u) < 1.e-12 + + +@pytest.mark.parametrize("j", [1, 2, 5]) +def test_restricted_function_space_coord_change(j): + mesh = UnitSquareMesh(1, 2) + x, y = SpatialCoordinate(mesh) + V = VectorFunctionSpace(mesh, "CG", 2) + new_mesh = Mesh(Function(V).interpolate(as_vector([x, y]))) + new_V = FunctionSpace(new_mesh, "CG", j) + bc = DirichletBC(new_V, 0, 1) + new_V_restricted = RestrictedFunctionSpace(new_V, name="Restricted", boundary_set=[1]) + + compare_function_space_assembly(new_V, new_V_restricted, [bc]) + + +@pytest.mark.parametrize(["i", "j"], [(1, 0), (2, 0), (2, 1)]) +def test_restricted_mixed_spaces(i, j): + mesh = UnitSquareMesh(1, 1) + DG = FunctionSpace(mesh, "DG", j) + CG = VectorFunctionSpace(mesh, "CG", i) + CG_res = RestrictedFunctionSpace(CG, "Restricted", boundary_set=[4]) + W = CG * DG + W_res = CG_res * DG + bc = DirichletBC(W.sub(0), 0, 4) + bc_res = DirichletBC(W_res.sub(0), 0, 4) + + sigma2, u2 = TrialFunctions(W_res) + tau2, v2 = TestFunctions(W_res) + + sigma, u = TrialFunctions(W) + tau, v = TestFunctions(W) + + a = (inner(sigma, tau) + inner(u, div(tau)) + inner(div(sigma), v)) * dx + a2 = (inner(sigma2, tau2) + inner(u2, div(tau2)) + inner(div(sigma2), v2)) * dx + + w = Function(W) + w2 = Function(W_res) + + f = Constant(1) + L = - inner(f, v) * dx + L2 = - inner(f, v2) * dx + + solve(a == L, w, bcs=[bc]) + solve(a2 == L2, w2, bcs=[bc_res]) + + assert errornorm(w.subfunctions[0], w2.subfunctions[0]) < 1.e-12 + assert errornorm(w.subfunctions[1], w2.subfunctions[1]) < 1.e-12