From 081789dccc75ca780b0df952b59329bcde6319d6 Mon Sep 17 00:00:00 2001 From: Daiane Iglesia Dolci <63597005+Ig-dolci@users.noreply.github.com> Date: Fri, 29 Nov 2024 10:07:28 +0000 Subject: [PATCH 1/5] Enable vtk write cofunction (#3892) --- firedrake/output/vtk_output.py | 17 ++++++++++++----- tests/firedrake/output/test_pvd_output.py | 22 ++++++++++++++++------ 2 files changed, 28 insertions(+), 11 deletions(-) diff --git a/firedrake/output/vtk_output.py b/firedrake/output/vtk_output.py index 2dc7ab0235..23072c5019 100644 --- a/firedrake/output/vtk_output.py +++ b/firedrake/output/vtk_output.py @@ -452,7 +452,7 @@ def __init__(self, filename, project_output=False, comm=None, mode="w", @no_annotations def _prepare_output(self, function, max_elem): from firedrake import FunctionSpace, VectorFunctionSpace, \ - TensorFunctionSpace, Function + TensorFunctionSpace, Function, Cofunction name = function.name() # Need to project/interpolate? @@ -477,8 +477,15 @@ def _prepare_output(self, function, max_elem): shape=shape) else: raise ValueError("Unsupported shape %s" % (shape, )) - output = Function(V) + if isinstance(function, Function): + output = Function(V) + else: + assert isinstance(function, Cofunction) + output = Function(V.dual()) + if self.project: + if isinstance(function, Cofunction): + raise ValueError("Can not project Cofunctions") output.project(function) else: output.interpolate(function) @@ -486,7 +493,7 @@ def _prepare_output(self, function, max_elem): return OFunction(array=get_array(output), name=name, function=output) def _write_vtu(self, *functions): - from firedrake.function import Function + from firedrake import Function, Cofunction # Check if the user has requested to write out a plain mesh if len(functions) == 1 and isinstance(functions[0], ufl.Mesh): @@ -496,8 +503,8 @@ def _write_vtu(self, *functions): functions = [Function(V)] for f in functions: - if not isinstance(f, Function): - raise ValueError("Can only output Functions or a single mesh, not %r" % type(f)) + if not isinstance(f, (Function, Cofunction)): + raise ValueError(f"Can only output Functions, Cofunctions or a single mesh, not {type(f).__name__}") meshes = tuple(extract_unique_domain(f) for f in functions) if not all(m == meshes[0] for m in meshes): raise ValueError("All functions must be on same mesh") diff --git a/tests/firedrake/output/test_pvd_output.py b/tests/firedrake/output/test_pvd_output.py index 6f17bcbe48..69c7243b7f 100644 --- a/tests/firedrake/output/test_pvd_output.py +++ b/tests/firedrake/output/test_pvd_output.py @@ -79,11 +79,16 @@ def test_bad_file_name(tmpdir): VTKFile(str(tmpdir.join("foo.vtu"))) -def test_different_functions(mesh, pvd): +@pytest.mark.parametrize("space", + ["primal", "dual"]) +def test_different_functions(mesh, pvd, space): V = FunctionSpace(mesh, "DG", 0) - - f = Function(V, name="foo") - g = Function(V, name="bar") + if space == "primal": + f = Function(V, name="foo") + g = Function(V, name="bar") + else: + f = Cofunction(V.dual(), name="foo") + g = Cofunction(V.dual(), name="bar") pvd.write(f) @@ -136,9 +141,14 @@ def test_not_function(mesh, pvd): pvd.write(grad(f)) -def test_append(mesh, tmpdir): +@pytest.mark.parametrize("space", + ["primal", "dual"]) +def test_append(mesh, tmpdir, space): V = FunctionSpace(mesh, "DG", 0) - g = Function(V) + if space == "primal": + g = Function(V) + else: + g = Cofunction(V.dual()) outfile = VTKFile(str(tmpdir.join("restart_test.pvd"))) outfile.write(g) From 2be9fb27f067a2b68bc86ff06b754bfd1c21a699 Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Fri, 29 Nov 2024 14:20:27 +0000 Subject: [PATCH 2/5] Function has fs.block_size components not fs.value_size (#3899) --- firedrake/function.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/function.py b/firedrake/function.py index 0eed4a5138..100c02be87 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -327,7 +327,7 @@ def split(self): @utils.cached_property def _components(self): - if self.function_space().value_size == 1: + if self.function_space().block_size == 1: return (self, ) else: return tuple(type(self)(self.function_space().sub(i), self.topological.sub(i)) From 4c9335439f5b89a51124e524d5d57381b9928e17 Mon Sep 17 00:00:00 2001 From: India Marsden <37078108+indiamai@users.noreply.github.com> Date: Mon, 2 Dec 2024 11:09:39 +0000 Subject: [PATCH 3/5] Modify dictionary to function on cellnames (#3900) --- firedrake/mesh.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index fd37e7b575..906547ae72 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -66,10 +66,10 @@ } -_supported_embedded_cell_types_and_gdims = [(ufl.Cell('interval'), 2), - (ufl.Cell('triangle'), 3), - (ufl.Cell("quadrilateral"), 3), - (ufl.TensorProductCell(ufl.Cell('interval'), ufl.Cell('interval')), 3)] +_supported_embedded_cell_types_and_gdims = [('interval', 2), + ('triangle', 3), + ("quadrilateral", 3), + ("interval * interval", 3)] unmarked = -1 @@ -2765,7 +2765,7 @@ def init_cell_orientations(self, expr): import firedrake.function as function import firedrake.functionspace as functionspace - if (self.ufl_cell(), self.geometric_dimension()) not in _supported_embedded_cell_types_and_gdims: + if (self.ufl_cell().cellname(), self.geometric_dimension()) not in _supported_embedded_cell_types_and_gdims: raise NotImplementedError('Only implemented for intervals embedded in 2d and triangles and quadrilaterals embedded in 3d') if hasattr(self, '_cell_orientations'): From 530243efb0a1dce6a4231198babafbdd24f621f1 Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Tue, 3 Dec 2024 17:05:12 +0000 Subject: [PATCH 4/5] `u.sub(0)` should still be in a `ComponentFunctionSpace` if the vector is 1 dimensional (#3902) Co-authored-by: Connor Ward --- firedrake/function.py | 36 ++++++++++++++++++++++------------ firedrake/functionspaceimpl.py | 19 +++++++++++------- 2 files changed, 35 insertions(+), 20 deletions(-) diff --git a/firedrake/function.py b/firedrake/function.py index 100c02be87..585a1511d1 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -124,12 +124,16 @@ def split(self): @utils.cached_property def _components(self): - if self.dof_dset.cdim == 1: - return (self, ) + if self.function_space().rank == 0: + return tuple((self, )) else: - return tuple(CoordinatelessFunction(self.function_space().sub(i), val=op2.DatView(self.dat, j), - name="view[%d](%s)" % (i, self.name())) - for i, j in enumerate(np.ndindex(self.dof_dset.dim))) + if self.dof_dset.cdim == 1: + return (CoordinatelessFunction(self.function_space().sub(0), val=self.dat, + name=f"view[0]({self.name()})"),) + else: + return tuple(CoordinatelessFunction(self.function_space().sub(i), val=op2.DatView(self.dat, j), + name=f"view[{i}]({self.name()})") + for i, j in enumerate(np.ndindex(self.dof_dset.dim))) @PETSc.Log.EventDecorator() def sub(self, i): @@ -143,9 +147,12 @@ def sub(self, i): rank-n :class:`~.FunctionSpace`, this returns a proxy object indexing the ith component of the space, suitable for use in boundary condition application.""" - if len(self.function_space()) == 1: - return self._components[i] - return self.subfunctions[i] + mixed = len(self.function_space()) != 1 + data = self.subfunctions if mixed else self._components + bound = len(data) + if i < 0 or i >= bound: + raise IndexError(f"Invalid component {i}, not in [0, {bound})") + return data[i] @property def cell_set(self): @@ -327,8 +334,8 @@ def split(self): @utils.cached_property def _components(self): - if self.function_space().block_size == 1: - return (self, ) + if self.function_space().rank == 0: + return tuple((self, )) else: return tuple(type(self)(self.function_space().sub(i), self.topological.sub(i)) for i in range(self.function_space().block_size)) @@ -345,9 +352,12 @@ def sub(self, i): :func:`~.VectorFunctionSpace` or :func:`~.TensorFunctionSpace` this returns a proxy object indexing the ith component of the space, suitable for use in boundary condition application.""" - if len(self.function_space()) == 1: - return self._components[i] - return self.subfunctions[i] + mixed = len(self.function_space()) != 1 + data = self.subfunctions if mixed else self._components + bound = len(data) + if i < 0 or i >= bound: + raise IndexError(f"Invalid component {i}, not in [0, {bound})") + return data[i] @PETSc.Log.EventDecorator() @FunctionMixin._ad_annotate_project diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py index 53f69e92ce..b8748ed535 100644 --- a/firedrake/functionspaceimpl.py +++ b/firedrake/functionspaceimpl.py @@ -182,10 +182,12 @@ def _components(self): @PETSc.Log.EventDecorator() def sub(self, i): - bound = len(self._components) + mixed = len(self) != 1 + data = self.subfunctions if mixed else self._components + bound = len(data) if i < 0 or i >= bound: - raise IndexError("Invalid component %d, not in [0, %d)" % (i, bound)) - return self._components[i] + raise IndexError(f"Invalid component {i}, not in [0, {bound})") + return data[i] @utils.cached_property def dm(self): @@ -654,13 +656,16 @@ def __getitem__(self, i): @utils.cached_property def _components(self): - return tuple(ComponentFunctionSpace(self, i) for i in range(self.block_size)) + if self.rank == 0: + return self.subfunctions + else: + return tuple(ComponentFunctionSpace(self, i) for i in range(self.block_size)) def sub(self, i): r"""Return a view into the ith component.""" - if self.rank == 0: - assert i == 0 - return self + bound = len(self._components) + if i < 0 or i >= bound: + raise IndexError(f"Invalid component {i}, not in [0, {bound})") return self._components[i] def __mul__(self, other): From 452b34243d011a3f8d1a5fb6ab89c096f760c555 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Wed, 4 Dec 2024 12:37:12 +0000 Subject: [PATCH 5/5] PETSc fixes (#3884) * Ravel swarm.getField values to work with updated petsc4py API * Update some install dependencies to work with latest PETSc/SLEPc. --------- Co-authored-by: Scott MacLachlan --- docker/Dockerfile.env | 2 +- docs/source/download.rst | 1 + firedrake/mesh.py | 32 +++++++++---------- scripts/firedrake-install | 4 ++- tests/firedrake/vertexonly/test_swarm.py | 18 +++++------ .../vertexonly/test_vertex_only_fs.py | 4 +-- .../test_vertex_only_mesh_generation.py | 2 +- 7 files changed, 33 insertions(+), 30 deletions(-) diff --git a/docker/Dockerfile.env b/docker/Dockerfile.env index a0ad86b989..ac15aa7e2f 100644 --- a/docker/Dockerfile.env +++ b/docker/Dockerfile.env @@ -13,7 +13,7 @@ RUN apt-get update \ cmake gfortran git libopenblas-serial-dev \ libtool python3-dev python3-pip python3-tk python3-venv \ python3-requests zlib1g-dev libboost-dev sudo gmsh \ - bison flex ninja-build \ + bison flex ninja-build pkg-config \ libocct-ocaf-dev libocct-data-exchange-dev \ swig graphviz \ libcurl4-openssl-dev libxml2-dev \ diff --git a/docs/source/download.rst b/docs/source/download.rst index 30c37e78b0..3865e8b175 100644 --- a/docs/source/download.rst +++ b/docs/source/download.rst @@ -135,6 +135,7 @@ they have the system dependencies: * zlib * flex, bison * Ninja +* pkg-config Firedrake has been successfully installed on Windows 10 using the Windows Subsystem for Linux. There are more detailed instructions for diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 906547ae72..7c99b47972 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -1966,7 +1966,7 @@ def _renumber_entities(self, reorder): if reorder: swarm = self.topology_dm parent = self._parent_mesh.topology_dm - swarm_parent_cell_nums = swarm.getField("DMSwarm_cellid") + swarm_parent_cell_nums = swarm.getField("DMSwarm_cellid").ravel() parent_renum = self._parent_mesh._dm_renumbering.getIndices() pStart, _ = parent.getChart() parent_renum_inv = np.empty_like(parent_renum) @@ -2063,7 +2063,7 @@ def cell_parent_cell_list(self): """Return a list of parent mesh cells numbers in vertex only mesh cell order. """ - cell_parent_cell_list = np.copy(self.topology_dm.getField("parentcellnum")) + cell_parent_cell_list = np.copy(self.topology_dm.getField("parentcellnum").ravel()) self.topology_dm.restoreField("parentcellnum") return cell_parent_cell_list[self.cell_closure[:, -1]] @@ -2082,7 +2082,7 @@ def cell_parent_base_cell_list(self): """ if not isinstance(self._parent_mesh, ExtrudedMeshTopology): raise AttributeError("Parent mesh is not extruded") - cell_parent_base_cell_list = np.copy(self.topology_dm.getField("parentcellbasenum")) + cell_parent_base_cell_list = np.copy(self.topology_dm.getField("parentcellbasenum").ravel()) self.topology_dm.restoreField("parentcellbasenum") return cell_parent_base_cell_list[self.cell_closure[:, -1]] @@ -2103,7 +2103,7 @@ def cell_parent_extrusion_height_list(self): """ if not isinstance(self._parent_mesh, ExtrudedMeshTopology): raise AttributeError("Parent mesh is not extruded.") - cell_parent_extrusion_height_list = np.copy(self.topology_dm.getField("parentcellextrusionheight")) + cell_parent_extrusion_height_list = np.copy(self.topology_dm.getField("parentcellextrusionheight").ravel()) self.topology_dm.restoreField("parentcellextrusionheight") return cell_parent_extrusion_height_list[self.cell_closure[:, -1]] @@ -2123,7 +2123,7 @@ def mark_entities(self, tf, label_value, label_name=None): @utils.cached_property # TODO: Recalculate if mesh moves def cell_global_index(self): """Return a list of unique cell IDs in vertex only mesh cell order.""" - cell_global_index = np.copy(self.topology_dm.getField("globalindex")) + cell_global_index = np.copy(self.topology_dm.getField("globalindex").ravel()) self.topology_dm.restoreField("globalindex") return cell_global_index @@ -2158,8 +2158,8 @@ def _make_input_ordering_sf(swarm, nroots, ilocal): # ilocal = None -> leaves are swarm points [0, 1, 2, ...). # ilocal can also be Firedrake cell numbers. sf = PETSc.SF().create(comm=swarm.comm) - input_ranks = swarm.getField("inputrank") - input_indices = swarm.getField("inputindex") + input_ranks = swarm.getField("inputrank").ravel() + input_indices = swarm.getField("inputindex").ravel() nleaves = len(input_ranks) if ilocal is not None and nleaves != len(ilocal): swarm.restoreField("inputrank") @@ -3652,7 +3652,7 @@ def _pic_swarm_in_mesh( # local_points[halo_indices] (it also updates local_points[~halo_indices]`, not changing any values there). # If some index of local_points_reduced corresponds to a missing point, local_points_reduced[index] is not updated # when we reduce and it does not update any leaf data, i.e., local_points, when we bcast. - owners = swarm.getField("DMSwarm_rank") + owners = swarm.getField("DMSwarm_rank").ravel() halo_indices, = np.where(owners != parent_mesh.comm.rank) halo_indices = halo_indices.astype(IntType) n = coords.shape[0] @@ -3868,13 +3868,13 @@ def _dmswarm_create( # NOTE ensure that swarm.restoreField is called for each field too! swarm_coords = swarm.getField("DMSwarmPIC_coor").reshape((num_vertices, gdim)) - swarm_parent_cell_nums = swarm.getField("DMSwarm_cellid") - field_parent_cell_nums = swarm.getField("parentcellnum") + swarm_parent_cell_nums = swarm.getField("DMSwarm_cellid").ravel() + field_parent_cell_nums = swarm.getField("parentcellnum").ravel() field_reference_coords = swarm.getField("refcoord").reshape((num_vertices, tdim)) - field_global_index = swarm.getField("globalindex") - field_rank = swarm.getField("DMSwarm_rank") - field_input_rank = swarm.getField("inputrank") - field_input_index = swarm.getField("inputindex") + field_global_index = swarm.getField("globalindex").ravel() + field_rank = swarm.getField("DMSwarm_rank").ravel() + field_input_rank = swarm.getField("inputrank").ravel() + field_input_index = swarm.getField("inputindex").ravel() swarm_coords[...] = coords swarm_parent_cell_nums[...] = plex_parent_cell_nums field_parent_cell_nums[...] = parent_cell_nums @@ -3895,8 +3895,8 @@ def _dmswarm_create( swarm.restoreField("DMSwarm_cellid") if extruded: - field_base_parent_cell_nums = swarm.getField("parentcellbasenum") - field_extrusion_heights = swarm.getField("parentcellextrusionheight") + field_base_parent_cell_nums = swarm.getField("parentcellbasenum").ravel() + field_extrusion_heights = swarm.getField("parentcellextrusionheight").ravel() field_base_parent_cell_nums[...] = base_parent_cell_nums field_extrusion_heights[...] = extrusion_heights swarm.restoreField("parentcellbasenum") diff --git a/scripts/firedrake-install b/scripts/firedrake-install index a83b04a451..1aae5fc3e2 100755 --- a/scripts/firedrake-install +++ b/scripts/firedrake-install @@ -742,7 +742,7 @@ def get_petsc_options(minimal=False): petsc_options.add("--download-suitesparse") petsc_options.add("--download-pastix") # Required by pastix - petsc_options.add("--download-hwloc") + petsc_options.update({"--download-hwloc", "--download-netlib-lapack", "--with-netlib-lapack-c-bindings"}) if osname == "Darwin": petsc_options.add("--download-hwloc-configure-arguments=--disable-opencl") # Serial mesh partitioner @@ -1517,6 +1517,7 @@ if mode == "install" or not args.update_script: if osname == "Darwin": package_manager = "brew" required_packages = ["gcc", + "pkg-config", "autoconf", "make", "automake", @@ -1535,6 +1536,7 @@ if mode == "install" or not args.update_script: "flex", # for ptscotch "cmake", "gfortran", + "pkg-config", "git", "libcurl4-openssl-dev", "pkgconf", # for p4est diff --git a/tests/firedrake/vertexonly/test_swarm.py b/tests/firedrake/vertexonly/test_swarm.py index 21c0b12eca..d3f5680041 100644 --- a/tests/firedrake/vertexonly/test_swarm.py +++ b/tests/firedrake/vertexonly/test_swarm.py @@ -208,7 +208,7 @@ def test_pic_swarm_in_mesh(parentmesh, redundant, exclude_halos): exclude_halos = True # Get point coords on current MPI rank - localpointcoords = np.copy(swarm.getField("DMSwarmPIC_coor")) + localpointcoords = np.copy(swarm.getField("DMSwarmPIC_coor").ravel()) swarm.restoreField("DMSwarmPIC_coor") if len(inputpointcoords.shape) > 1: localpointcoords = np.reshape(localpointcoords, (-1, inputpointcoords.shape[1])) @@ -218,11 +218,11 @@ def test_pic_swarm_in_mesh(parentmesh, redundant, exclude_halos): nptslocal = len(localpointcoords) nptsglobal = MPI.COMM_WORLD.allreduce(nptslocal, op=MPI.SUM) # Get parent PETSc cell indices on current MPI rank - localparentcellindices = np.copy(swarm.getField("DMSwarm_cellid")) + localparentcellindices = np.copy(swarm.getField("DMSwarm_cellid").ravel()) swarm.restoreField("DMSwarm_cellid") # also get the global coordinate numbering - globalindices = np.copy(swarm.getField("globalindex")) + globalindices = np.copy(swarm.getField("globalindex").ravel()) swarm.restoreField("globalindex") # Tests @@ -233,7 +233,7 @@ def test_pic_swarm_in_mesh(parentmesh, redundant, exclude_halos): # get custom fields on swarm - will fail if didn't get created for name, size, dtype in other_fields: - f = swarm.getField(name) + f = swarm.getField(name).ravel() assert len(f) == size*nptslocal assert f.dtype == dtype swarm.restoreField(name) @@ -330,7 +330,7 @@ def test_pic_swarm_in_mesh(parentmesh, redundant, exclude_halos): # Check that the rank numbering is correct. Since we know all points are at # the midpoints of cells, there should be no disagreement about cell # ownership and the voting algorithm should have no effect. - owned_ranks = np.copy(swarm.getField("DMSwarm_rank")) + owned_ranks = np.copy(swarm.getField("DMSwarm_rank").ravel()) swarm.restoreField("DMSwarm_rank") if exclude_halos: assert np.array_equal(owned_ranks, inputlocalpointcoordranks) @@ -339,7 +339,7 @@ def test_pic_swarm_in_mesh(parentmesh, redundant, exclude_halos): assert np.all(np.isin(inputlocalpointcoordranks, owned_ranks)) # check that the input rank is correct - input_ranks = np.copy(swarm.getField("inputrank")) + input_ranks = np.copy(swarm.getField("inputrank").ravel()) swarm.restoreField("inputrank") if exclude_halos: assert np.all(input_ranks == input_rank) @@ -351,7 +351,7 @@ def test_pic_swarm_in_mesh(parentmesh, redundant, exclude_halos): assert np.all(input_ranks < parentmesh.comm.size) # check that the input index is correct - input_indices = np.copy(swarm.getField("inputindex")) + input_indices = np.copy(swarm.getField("inputindex").ravel()) swarm.restoreField("inputindex") if exclude_halos: assert np.array_equal(input_indices, input_local_coord_indices) @@ -365,7 +365,7 @@ def test_pic_swarm_in_mesh(parentmesh, redundant, exclude_halos): # check we have unique parent cell numbers, which we should since we have # points at cell midpoints - parentcellnums = np.copy(swarm.getField("parentcellnum")) + parentcellnums = np.copy(swarm.getField("parentcellnum").ravel()) swarm.restoreField("parentcellnum") assert len(np.unique(parentcellnums)) == len(parentcellnums) @@ -378,7 +378,7 @@ def test_pic_swarm_in_mesh(parentmesh, redundant, exclude_halos): ): swarm.setPointCoordinates(localpointcoords, redundant=False, mode=PETSc.InsertMode.INSERT_VALUES) - petsclocalparentcellindices = np.copy(swarm.getField("DMSwarm_cellid")) + petsclocalparentcellindices = np.copy(swarm.getField("DMSwarm_cellid").ravel()) swarm.restoreField("DMSwarm_cellid") if exclude_halos: assert np.all(petsclocalparentcellindices == localparentcellindices) diff --git a/tests/firedrake/vertexonly/test_vertex_only_fs.py b/tests/firedrake/vertexonly/test_vertex_only_fs.py index a44791d0ea..a900aff568 100644 --- a/tests/firedrake/vertexonly/test_vertex_only_fs.py +++ b/tests/firedrake/vertexonly/test_vertex_only_fs.py @@ -116,7 +116,7 @@ def functionspace_tests(vm): h.dat.data_wo_with_halos[:] = -1 h.interpolate(g) # Exclude points which we know are missing - these should all be equal to -1 - input_ordering_parent_cell_nums = vm.input_ordering.topology_dm.getField("parentcellnum") + input_ordering_parent_cell_nums = vm.input_ordering.topology_dm.getField("parentcellnum").ravel() vm.input_ordering.topology_dm.restoreField("parentcellnum") idxs_to_include = input_ordering_parent_cell_nums != -1 assert np.allclose(h.dat.data_ro_with_halos[idxs_to_include], np.prod(vm.input_ordering.coordinates.dat.data_ro_with_halos[idxs_to_include].reshape(-1, vm.input_ordering.geometric_dimension()), axis=1)) @@ -221,7 +221,7 @@ def vectorfunctionspace_tests(vm): h.dat.data_wo_with_halos[:] = -1 h.interpolate(g) # Exclude points which we know are missing - these should all be equal to -1 - input_ordering_parent_cell_nums = vm.input_ordering.topology_dm.getField("parentcellnum") + input_ordering_parent_cell_nums = vm.input_ordering.topology_dm.getField("parentcellnum").ravel() vm.input_ordering.topology_dm.restoreField("parentcellnum") idxs_to_include = input_ordering_parent_cell_nums != -1 assert np.allclose(h.dat.data_ro[idxs_to_include], 2*vm.input_ordering.coordinates.dat.data_ro_with_halos[idxs_to_include]) diff --git a/tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py b/tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py index ac33b461de..8031173b53 100644 --- a/tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py +++ b/tests/firedrake/vertexonly/test_vertex_only_mesh_generation.py @@ -199,7 +199,7 @@ def verify_vertexonly_mesh(m, vm, inputvertexcoords, name): # Correct parent cell numbers stored_vertex_coords = np.copy(vm.topology_dm.getField("DMSwarmPIC_coor")).reshape((vm.num_cells(), gdim)) vm.topology_dm.restoreField("DMSwarmPIC_coor") - stored_parent_cell_nums = np.copy(vm.topology_dm.getField("parentcellnum")) + stored_parent_cell_nums = np.copy(vm.topology_dm.getField("parentcellnum").ravel()) vm.topology_dm.restoreField("parentcellnum") assert len(stored_vertex_coords) == len(stored_parent_cell_nums) if MPI.COMM_WORLD.size == 1: