diff --git a/firedrake/cython/mgimpl.pyx b/firedrake/cython/mgimpl.pyx index c6789e9ac1..811a6a6eed 100644 --- a/firedrake/cython/mgimpl.pyx +++ b/firedrake/cython/mgimpl.pyx @@ -352,6 +352,99 @@ def adaptive_parent_child_cell_maps(PETSc.DM coarse_dm, return np.asarray(coarse_to_fine), np.asarray(fine_to_coarse) +@cython.boundscheck(False) +@cython.wraparound(False) +def preserved_points(PETSc.DM coarse_dm, + PETSc.Section coarse_cell_numbering, + PETSc.DM fine_dm, + PETSc.Section fine_cell_numbering, + np.ndarray coarse_to_fine_cells): + """Pair unrefined fine points with their coarse originals. + + Adaptive refinement copies an untouched coarse cell into the fine mesh + without change. It therefore preserves the cone of every point in that + cell, and so preserves the whole plex closure, point for point. Such a + cell has exactly one child. The right-padding of ``coarse_to_fine_cells`` + with -1 identifies which cells these are. + + Parameters + ---------- + coarse_dm : PETSc.DM + The coarse mesh DMPlex. + coarse_cell_numbering : PETSc.Section + The cell numbering section of the coarse mesh. + fine_dm : PETSc.DM + The adaptively refined DMPlex. + fine_cell_numbering : PETSc.Section + The cell numbering section of the fine mesh. + coarse_to_fine_cells : numpy.ndarray + The Firedrake-numbered coarse-to-fine cell map. + + Returns + ------- + numpy.ndarray + An array over the chart of ``fine_dm``. For each fine point, it + holds the coarse point it was copied from, or -1 if refinement + changed that point. + + """ + cdef: + PetscInt ncoarse, nfine, max_children, c, i, off, child + PetscInt cStart, cEnd, pStart, pEnd, coarse_size, fine_size + PetscInt *coarse_closure = NULL + PetscInt *fine_closure = NULL + PetscInt[::1] coarse_point, fine_point, fine_to_coarse + PetscInt[:, ::1] coarse_to_fine + + coarse_to_fine = coarse_to_fine_cells + ncoarse = num_owned_cells(coarse_dm) + assert ncoarse == coarse_to_fine.shape[0] + max_children = coarse_to_fine.shape[1] + nfine = num_owned_cells(fine_dm) + + # Both cell maps are in Firedrake numbering, so invert each mesh's cell + # numbering section to get back to the plex points the closures live on. + coarse_point = np.full(ncoarse, -1, dtype=IntType) + cStart, cEnd = coarse_dm.getHeightStratum(0) + for c in range(cStart, cEnd): + CHKERR(PetscSectionGetOffset(coarse_cell_numbering.sec, c, &off)) + if 0 <= off < ncoarse: + coarse_point[off] = c + fine_point = np.full(nfine, -1, dtype=IntType) + cStart, cEnd = fine_dm.getHeightStratum(0) + for c in range(cStart, cEnd): + CHKERR(PetscSectionGetOffset(fine_cell_numbering.sec, c, &off)) + if 0 <= off < nfine: + fine_point[off] = c + + pStart, pEnd = fine_dm.getChart() + fine_to_coarse = np.full(pEnd - pStart, -1, dtype=IntType) + for c in range(ncoarse): + child = coarse_to_fine[c, 0] + if child < 0 or (max_children > 1 and coarse_to_fine[c, 1] >= 0): + continue + if coarse_point[c] < 0 or fine_point[child] < 0: + continue + CHKERR(DMPlexGetTransitiveClosure(coarse_dm.dm, coarse_point[c], PETSC_TRUE, + &coarse_size, &coarse_closure)) + CHKERR(DMPlexGetTransitiveClosure(fine_dm.dm, fine_point[child], PETSC_TRUE, + &fine_size, &fine_closure)) + # A one-child cell that refinement did change would have a closure + # of a different size. Skip it and let the transfer kernel handle it. + if coarse_size == fine_size: + for i in range(coarse_size): + # Each closure interleaves a point with its orientation. Copy + # a point only when its orientation matches in both meshes: + # only then do the two cells order their nodes the same way. + if coarse_closure[2*i + 1] == fine_closure[2*i + 1]: + fine_to_coarse[fine_closure[2*i] - pStart] = coarse_closure[2*i] + CHKERR(DMPlexRestoreTransitiveClosure(coarse_dm.dm, coarse_point[c], PETSC_TRUE, + &coarse_size, &coarse_closure)) + CHKERR(DMPlexRestoreTransitiveClosure(fine_dm.dm, fine_point[child], PETSC_TRUE, + &fine_size, &fine_closure)) + return np.asarray(fine_to_coarse) + + # Exposition: # # These next functions compute maps from coarse mesh cells to fine diff --git a/firedrake/mg/interface.py b/firedrake/mg/interface.py index 186e0acddd..6c42d3c434 100644 --- a/firedrake/mg/interface.py +++ b/firedrake/mg/interface.py @@ -102,7 +102,11 @@ def prolong(coarse, fine): for d in [coarse, coarse_coords]: d.dat.global_to_local_begin(op2.READ) d.dat.global_to_local_end(op2.READ) - op2.par_loop(kernel, fine.node_set, *kernel_args) + # Adaptive refinement leaves most of the mesh unchanged. Copy the + # value at those nodes instead of evaluating it there. + node_subset = utils.transfer_node_subset(Vc, Vf) + op2.par_loop(kernel, node_subset, *kernel_args) + utils.prolong_preserved_nodes(coarse, fine) if needs_quadrature: # Transfer to the actual target space @@ -184,7 +188,11 @@ def restrict(fine_dual, coarse_dual): for d in [coarse_coords]: d.dat.global_to_local_begin(op2.READ) d.dat.global_to_local_end(op2.READ) - op2.par_loop(kernel, fine_dual.node_set, *kernel_args) + # Restriction is the transpose of prolongation. It skips the same + # fine nodes and adds their value to the matching coarse node. + node_subset = utils.transfer_node_subset(Vc, Vf) + op2.par_loop(kernel, node_subset, *kernel_args) + utils.restrict_preserved_nodes(fine_dual, coarse_dual) fine_dual = coarse_dual return coarse_dual diff --git a/firedrake/mg/utils.py b/firedrake/mg/utils.py index 3bbd52c77d..042bf5edaf 100644 --- a/firedrake/mg/utils.py +++ b/firedrake/mg/utils.py @@ -2,12 +2,15 @@ import numpy from fractions import Fraction +from mpi4py import MPI from pyop2 import op2 +from firedrake.petsc import PETSc from firedrake.utils import IntType from firedrake.functionspacedata import entity_dofs_key import finat.ufl import firedrake from firedrake.cython import mgimpl as impl +from firedrake.halo import _get_mtype def fine_node_to_coarse_node_map(Vf, Vc): @@ -210,6 +213,229 @@ def coarse_cell_child_count( return cache.setdefault(key, op2.Dat(dset, counts, dtype=IntType)) +def _preserved_point_sf(coarse_mesh, fine_mesh, coarse_to_fine): + """Create the SF that pairs unrefined points with their coarse originals. + + Adaptive refinement leaves some cells untouched. This SF maps each + unrefined point in ``fine_mesh`` back to the coarse point it came from. + + Parameters + ---------- + coarse_mesh : firedrake.mesh.AbstractMeshTopology + The mesh before refinement. + fine_mesh : firedrake.mesh.AbstractMeshTopology + The mesh after refinement. + coarse_to_fine : numpy.ndarray + The coarse-to-fine cell map that relates the two meshes. + + Returns + ------- + PETSc.SF + An SF with roots on the points of ``coarse_mesh`` and leaves on the + unrefined points of ``fine_mesh``. It has no leaves where refinement + changed every cell, as a uniform refinement does. + + """ + coarse_plex = coarse_mesh.topology_dm + fine_plex = fine_mesh.topology_dm + fine_to_coarse_points = impl.preserved_points( + coarse_plex, coarse_mesh._cell_numbering, + fine_plex, fine_mesh._cell_numbering, + coarse_to_fine, + ) + leaves = numpy.nonzero(fine_to_coarse_points >= 0)[0].astype(IntType) + # Refinement acts on each rank's own plex. A fine point and the coarse + # point it was copied from always live on the same rank. + remote = numpy.empty((len(leaves), 2), dtype=IntType) + remote[:, 0] = coarse_plex.comm.rank + remote[:, 1] = fine_to_coarse_points[leaves] + pStart, pEnd = coarse_plex.getChart() + point_sf = PETSc.SF().create(comm=coarse_plex.comm) + point_sf.setGraph(pEnd - pStart, leaves, remote) + return point_sf + + +def preserved_node_sf( + Vc: firedrake.functionspaceimpl.WithGeometry, + Vf: firedrake.functionspaceimpl.WithGeometry, +) -> PETSc.SF | None: + """Find the nodes that adaptive refinement leaves unchanged. + + An unrefined cell has the same nodes in both spaces. The transfer + operators can then copy values between them instead of evaluating them. + This is cheaper, and exact. + + Parameters + ---------- + Vc : firedrake.functionspaceimpl.WithGeometry + The coarse function space. + Vf : firedrake.functionspaceimpl.WithGeometry + The fine function space, on the next level of the same hierarchy. + + Returns + ------- + PETSc.SF or None + An SF with roots on the nodes of ``Vc`` and leaves on the owned nodes + of ``Vf`` that match. `None` says that no node matches anywhere, so + that a caller can skip the copy and evaluate every node instead. + + Notes + ----- + Four things stop any node from matching, and every rank reaches the same + verdict on each of them: + + * the two spaces lay their nodes out differently, because their elements + or their boundary sets differ; + * one of the meshes is extruded. Its DMPlex holds the base mesh alone, so + a point there stands for a whole vertical column of nodes and a + `PETSc.Section` cannot address one node within that column; + * the two meshes are not consecutive levels of one hierarchy; + * the refinement is uniform, and so rebuilt every cell. + + """ + if Vc.ufl_element() != Vf.ufl_element() or Vc.boundary_set != Vf.boundary_set: + return None + if Vc.extruded or Vf.extruded: + return None + hierarchy, levelc = get_level(Vc.mesh()) + _, levelf = get_level(Vf.mesh()) + if hierarchy is None or levelc + Fraction(1, hierarchy.refinements_per_level) != levelf: + return None + cache = Vf.mesh().topology._shared_data_cache["hierarchy_preserved_node_sf"] + key = _cache_key(Vc, Vf) + try: + return cache[key] + except KeyError: + coarse_to_fine = hierarchy.coarse_to_fine_cells[levelc] + point_sf = _preserved_point_sf(Vc.mesh().topology, Vf.mesh().topology, + coarse_to_fine) + root_section = Vc.dm.getSection() + leaf_section = Vf.dm.getSection() + # `distributeSection` builds its own section over the range of points + # that the SF touches. Only the broadcast root offsets are needed + # here. Pad them back out to the full chart that `createSectionSF` + # expects. + remote_offsets, distributed_section = point_sf.distributeSection(root_section) + pStart, pEnd = leaf_section.getChart() + lpStart, lpEnd = distributed_section.getChart() + offsets = numpy.zeros(pEnd - pStart, dtype=IntType) + offsets[lpStart - pStart:lpEnd - pStart] = remote_offsets + section_sf = point_sf.createSectionSF(root_section, offsets, leaf_section) + # The transfer kernels compute only the owned fine nodes and leave + # the halo to a later exchange. Keep only the owned leaves here too: + # a ghost fine node reduced onto its coarse node would count twice. + nroots, ilocal, iremote = section_sf.getGraph() + owned = ilocal < Vf.node_set.size + # Every rank must agree on whether to copy or to evaluate, or they + # generate different code. A rank with no owned leaf of its own still + # takes part where another rank has one. + if not Vf.mesh().comm.allreduce(bool(owned.any()), op=MPI.LOR): + return cache.setdefault(key, None) + trimmed = PETSc.SF().create(comm=section_sf.comm) + trimmed.setGraph(nroots, ilocal[owned], iremote[owned]) + return cache.setdefault(key, trimmed) + + +def transfer_node_subset( + Vc: firedrake.functionspaceimpl.WithGeometry, + Vf: firedrake.functionspaceimpl.WithGeometry, +) -> op2.Set: + """Find the fine nodes that the transfer kernels must evaluate. + + These are the nodes of ``Vf`` that :func:`preserved_node_sf` does not + already account for. Prolongation and restriction can copy the rest. + + Parameters + ---------- + Vc : firedrake.functionspaceimpl.WithGeometry + The coarse function space. + Vf : firedrake.functionspaceimpl.WithGeometry + The fine function space, on the next level of the same hierarchy. + + Returns + ------- + pyop2.types.set.Set or pyop2.types.set.Subset + A subset of the nodes of ``Vf``. Where nothing is preserved this is + ``Vf.node_set`` itself, which spares the kernel a level of + indirection that would index every node anyway. + + """ + section_sf = preserved_node_sf(Vc, Vf) + if section_sf is None: + return Vf.node_set + cache = Vf.mesh().topology._shared_data_cache["hierarchy_transfer_node_subset"] + key = _cache_key(Vc, Vf) + try: + return cache[key] + except KeyError: + _, preserved, _ = section_sf.getGraph() + nodes = numpy.setdiff1d(numpy.arange(Vf.node_set.size, dtype=IntType), + preserved) + return cache.setdefault(key, op2.Subset(Vf.node_set, nodes)) + + +def prolong_preserved_nodes(coarse, fine): + """Copy coarse values onto the fine nodes that adaptive refinement preserved. + + Parameters + ---------- + coarse : firedrake.function.Function + The function on the coarse mesh. + fine : firedrake.function.Function + The function on the refined mesh. The transfer kernel has already + computed its other nodes. + + """ + section_sf = preserved_node_sf(coarse.function_space(), fine.function_space()) + if section_sf is None: + return + mtype, _ = _get_mtype(fine.dat) + # The source coarse node can be a ghost node. Only owned fine nodes are + # written here, the same as the transfer kernel writes. + source = coarse.dat.data_ro_with_halos + target = fine.dat.data_wo + section_sf.bcastBegin(mtype, source, target, MPI.REPLACE) + section_sf.bcastEnd(mtype, source, target, MPI.REPLACE) + + +def restrict_preserved_nodes(fine_dual, coarse_dual): + """Add the contribution of preserved nodes to the coarse dual. + + Prolongation copies a preserved node's value without change. Restriction + is its transpose, so it adds the fine value to the coarse node unchanged. + + Parameters + ---------- + fine_dual : firedrake.cofunction.Cofunction + The cofunction on the refined mesh. + coarse_dual : firedrake.cofunction.Cofunction + The cofunction on the coarse mesh. It already holds the contribution + that the transfer kernel accumulated from the other fine nodes. + + Notes + ----- + This adds rather than copies, at every stage. A coarse basis function + does not vanish on the cells around a preserved node that refinement did + split, so the same coarse node also collects a contribution from the + kernel. Several preserved fine nodes can likewise reduce onto one coarse + node, and the coarse node can be a ghost on the rank that owns them. + + """ + coarse_V = coarse_dual.function_space() + section_sf = preserved_node_sf(coarse_V, fine_dual.function_space()) + if section_sf is None: + return + buffer = firedrake.Function(coarse_V) + mtype, _ = _get_mtype(buffer.dat) + source = fine_dual.dat.data_ro + target = buffer.dat.data_wo_with_halos + section_sf.reduceBegin(mtype, source, target, MPI.SUM) + section_sf.reduceEnd(mtype, source, target, MPI.SUM) + buffer.dat.local_to_global_begin(op2.INC) + buffer.dat.local_to_global_end(op2.INC) + coarse_dual.dat.data[...] += buffer.dat.data_ro + + def physical_node_locations(V): element = V.ufl_element() if V.value_shape: diff --git a/tests/firedrake/multigrid/test_adaptive_multigrid.py b/tests/firedrake/multigrid/test_adaptive_multigrid.py index ac9e557d0a..bfdc63abe1 100644 --- a/tests/firedrake/multigrid/test_adaptive_multigrid.py +++ b/tests/firedrake/multigrid/test_adaptive_multigrid.py @@ -2,6 +2,7 @@ import numpy as np from mpi4py import MPI from firedrake import * +from firedrake.mg.utils import preserved_node_sf from firedrake.utils import complex_mode @@ -34,8 +35,9 @@ def _linear_expr(mesh): def coarse_mesh(request): dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)} mesher = request.param - # Big enough that refining part of it leaves untouched cells behind, and - # that a coarse cell's child count varies widely across the mesh. + # Big enough that refining part of it leaves untouched cells behind. + # The transfers copy those cells' nodes instead of evaluating them. + # It also gives a coarse cell's child count a wide range. if mesher == "firedrake-square": return UnitSquareMesh(4, 4, distribution_parameters=dparams) elif mesher == "firedrake-cube": @@ -116,8 +118,10 @@ def test_refine_marked_elements_is_local(): @pytest.mark.parallel([1, 2]) def test_refine_marked_elements_repeats(coarse_mesh): - """A marker value of n refines the marked cells n times, and the cell maps - reach all the way from the original mesh to the n-times-refined one.""" + """A marker value of n refines the marked cells n times. + + The cell maps reach all the way from the original mesh to the + n-times-refined one.""" mesh = coarse_mesh ncells = {} max_children = {} @@ -148,8 +152,9 @@ def test_refine_marked_elements_repeats(coarse_mesh): def test_add_mesh_rejects_unrelated_mesh(): """Cell maps are only meaningful relative to the mesh they were built - against, so a mesh refined from anything but the finest level is refused - rather than silently recorded with somebody else's maps.""" + against. Refuse a mesh refined from anything but the finest level, + instead of silently recording it with maps that belong to another + mesh.""" mh = MeshHierarchy(UnitSquareMesh(2, 2)) other = UnitSquareMesh(4, 4) @@ -215,10 +220,11 @@ def test_CG1_native_transfers_use_adaptive_cell_maps(coarse_mesh): def _assert_adapt_after_uniform_refinement(mh): - """Adaptively refine the finest level of the uniformly-refined hierarchy - ``mh`` by marking a single cell, and check that the cell maps of the level - this adds are sane. Shared by the ``test_adapt_after_uniform_*refinement`` - tests, which only differ in how ``mh`` itself was built. + """Adaptively refine the finest level of the hierarchy ``mh``. + + Mark a single cell and check that the cell maps of the new level are + sane. The ``test_adapt_after_uniform_*refinement`` tests share this + helper; they differ only in how they build ``mh``. """ mesh = mh[-1] level = len(mh) @@ -264,9 +270,10 @@ def test_adapt_after_uniform_netgen_refinement(): @pytest.mark.parallel([1, 2]) @pytest.mark.parametrize("degree", [1, 2]) def test_adapt_preserves_mesh_metadata(degree): - """Adaptive refinement carries the Netgen geometry and flags, and the mesh - construction parameters, over to the refined mesh, so that the refined - mesh can itself be refined again.""" + """Adaptive refinement carries mesh metadata to the refined mesh. + + It copies the Netgen geometry, flags, and construction parameters, so + the refined mesh can itself be refined again.""" from netgen.geom2d import CSG2d, Circle geo = CSG2d() geo.Add(Circle(center=(0, 0), radius=1.0, bc="circle")) @@ -304,8 +311,9 @@ def test_adapt_after_uniform_refinement(coarse_mesh, refine): @pytest.mark.parametrize("refine", [1, 2]) def test_adapt_before_uniform_refinement(coarse_mesh, refine): """An adaptively refined mesh can be uniformly refined into a hierarchy. - Its plex numbers cells by refinement case, so its owned cells are - interleaved with its halo cells, which the cell maps must not assume away. + + Its plex numbers cells by refinement case. This interleaves owned cells + with halo cells, and the cell maps must not assume otherwise. """ netgen_flags = {} if hasattr(coarse_mesh, "netgen_mesh") else None @@ -332,36 +340,27 @@ def test_adapt_before_uniform_refinement(coarse_mesh, refine): assert (fine_to_coarse[coarse_to_fine, 0] == parents).all() -@pytest.mark.parallel([1, 2, 4]) -@pytest.mark.parametrize("operator", ["prolong", "inject"]) -def test_DG0(mh, operator): - """Prolongation & Injection test for DG0""" - V_coarse = FunctionSpace(mh[0], "DG", 0) - V_fine = FunctionSpace(mh[-1], "DG", 0) - u_coarse = Function(V_coarse) - u_fine = Function(V_fine) - xc, *_ = SpatialCoordinate(V_coarse.mesh()) - stepc = conditional(ge(xc, 0), 1, 0) - xf, *_ = SpatialCoordinate(V_fine.mesh()) - stepf = conditional(ge(xf, 0), 1, 0) - - if operator == "prolong": - u_coarse.interpolate(stepc) - assert errornorm(stepc, u_coarse) <= 1e-12 - - prolong(u_coarse, u_fine) - assert errornorm(stepf, u_fine) <= 1e-12 - if operator == "inject": - u_fine.interpolate(stepf) - assert errornorm(stepf, u_fine) <= 1e-12 - - if complex_mode: - with pytest.raises(NotImplementedError): - inject(u_fine, u_coarse) - return - else: - inject(u_fine, u_coarse) - assert errornorm(stepc, u_coarse) <= 1e-12 +def _representable_expr(mesh, degree): + """An expression that a space of the given degree holds exactly on any mesh.""" + x = SpatialCoordinate(mesh) + if degree == 0: + # A DG0 space holds an expression exactly on every mesh of the + # hierarchy only if it is constant on each coarsest cell. + return conditional(ge(x[0], 0), 1, 0) + return sum(xi ** degree for xi in x) + + +def _copied_nodes(mh, V): + """Count the nodes of ``V`` that the transfers copy rather than evaluate.""" + copied = 0 + for level in range(len(mh) - 1): + V_coarse = V.reconstruct(mesh=mh[level]) + V_fine = V.reconstruct(mesh=mh[level + 1]) + section_sf = preserved_node_sf(V_coarse, V_fine) + if section_sf is not None: + _, preserved, _ = section_sf.getGraph() + copied += len(preserved) + return mh[0].comm.allreduce(copied, MPI.SUM) def _coarse_cell_integrals(mh, level, u_coarse, u_fine): @@ -432,74 +431,49 @@ def test_dg_injection_conserves_mass(mh, family, degree): @pytest.mark.parallel([1, 2, 4]) -@pytest.mark.parametrize("operator", ["prolong", "inject"]) -def test_CG1(mh, operator): - """Prolongation & Injection test for CG1""" - V_coarse = FunctionSpace(mh[0], "CG", 1) - V_fine = FunctionSpace(mh[-1], "CG", 1) - u_coarse = Function(V_coarse) - u_fine = Function(V_fine) - xc, *_ = SpatialCoordinate(V_coarse.mesh()) - xf, *_ = SpatialCoordinate(V_fine.mesh()) - - if operator == "prolong": - u_coarse.interpolate(xc) - assert errornorm(xc, u_coarse) <= 1e-12 +@pytest.mark.parametrize("family, degree", [("DG", 0), ("CG", 1), ("CG", 2), ("CG", 3)]) +def test_transfers(mh, family, degree): + """Prolongation, injection and restriction on an adaptive hierarchy. - prolong(u_coarse, u_fine) - assert errornorm(xf, u_fine) <= 1e-12 - if operator == "inject": - u_fine.interpolate(xf) - assert errornorm(xf, u_fine) <= 1e-12 - - inject(u_fine, u_coarse) - assert errornorm(xc, u_coarse) <= 1e-12 + Degree 3 puts more than one node on an edge. This catches a copied + entity whose nodes come out in the wrong order. + """ + V_coarse = FunctionSpace(mh[0], family, degree) + V_fine = FunctionSpace(mh[-1], family, degree) + expr_coarse = _representable_expr(mh[0], degree) + expr_fine = _representable_expr(mh[-1], degree) + # The cells that refinement leaves alone are transferred by copying + # their nodes. A hierarchy that refines everything says nothing about + # them. + assert _copied_nodes(mh, V_coarse) > 0 -@pytest.mark.parallel([1, 2, 4]) -def test_restrict_CG1(mh): - """Test restriction with CG1""" - V_coarse = FunctionSpace(mh[0], "CG", 1) - V_fine = FunctionSpace(mh[-1], "CG", 1) - u_coarse = Function(V_coarse) + u_coarse = Function(V_coarse).interpolate(expr_coarse) u_fine = Function(V_fine) - xc, *_ = SpatialCoordinate(V_coarse.mesh()) - - u_coarse.interpolate(xc) prolong(u_coarse, u_fine) + assert errornorm(expr_fine, u_fine) <= 1e-12 - rf = assemble(conj(TestFunction(V_fine)) * dx) - rc = Cofunction(V_coarse.dual()) - restrict(rf, rc) - + # Restriction is the transpose of prolongation. This ties the two + # together: restriction must not also accumulate, through the kernel, + # the same nodes that prolongation copies. + r_fine = assemble(conj(TestFunction(V_fine)) * dx) + r_coarse = Cofunction(V_coarse.dual()) + restrict(r_fine, r_coarse) assert np.allclose( - assemble(action(rc, u_coarse)), - assemble(action(rf, u_fine)), + assemble(action(r_coarse, u_coarse)), + assemble(action(r_fine, u_fine)), rtol=1e-12 ) - -@pytest.mark.parallel([1, 2, 4]) -def test_restrict_DG0(mh): - """Test restriction with DG0""" - V_coarse = FunctionSpace(mh[0], "DG", 0) - V_fine = FunctionSpace(mh[-1], "DG", 0) - u_coarse = Function(V_coarse) - u_fine = Function(V_fine) - xc, *_ = SpatialCoordinate(V_coarse.mesh()) - - u_coarse.interpolate(xc) - prolong(u_coarse, u_fine) - - rf = assemble(conj(TestFunction(V_fine)) * dx) - rc = Cofunction(V_coarse.dual()) - restrict(rf, rc) - - assert np.allclose( - assemble(action(rc, u_coarse)), - assemble(action(rf, u_fine)), - rtol=1e-12 - ) + # Injection + u_fine = Function(V_fine).interpolate(expr_fine) + u_injected = Function(V_coarse) + if family in {"DG", "DQ"} and complex_mode: + with pytest.raises(NotImplementedError): + inject(u_fine, u_injected) + else: + inject(u_fine, u_injected) + assert errornorm(expr_coarse, u_injected) <= 1e-12 @pytest.mark.parallel([1, 2])