From 226239896610e495b514e2f9b2cc2c76711ac9ad Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Wed, 2 Oct 2024 08:03:58 +0100 Subject: [PATCH] Pass graph partitioner to `mesh::refine` (#3444) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Remove refine option `distribute` and `ghost_mode` in favor of an explicit partitioner * Remove wrongly added transfer_matrix -> later PR * Add author name * Doc fix * Fix typi * Move maintain_coarse_partitioner to refine.cpp * Works better with the actual file * Add doc string to partitioner * Make doxygen happy * Switch to general partitioner naming * Small updates * Fix test failure * Simplifications * Small fix * Small fix * Work on docs * Simplify * Improve function name * Namespace fix * Compile fix * Add nanobind include * Small re-factor * Fixes * More updates * Fix broken tests * Doc update * Update test * Test updates * Work on docs * Lint * Minor doc updates --------- Co-authored-by: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Co-authored-by: Jørgen Schartum Dokken --- cpp/dolfinx/mesh/MeshTags.h | 4 +- cpp/dolfinx/refinement/refine.h | 100 +++++++------- cpp/dolfinx/refinement/utils.h | 95 +++++--------- cpp/test/mesh/refinement/interval.cpp | 12 +- cpp/test/mesh/refinement/rectangle.cpp | 6 +- python/dolfinx/mesh.py | 122 ++++++++++-------- python/dolfinx/wrappers/mesh.cpp | 96 +++++--------- python/dolfinx/wrappers/mesh.h | 59 +++++++++ python/dolfinx/wrappers/refinement.cpp | 45 +++++-- python/test/unit/refinement/test_interval.py | 12 +- .../test/unit/refinement/test_refinement.py | 113 ++++++++-------- 11 files changed, 330 insertions(+), 334 deletions(-) create mode 100644 python/dolfinx/wrappers/mesh.h diff --git a/cpp/dolfinx/mesh/MeshTags.h b/cpp/dolfinx/mesh/MeshTags.h index 32cb2de6a31..3026c258158 100644 --- a/cpp/dolfinx/mesh/MeshTags.h +++ b/cpp/dolfinx/mesh/MeshTags.h @@ -127,8 +127,8 @@ class MeshTags }; /// @brief Create MeshTags from arrays -/// @param[in] topology Mesh topology that the tags are associated with -/// @param[in] dim Topological dimension of tagged entities +/// @param[in] topology Mesh topology that the tags are associated with. +/// @param[in] dim Topological dimension of tagged entities. /// @param[in] entities Local vertex indices for tagged entities. /// @param[in] values Tag values for each entity in `entities`. The /// length of `values` must be equal to number of rows in `entities`. diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index 32f18d89827..430bcc495bc 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -1,4 +1,4 @@ -// Copyright (C) 2010-2023 Garth N. Wells +// Copyright (C) 2010-2024 Garth N. Wells and Paul T. Kühner // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -6,6 +6,7 @@ #pragma once +#include "dolfinx/graph/AdjacencyList.h" #include "dolfinx/mesh/Mesh.h" #include "dolfinx/mesh/Topology.h" #include "dolfinx/mesh/cell_types.h" @@ -15,85 +16,72 @@ #include #include #include +#include #include namespace dolfinx::refinement { -namespace impl -{ -/// @brief create the refined mesh by optionally redistributing it -template -mesh::Mesh -create_refined_mesh(const mesh::Mesh& mesh, - const graph::AdjacencyList& cell_adj, - std::span new_vertex_coords, - std::array xshape, bool redistribute, - mesh::GhostMode ghost_mode) -{ - if (dolfinx::MPI::size(mesh.comm()) == 1) - { - // No parallel construction necessary, i.e. redistribute also has no - // effect - return mesh::create_mesh(mesh.comm(), cell_adj.array(), - mesh.geometry().cmap(), new_vertex_coords, xshape, - ghost_mode); - } - else - { - // Let partition handle the parallel construction of the mesh - return partition(mesh, cell_adj, new_vertex_coords, xshape, redistribute, - ghost_mode); - } -} -} // namespace impl - -/// @brief Refine with markers, optionally redistributing, and -/// optionally calculating the parent-child relationships. +/// @brief Refine a mesh with markers. +/// +/// The refined mesh can be optionally re-partitioned across processes. +/// Passing `nullptr` for `partitioner`, refined cells will be on the +/// same process as the parent cell. /// -/// @param[in] mesh Input mesh to be refined -/// @param[in] edges Optional indices of the edges that should be split -/// by this refinement, if optional is not set, a uniform refinement is -/// performed, same behavior as passing a list of all indices. -/// @param[in] redistribute Flag to call the Mesh Partitioner to -/// redistribute after refinement. -/// @param[in] ghost_mode Ghost mode of the refined mesh. +/// Parent-child relationships can be optionally computed. Parent-child +/// relationships can be used to create MeshTags on the refined mesh +/// from MeshTags on the parent mesh. +/// +/// @warning Using the default partitioner for a refined mesh, the +/// refined mesh will **not** include ghosts cells (cells connected by +/// facet to an owned cell) even if the parent mesh is ghosted. +/// +/// @warning Passing `nullptr` for `partitioner`, the refined mesh will +/// **not** have ghosts cells even if the parent mesh is ghosted. The +/// possibility to not re-partition the refined mesh and include ghost +/// cells in the refined mesh will be added in a future release. +/// +/// @param[in] mesh Input mesh to be refined. +/// @param[in] edges Indices of the edges that should be split in the +/// refinement. If not provided (`std::nullopt`), uniform refinement is +/// performed. +/// @param[in] partitioner Partitioner to be used to distribute the +/// refined mesh. If not callable, refined cells will be on the same +/// process as the parent cell. /// @param[in] option Control the computation of parent facets, parent -/// cells. If an option is unselected, an empty list is returned. -/// @return New Mesh and optional parent cell index, parent facet +/// cells. If an option is not selected, an empty list is returned. +/// @return New mesh, and optional parent cell indices and parent facet /// indices. template std::tuple, std::optional>, std::optional>> refine(const mesh::Mesh& mesh, - std::optional> edges, bool redistribute, - mesh::GhostMode ghost_mode = mesh::GhostMode::shared_facet, + std::optional> edges, + const mesh::CellPartitionFunction& partitioner + = mesh::create_cell_partitioner(mesh::GhostMode::none), Option option = Option::none) { auto topology = mesh.topology(); assert(topology); - - mesh::CellType cell_t = topology->cell_type(); - if (!mesh::is_simplex(cell_t)) + if (!mesh::is_simplex(topology->cell_type())) throw std::runtime_error("Refinement only defined for simplices"); - bool oned = topology->cell_type() == mesh::CellType::interval; + auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet] - = oned ? interval::compute_refinement_data(mesh, edges, option) - : plaza::compute_refinement_data(mesh, edges, option); + = (topology->cell_type() == mesh::CellType::interval) + ? interval::compute_refinement_data(mesh, edges, option) + : plaza::compute_refinement_data(mesh, edges, option); - mesh::Mesh refined_mesh = impl::create_refined_mesh( - mesh, cell_adj, std::span(new_vertex_coords), xshape, - redistribute, ghost_mode); + mesh::Mesh mesh1 = mesh::create_mesh( + mesh.comm(), mesh.comm(), cell_adj.array(), mesh.geometry().cmap(), + mesh.comm(), new_vertex_coords, xshape, partitioner); - // Report the number of refined cellse + // Report the number of refined cells const int D = topology->dim(); const std::int64_t n0 = topology->index_map(D)->size_global(); - const std::int64_t n1 = refined_mesh.topology()->index_map(D)->size_global(); + const std::int64_t n1 = mesh1.topology()->index_map(D)->size_global(); spdlog::info( "Number of cells increased from {} to {} ({}% increase).", n0, n1, 100.0 * (static_cast(n1) / static_cast(n0) - 1.0)); - return {std::move(refined_mesh), std::move(parent_cell), - std::move(parent_facet)}; + return {std::move(mesh1), std::move(parent_cell), std::move(parent_facet)}; } - } // namespace dolfinx::refinement diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h index 275b40d7cc2..ec550cb3575 100644 --- a/cpp/dolfinx/refinement/utils.h +++ b/cpp/dolfinx/refinement/utils.h @@ -264,85 +264,48 @@ create_new_vertices(MPI_Comm comm, xshape}; } -/// Use vertex and topology data to partition new mesh across -/// processes -/// @param[in] old_mesh -/// @param[in] cell_topology Topology of cells, (vertex indices) -/// @param[in] new_coords New coordinates, row-major storage -/// @param[in] xshape The shape of `new_coords` -/// @param[in] redistribute Call graph partitioner if true -/// @param[in] ghost_mode None or shared_facet -/// @return New mesh -template -mesh::Mesh partition(const mesh::Mesh& old_mesh, - const graph::AdjacencyList& cell_topology, - std::span new_coords, - std::array xshape, bool redistribute, - mesh::GhostMode ghost_mode) -{ - if (redistribute) - { - return mesh::create_mesh(old_mesh.comm(), cell_topology.array(), - old_mesh.geometry().cmap(), new_coords, xshape, - ghost_mode); - } - else - { - auto partitioner - = [](MPI_Comm comm, int, const std::vector& cell_types, - const std::vector>& cell_topology) - { - std::int32_t mpi_rank = MPI::rank(comm); - std::int32_t num_cell_vertices - = mesh::num_cell_vertices(cell_types.front()); - std::int32_t num_cells = cell_topology.front().size() / num_cell_vertices; - std::vector destinations(num_cells, mpi_rank); - std::vector dest_offsets(num_cells + 1); - std::iota(dest_offsets.begin(), dest_offsets.end(), 0); - return graph::AdjacencyList(std::move(destinations), - std::move(dest_offsets)); - }; - - return mesh::create_mesh(old_mesh.comm(), old_mesh.comm(), - cell_topology.array(), old_mesh.geometry().cmap(), - old_mesh.comm(), new_coords, xshape, partitioner); - } -} - -/// @brief Given an index map, add "n" extra indices at the end of local range +/// @todo Improve docstring. +/// +/// @brief Given an index map, add "n" extra indices at the end of local +/// range. /// -/// @note The returned global indices (local and ghosts) are adjust for the -/// addition of new local indices. -/// @param[in] map Index map -/// @param[in] n Number of new local indices -/// @return New global indices for both owned and ghosted indices in input -/// index map. +/// @note The returned global indices (local and ghosts) are adjust for +/// the addition of new local indices. +/// +/// @param[in] map Index map. +/// @param[in] n Number of new local indices. +/// @return New global indices for both owned and ghosted indices in +/// input index map. std::vector adjust_indices(const common::IndexMap& map, std::int32_t n); -/// @brief Transfer facet MeshTags from coarse mesh to refined mesh -/// @note The refined mesh must not have been redistributed during -/// refinement -/// @note GhostMode must be GhostMode.none -/// @param[in] tags0 Tags on the parent mesh -/// @param[in] topology1 Refined mesh topology +/// @brief Transfer facet MeshTags from coarse mesh to refined mesh. +/// +/// @warning The refined mesh must not have been redistributed during +/// refinement. +/// +/// @warning Mesh/topology `GhostMode` must be mesh::GhostMode::none. +/// +/// @param[in] tags0 Tags on the parent mesh. +/// @param[in] topology1 Refined mesh topology. /// @param[in] cell Parent cell of each cell in refined mesh -/// @param[in] facet Local facets of parent in each cell in refined mesh -/// @return (0) entities and (1) values on the refined topology +/// @param[in] facet Local facets of parent in each cell in refined mesh. +/// @return (0) entities and (1) values on the refined topology. std::array, 2> transfer_facet_meshtag( const mesh::MeshTags& tags0, const mesh::Topology& topology1, std::span cell, std::span facet); /// @brief Transfer cell MeshTags from coarse mesh to refined mesh. /// -/// @note The refined mesh must not have been redistributed during +/// @warning The refined mesh must not have been redistributed during /// refinement. -/// @note GhostMode must be GhostMode.none /// -/// @param[in] tags0 Tags on the parent mesh -/// @param[in] topology1 Refined mesh topology -/// @param[in] parent_cell Parent cell of each cell in refined mesh -/// @return (0) entities and (1) values on the refined topology +/// @warning Mesh/topology `GhostMode` must be mesh::GhostMode::none. +/// +/// @param[in] tags0 Tags on the parent mesh. +/// @param[in] topology1 Refined mesh topology. +/// @param[in] parent_cell Parent cell of each cell in refined mesh. +/// @return (0) entities and (1) values on the refined topology. std::array, 2> transfer_cell_meshtag(const mesh::MeshTags& tags0, const mesh::Topology& topology1, diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp index 763418ac4b3..1253fe00b04 100644 --- a/cpp/test/mesh/refinement/interval.cpp +++ b/cpp/test/mesh/refinement/interval.cpp @@ -73,8 +73,7 @@ TEMPLATE_TEST_CASE("Interval uniform refinement", // TODO: parent_facet auto [refined_mesh, parent_edge, parent_facet] = refinement::refine( - mesh, std::nullopt, false, mesh::GhostMode::shared_facet, - refinement::Option::parent_cell); + mesh, std::nullopt, nullptr, refinement::Option::parent_cell); std::vector expected_x = { /* v_0 */ 0.0, 0.0, 0.0, @@ -114,7 +113,8 @@ TEMPLATE_TEST_CASE("Interval adaptive refinement", std::vector edges{1}; // TODO: parent_facet auto [refined_mesh, parent_edge, parent_facet] = refinement::refine( - mesh, std::span(edges), false, mesh::GhostMode::shared_facet, + mesh, std::span(edges), + mesh::create_cell_partitioner(mesh::GhostMode::shared_facet), refinement::Option::parent_cell); std::vector expected_x = { @@ -190,15 +190,13 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)", mesh::Mesh mesh = create_mesh(); mesh.topology()->create_connectivity(1, 0); - // TODO: parent_facet auto [refined_mesh, parent_edges, parent_facet] = refinement::refine( - mesh, std::nullopt, false, mesh::GhostMode::shared_facet, - refinement::Option::parent_cell); + mesh, std::nullopt, nullptr, refinement::Option::parent_cell); T rank_d = static_cast(rank); T comm_size_d = static_cast(comm_size); - auto x = refined_mesh.geometry().x(); + std::span x = refined_mesh.geometry().x(); std::ranges::sort(x); std::vector expected_x = {rank_d / comm_size_d, diff --git a/cpp/test/mesh/refinement/rectangle.cpp b/cpp/test/mesh/refinement/rectangle.cpp index 954879abb1d..ff2a67979df 100644 --- a/cpp/test/mesh/refinement/rectangle.cpp +++ b/cpp/test/mesh/refinement/rectangle.cpp @@ -99,9 +99,9 @@ plotter.show() // plaza requires the edges to be pre initialized! mesh.topology()->create_entities(1); - auto [mesh_fine, parent_cell, parent_facet] - = refinement::refine(mesh, std::nullopt, false, mesh::GhostMode::none, - refinement::Option::parent_cell_and_facet); + auto [mesh_fine, parent_cell, parent_facet] = refinement::refine( + mesh, std::nullopt, mesh::create_cell_partitioner(mesh::GhostMode::none), + refinement::Option::parent_cell_and_facet); // vertex layout: // 8---7---5 diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index 28b9eb759ea..2203ca8d961 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -276,11 +276,11 @@ class Mesh: _geometry: Geometry _ufl_domain: typing.Optional[ufl.Mesh] - def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]): + def __init__(self, msh, domain: typing.Optional[ufl.Mesh]): """Initialize mesh from a C++ mesh. Args: - mesh: A C++ mesh object. + msh: A C++ mesh object. domain: A UFL domain. Note: @@ -289,7 +289,7 @@ def __init__(self, mesh, domain: typing.Optional[ufl.Mesh]): This class is combined with different base classes that depend on the scalar type used in the Mesh. """ - self._cpp_object = mesh + self._cpp_object = msh self._topology = Topology(self._cpp_object.topology) self._geometry = Geometry(self._cpp_object.geometry) self._ufl_domain = domain @@ -439,25 +439,25 @@ def compute_incident_entities( return _cpp.mesh.compute_incident_entities(topology._cpp_object, entities, d0, d1) -def compute_midpoints(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32]): +def compute_midpoints(msh: Mesh, dim: int, entities: npt.NDArray[np.int32]): """Compute the midpoints of a set of mesh entities. Args: - mesh: The mesh. + msh: The mesh. dim: Topological dimension of the mesh entities to consider. entities: Indices of entities in ``mesh`` to consider. Returns: Midpoints of the entities, shape ``(num_entities, 3)``. """ - return _cpp.mesh.compute_midpoints(mesh._cpp_object, dim, entities) + return _cpp.mesh.compute_midpoints(msh._cpp_object, dim, entities) -def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray: +def locate_entities(msh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray: """Compute mesh entities satisfying a geometric marking function. Args: - mesh: Mesh to locate entities on. + msh: Mesh to locate entities on. dim: Topological dimension of the mesh entities to consider. marker: A function that takes an array of points ``x`` with shape ``(gdim, num_points)`` and returns an array of @@ -467,10 +467,10 @@ def locate_entities(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray Returns: Indices (local to the process) of marked mesh entities. """ - return _cpp.mesh.locate_entities(mesh._cpp_object, dim, marker) + return _cpp.mesh.locate_entities(msh._cpp_object, dim, marker) -def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray: +def locate_entities_boundary(msh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray: """Compute mesh entities that are connected to an owned boundary facet and satisfy a geometric marking function. @@ -484,7 +484,7 @@ def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> n parallel communication. Args: - mesh: Mesh to locate boundary entities on. + msh: Mesh to locate boundary entities on. dim: Topological dimension of the mesh entities to consider marker: Function that takes an array of points ``x`` with shape ``(gdim, num_points)`` and returns an array of booleans of @@ -494,12 +494,12 @@ def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> n Returns: Indices (local to the process) of marked mesh entities. """ - return _cpp.mesh.locate_entities_boundary(mesh._cpp_object, dim, marker) + return _cpp.mesh.locate_entities_boundary(msh._cpp_object, dim, marker) def transfer_meshtag( meshtag: MeshTags, - mesh1: Mesh, + msh1: Mesh, parent_cell: npt.NDArray[np.int32], parent_facet: typing.Optional[npt.NDArray[np.int8]] = None, ) -> MeshTags: @@ -507,7 +507,7 @@ def transfer_meshtag( Args: meshtag: Mesh tags on the coarse, parent mesh. - mesh1: The refined mesh. + msh1: The refined mesh. parent_cell: Index of the parent cell for each cell in the refined mesh. parent_facet: Index of the local parent facet for each cell @@ -519,13 +519,13 @@ def transfer_meshtag( """ if meshtag.dim == meshtag.topology.dim: mt = _cpp.refinement.transfer_cell_meshtag( - meshtag._cpp_object, mesh1.topology._cpp_object, parent_cell + meshtag._cpp_object, msh1.topology._cpp_object, parent_cell ) return MeshTags(mt) elif meshtag.dim == meshtag.topology.dim - 1: assert parent_facet is not None mt = _cpp.refinement.transfer_facet_meshtag( - meshtag._cpp_object, mesh1.topology._cpp_object, parent_cell, parent_facet + meshtag._cpp_object, msh1.topology._cpp_object, parent_cell, parent_facet ) return MeshTags(mt) else: @@ -533,32 +533,47 @@ def transfer_meshtag( def refine( - mesh: Mesh, + msh: Mesh, edges: typing.Optional[np.ndarray] = None, - redistribute: bool = True, - ghost_mode: GhostMode = GhostMode.shared_facet, + partitioner: typing.Optional[typing.Callable] = create_cell_partitioner(GhostMode.none), option: RefinementOption = RefinementOption.none, ) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]: """Refine a mesh. + Passing ``None`` for ``partitioner``, refined cells will be on the + same process as the parent cell. + + Note: + Using the default partitioner for the refined mesh, the refined + mesh will **not** include ghosts cells (cells connected by facet + to an owned cells) even if the parent mesh is ghosted. + + Warning: + Passing ``None`` for ``partitioner``, the refined mesh will + **not** have ghosts cells even if the parent mesh has ghost + cells. The possibility to not re-partition the refined mesh and + include ghost cells in the refined mesh will be added in a + future release. + Args: - mesh: Mesh from which to create the refined mesh. + msh: Mesh from which to create the refined mesh. edges: Indices of edges to split during refinement. If ``None``, mesh refinement is uniform. - redistribute: - Refined mesh is re-partitioned if ``True``. - ghost_mode: Ghost mode to use for the refined mesh. - option: Controls whether parent cells and/or parent facets are computed. - + partitioner: Partitioner to distribute the refined mesh. If + ``None`` no redistribution is performed, i.e. refined cells + remain on the same process as the parent cell. + option: Controls whether parent cells and/or parent facets are + computed. Returns: Refined mesh, (optional) parent cells, (optional) parent facets """ mesh1, parent_cell, parent_facet = _cpp.refinement.refine( - mesh._cpp_object, edges, redistribute, ghost_mode, option + msh._cpp_object, edges, partitioner, option ) - # Create new ufl domain as it will carry a reference to the C++ mesh in the ufl_cargo - ufl_domain = ufl.Mesh(mesh._ufl_domain.ufl_coordinate_element()) # type: ignore + # Create new ufl domain as it will carry a reference to the C++ mesh + # in the ufl_cargo + ufl_domain = ufl.Mesh(msh._ufl_domain.ufl_coordinate_element()) # type: ignore return Mesh(mesh1, ufl_domain), parent_cell, parent_facet @@ -582,8 +597,8 @@ def create_mesh( x: Mesh geometry ('node' coordinates), with shape ``(num_nodes, gdim)``. e: UFL mesh. The mesh scalar type is determined by the scalar type of ``e``. - partitioner: Function that computes the parallel distribution of - cells across MPI ranks. + partitioner: Function that determines the parallel distribution + of cells across MPI ranks. Note: If required, the coordinates ``x`` will be cast to the same @@ -593,7 +608,7 @@ def create_mesh( A mesh. """ if partitioner is None and comm.size > 1: - partitioner = _cpp.mesh.create_cell_partitioner(GhostMode.none) + partitioner = create_cell_partitioner(GhostMode.none) x = np.asarray(x, order="C") if x.ndim == 1: @@ -635,9 +650,9 @@ def create_mesh( x = np.asarray(x, dtype=dtype, order="C") cells = np.asarray(cells, dtype=np.int64, order="C") - mesh = _cpp.mesh.create_mesh(comm, cells, cmap._cpp_object, x, partitioner) + msh = _cpp.mesh.create_mesh(comm, cells, cmap._cpp_object, x, partitioner) - return Mesh(mesh, domain) + return Mesh(msh, domain) def create_submesh( @@ -646,7 +661,7 @@ def create_submesh( """Create a mesh with specified entities from another mesh. Args: - mesh: Mesh to create the sub-mesh from. + msh: Mesh to create the sub-mesh from. dim: Topological dimension of the entities in ``msh`` to include in the sub-mesh. entities: Indices of entities in ``msh`` to include in the @@ -675,7 +690,7 @@ def create_submesh( def meshtags( - mesh: Mesh, + msh: Mesh, dim: int, entities: npt.NDArray[np.int32], values: typing.Union[np.ndarray, int, float], @@ -683,7 +698,7 @@ def meshtags( """Create a MeshTags object that associates data with a subset of mesh entities. Args: - mesh: The mesh. + msh: The mesh. dim: Topological dimension of the mesh entity. entities: Indices(local to process) of entities to associate values with . The array must be sorted and must not contain @@ -716,19 +731,19 @@ def meshtags( raise NotImplementedError(f"Type {values.dtype} not supported.") return MeshTags( - ftype(mesh.topology._cpp_object, dim, np.asarray(entities, dtype=np.int32), values) + ftype(msh.topology._cpp_object, dim, np.asarray(entities, dtype=np.int32), values) ) def meshtags_from_entities( - mesh: Mesh, dim: int, entities: _cpp.graph.AdjacencyList_int32, values: npt.NDArray[typing.Any] + msh: Mesh, dim: int, entities: _cpp.graph.AdjacencyList_int32, values: npt.NDArray[typing.Any] ): """Create a :class:dolfinx.mesh.MeshTags` object that associates data with a subset of mesh entities, where the entities are defined by their vertices. Args: - mesh: The mesh. + msh: The mesh. dim: Topological dimension of the mesh entity. entities: Entities to associated values with, with entities defined by their vertices. @@ -748,7 +763,7 @@ def meshtags_from_entities( elif isinstance(values, float): values = np.full(entities.num_nodes, values, dtype=np.double) values = np.asarray(values) - return MeshTags(_cpp.mesh.create_meshtags(mesh.topology._cpp_object, dim, entities, values)) + return MeshTags(_cpp.mesh.create_meshtags(msh.topology._cpp_object, dim, entities, values)) def create_interval( @@ -779,12 +794,13 @@ def create_interval( partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode) domain = ufl.Mesh(basix.ufl.element("Lagrange", "interval", 1, shape=(1,), dtype=dtype)) # type: ignore if np.issubdtype(dtype, np.float32): - mesh = _cpp.mesh.create_interval_float32(comm, nx, points, ghost_mode, partitioner) + msh = _cpp.mesh.create_interval_float32(comm, nx, points, ghost_mode, partitioner) elif np.issubdtype(dtype, np.float64): - mesh = _cpp.mesh.create_interval_float64(comm, nx, points, ghost_mode, partitioner) + msh = _cpp.mesh.create_interval_float64(comm, nx, points, ghost_mode, partitioner) else: raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}") - return Mesh(mesh, domain) + + return Mesh(msh, domain) def create_unit_interval( @@ -847,12 +863,13 @@ def create_rectangle( partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode) domain = ufl.Mesh(basix.ufl.element("Lagrange", cell_type.name, 1, shape=(2,), dtype=dtype)) # type: ignore if np.issubdtype(dtype, np.float32): - mesh = _cpp.mesh.create_rectangle_float32(comm, points, n, cell_type, partitioner, diagonal) + msh = _cpp.mesh.create_rectangle_float32(comm, points, n, cell_type, partitioner, diagonal) elif np.issubdtype(dtype, np.float64): - mesh = _cpp.mesh.create_rectangle_float64(comm, points, n, cell_type, partitioner, diagonal) + msh = _cpp.mesh.create_rectangle_float64(comm, points, n, cell_type, partitioner, diagonal) else: raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}") - return Mesh(mesh, domain) + + return Mesh(msh, domain) def create_unit_square( @@ -925,12 +942,13 @@ def create_box( partitioner = _cpp.mesh.create_cell_partitioner(ghost_mode) domain = ufl.Mesh(basix.ufl.element("Lagrange", cell_type.name, 1, shape=(3,), dtype=dtype)) # type: ignore if np.issubdtype(dtype, np.float32): - mesh = _cpp.mesh.create_box_float32(comm, points, n, cell_type, partitioner) + msh = _cpp.mesh.create_box_float32(comm, points, n, cell_type, partitioner) elif np.issubdtype(dtype, np.float64): - mesh = _cpp.mesh.create_box_float64(comm, points, n, cell_type, partitioner) + msh = _cpp.mesh.create_box_float64(comm, points, n, cell_type, partitioner) else: raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}") - return Mesh(mesh, domain) + + return Mesh(msh, domain) def create_unit_cube( @@ -973,12 +991,12 @@ def create_unit_cube( def entities_to_geometry( - mesh: Mesh, dim: int, entities: npt.NDArray[np.int32], permute=False + msh: Mesh, dim: int, entities: npt.NDArray[np.int32], permute=False ) -> npt.NDArray[np.int32]: """Compute the geometric DOFs associated with the closure of the given mesh entities. Args: - mesh: The mesh. + msh: The mesh. dim: Topological dimension of the entities of interest. entities: Entity indices (local to the process). permute: Permute the DOFs such that they are consistent with the @@ -989,7 +1007,7 @@ def entities_to_geometry( The geometric DOFs associated with the closure of the entities in `entities`. """ - return _cpp.mesh.entities_to_geometry(mesh._cpp_object, dim, entities, permute) + return _cpp.mesh.entities_to_geometry(msh._cpp_object, dim, entities, permute) def exterior_facet_indices(topology: Topology) -> npt.NDArray[np.int32]: diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index f8c3679974d..5cf66c9c6ab 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -4,11 +4,11 @@ // // SPDX-License-Identifier: LGPL-3.0-or-later +#include "mesh.h" #include "MPICommWrapper.h" #include "array.h" #include "caster_mpi.h" #include "numpy_dtype.h" -#include #include #include #include @@ -21,7 +21,6 @@ #include #include #include -#include #include #include #include @@ -36,49 +35,10 @@ namespace nb = nanobind; -namespace -{ -/// Wrap a Python graph partitioning function as a C++ function -template -auto create_partitioner_cpp(Functor p) +namespace dolfinx_wrappers { - return [p](MPI_Comm comm, int nparts, - const dolfinx::graph::AdjacencyList& local_graph, - bool ghosting) - { - return p(dolfinx_wrappers::MPICommWrapper(comm), nparts, local_graph, - ghosting); - }; -} - -/// Wrap a C++ cell partitioning function as a Python function -template -auto create_cell_partitioner_py(Functor p) +namespace part::impl { - return [p](dolfinx_wrappers::MPICommWrapper comm, int n, - const std::vector& cell_types, - std::vector> cells_nb) - { - std::vector> cells; - std::ranges::transform( - cells_nb, std::back_inserter(cells), [](auto c) - { return std::span(c.data(), c.size()); }); - return p(comm.get(), n, cell_types, cells); - }; -} - -using PythonCellPartitionFunction - = std::function( - dolfinx_wrappers::MPICommWrapper, int, - const std::vector&, - std::vector>)>; - -using CppCellPartitionFunction - = std::function( - MPI_Comm, int, const std::vector& q, - const std::vector>&)>; - -/// Wrap a Python cell graph partitioning function as a C++ function CppCellPartitionFunction create_cell_partitioner_cpp(const PythonCellPartitionFunction& p) { @@ -103,10 +63,8 @@ create_cell_partitioner_cpp(const PythonCellPartitionFunction& p) else return nullptr; } -} // namespace +} // namespace part::impl -namespace dolfinx_wrappers -{ template void declare_meshtags(nb::module_& m, std::string type) { @@ -260,7 +218,8 @@ void declare_mesh(nb::module_& m, std::string type) "__init__", [](dolfinx::mesh::Mesh* mesh, MPICommWrapper comm, std::shared_ptr topology, - dolfinx::mesh::Geometry& geometry) { + dolfinx::mesh::Geometry& geometry) + { new (mesh) dolfinx::mesh::Mesh(comm.get(), topology, geometry); }, nb::arg("comm"), nb::arg("topology"), nb::arg("geometry")) @@ -279,11 +238,12 @@ void declare_mesh(nb::module_& m, std::string type) m.def( create_interval.c_str(), [](MPICommWrapper comm, std::int64_t n, std::array p, - dolfinx::mesh::GhostMode ghost_mode, - const PythonCellPartitionFunction& part) + dolfinx::mesh::GhostMode mode, + const part::impl::PythonCellPartitionFunction& part) { return dolfinx::mesh::create_interval( - comm.get(), n, p, ghost_mode, create_cell_partitioner_cpp(part)); + comm.get(), n, p, mode, + part::impl::create_cell_partitioner_cpp(part)); }, nb::arg("comm"), nb::arg("n"), nb::arg("p"), nb::arg("ghost_mode"), nb::arg("partitioner").none()); @@ -293,12 +253,12 @@ void declare_mesh(nb::module_& m, std::string type) create_rectangle.c_str(), [](MPICommWrapper comm, std::array, 2> p, std::array n, dolfinx::mesh::CellType celltype, - const PythonCellPartitionFunction& part, + const part::impl::PythonCellPartitionFunction& part, dolfinx::mesh::DiagonalType diagonal) { return dolfinx::mesh::create_rectangle( - comm.get(), p, n, celltype, create_cell_partitioner_cpp(part), - diagonal); + comm.get(), p, n, celltype, + part::impl::create_cell_partitioner_cpp(part), diagonal); }, nb::arg("comm"), nb::arg("p"), nb::arg("n"), nb::arg("celltype"), nb::arg("partitioner").none(), nb::arg("diagonal")); @@ -308,11 +268,12 @@ void declare_mesh(nb::module_& m, std::string type) create_box.c_str(), [](MPICommWrapper comm, std::array, 2> p, std::array n, dolfinx::mesh::CellType celltype, - const PythonCellPartitionFunction& part) + const part::impl::PythonCellPartitionFunction& part) { MPI_Comm _comm = comm.get(); - return dolfinx::mesh::create_box(_comm, _comm, p, n, celltype, - create_cell_partitioner_cpp(part)); + return dolfinx::mesh::create_box( + _comm, _comm, p, n, celltype, + part::impl::create_cell_partitioner_cpp(part)); }, nb::arg("comm"), nb::arg("p"), nb::arg("n"), nb::arg("celltype"), nb::arg("partitioner").none()); @@ -323,7 +284,7 @@ void declare_mesh(nb::module_& m, std::string type) nb::c_contig>>& cells_nb, const std::vector>& elements, nb::ndarray x, - const PythonCellPartitionFunction& p) + const part::impl::PythonCellPartitionFunction& p) { std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape(1); @@ -365,7 +326,7 @@ void declare_mesh(nb::module_& m, std::string type) nb::ndarray, nb::c_contig> cells, const dolfinx::fem::CoordinateElement& element, nb::ndarray x, - const PythonCellPartitionFunction& p) + const part::impl::PythonCellPartitionFunction& p) { std::size_t shape1 = x.ndim() == 1 ? 1 : x.shape(1); if (p) @@ -599,7 +560,8 @@ void mesh(nb::module_& m) m.def( "compute_entities", [](MPICommWrapper comm, const dolfinx::mesh::Topology& topology, int dim, - int index) { + int index) + { return dolfinx::mesh::compute_entities(comm.get(), topology, dim, index); }, @@ -721,8 +683,6 @@ void mesh(nb::module_& m) ghost_owners_span, boundary_vertices_span); }); - // dolfinx::mesh::MeshTags - declare_meshtags(m, "int8"); declare_meshtags(m, "int32"); declare_meshtags(m, "int64"); @@ -733,10 +693,11 @@ void mesh(nb::module_& m) m.def( "create_cell_partitioner", - [](dolfinx::mesh::GhostMode gm) -> PythonCellPartitionFunction + [](dolfinx::mesh::GhostMode mode) + -> part::impl::PythonCellPartitionFunction { - return create_cell_partitioner_py( - dolfinx::mesh::create_cell_partitioner(gm)); + return part::impl::create_cell_partitioner_py( + dolfinx::mesh::create_cell_partitioner(mode)); }, "Create default cell partitioner."); m.def( @@ -746,11 +707,12 @@ void mesh(nb::module_& m) const dolfinx::graph::AdjacencyList& local_graph, bool ghosting)> part, - dolfinx::mesh::GhostMode ghost_mode) -> PythonCellPartitionFunction + dolfinx::mesh::GhostMode mode) + -> part::impl::PythonCellPartitionFunction { - return create_cell_partitioner_py( + return part::impl::create_cell_partitioner_py( dolfinx::mesh::create_cell_partitioner( - ghost_mode, create_partitioner_cpp(part))); + mode, part::impl::create_partitioner_cpp(part))); }, nb::arg("part"), nb::arg("ghost_mode") = dolfinx::mesh::GhostMode::none, "Create a cell partitioner from a graph partitioning function."); diff --git a/python/dolfinx/wrappers/mesh.h b/python/dolfinx/wrappers/mesh.h new file mode 100644 index 00000000000..5471dc1ab24 --- /dev/null +++ b/python/dolfinx/wrappers/mesh.h @@ -0,0 +1,59 @@ +// Copyright (C) 2017-2024 Chris N. Richardson and Garth N. Wells +// +// This file is part of DOLFINx (https://www.fenicsproject.org) +// +// SPDX-License-Identifier: LGPL-3.0-or-later + +#include "MPICommWrapper.h" +#include +#include +#include + +namespace nb = nanobind; + +namespace dolfinx_wrappers::part::impl +{ +/// Wrap a Python graph partitioning function as a C++ function +template +auto create_partitioner_cpp(Functor p) +{ + return [p](MPI_Comm comm, int nparts, + const dolfinx::graph::AdjacencyList& local_graph, + bool ghosting) + { + return p(dolfinx_wrappers::MPICommWrapper(comm), nparts, local_graph, + ghosting); + }; +} + +/// Wrap a C++ cell partitioning function as a Python function +template +auto create_cell_partitioner_py(Functor p) +{ + return [p](dolfinx_wrappers::MPICommWrapper comm, int n, + const std::vector& cell_types, + std::vector> cells_nb) + { + std::vector> cells; + std::ranges::transform( + cells_nb, std::back_inserter(cells), [](auto c) + { return std::span(c.data(), c.size()); }); + return p(comm.get(), n, cell_types, cells); + }; +} + +using PythonCellPartitionFunction + = std::function( + dolfinx_wrappers::MPICommWrapper, int, + const std::vector&, + std::vector>)>; + +using CppCellPartitionFunction + = std::function( + MPI_Comm, int, const std::vector& q, + const std::vector>&)>; + +/// Wrap a Python cell graph partitioning function as a C++ function +CppCellPartitionFunction +create_cell_partitioner_cpp(const PythonCellPartitionFunction& p); +} // namespace dolfinx_wrappers::part::impl diff --git a/python/dolfinx/wrappers/refinement.cpp b/python/dolfinx/wrappers/refinement.cpp index 37ca3937992..d9b86fcf781 100644 --- a/python/dolfinx/wrappers/refinement.cpp +++ b/python/dolfinx/wrappers/refinement.cpp @@ -1,21 +1,24 @@ -// Copyright (C) 2018 Chris N. Richardson and Garth N. Wells +// Copyright (C) 2018-2024 Chris N. Richardson and Garth N. Wells // // This file is part of DOLFINx (https://www.fenicsproject.org) // // SPDX-License-Identifier: LGPL-3.0-or-later +#include "MPICommWrapper.h" #include "array.h" +#include "mesh.h" #include #include #include #include #include #include -#include #include #include +#include #include #include +#include #include #include #include @@ -25,11 +28,10 @@ namespace nb = nanobind; -namespace dolfinx_wrappers +namespace { - template -void export_refinement_with_variable_mesh_type(nb::module_& m) +void export_refinement(nb::module_& m) { m.def( "refine", @@ -37,7 +39,9 @@ void export_refinement_with_variable_mesh_type(nb::module_& m) std::optional< nb::ndarray, nb::c_contig>> edges, - bool redistribute, dolfinx::mesh::GhostMode ghost_mode, + std::optional< + dolfinx_wrappers::part::impl::PythonCellPartitionFunction> + partitioner, dolfinx::refinement::Option option) { std::optional> cpp_edges(std::nullopt); @@ -47,30 +51,45 @@ void export_refinement_with_variable_mesh_type(nb::module_& m) std::span(edges.value().data(), edges.value().size())); } + dolfinx_wrappers::part::impl::CppCellPartitionFunction cpp_partitioner + = partitioner.has_value() + ? dolfinx_wrappers::part::impl::create_cell_partitioner_cpp( + partitioner.value()) + : nullptr; auto [mesh1, cell, facet] = dolfinx::refinement::refine( - mesh, cpp_edges, redistribute, ghost_mode, option); + mesh, cpp_edges, cpp_partitioner, option); std::optional> python_cell( std::nullopt); if (cell.has_value()) - python_cell.emplace(as_nbarray(std::move(cell.value()))); + { + python_cell.emplace( + dolfinx_wrappers::as_nbarray(std::move(cell.value()))); + } std::optional> python_facet( std::nullopt); if (facet.has_value()) - python_facet.emplace(as_nbarray(std::move(facet.value()))); + { + python_facet.emplace( + dolfinx_wrappers::as_nbarray(std::move(facet.value()))); + } return std::tuple{std::move(mesh1), std::move(python_cell), std::move(python_facet)}; }, - nb::arg("mesh"), nb::arg("edges") = nb::none(), nb::arg("redistribute"), - nb::arg("ghost_mode"), nb::arg("option")); + nb::arg("mesh"), nb::arg("edges").none(), nb::arg("partitioner").none(), + nb::arg("option")); } +} // namespace + +namespace dolfinx_wrappers +{ void refinement(nb::module_& m) { - export_refinement_with_variable_mesh_type(m); - export_refinement_with_variable_mesh_type(m); + export_refinement(m); + export_refinement(m); nb::enum_(m, "RefinementOption") .value("none", dolfinx::refinement::Option::none) diff --git a/python/test/unit/refinement/test_interval.py b/python/test/unit/refinement/test_interval.py index 28779644840..27040fb0b90 100644 --- a/python/test/unit/refinement/test_interval.py +++ b/python/test/unit/refinement/test_interval.py @@ -20,16 +20,14 @@ "ghost_mode_refined", [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet], ) -@pytest.mark.parametrize("redistribute", [True, False]) @pytest.mark.parametrize("option", [mesh.RefinementOption.none, mesh.RefinementOption.parent_cell]) -def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option): +def test_refine_interval(n, ghost_mode, ghost_mode_refined, option): msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode) msh_refined, edges, vertices = mesh.refine( msh, - redistribute=redistribute, - ghost_mode=ghost_mode_refined, option=option, ) + # TODO: add create_cell_partitioner(ghost_mode) when works # vertex count assert msh_refined.topology.index_map(0).size_global == 2 * n + 1 @@ -52,16 +50,14 @@ def test_refine_interval(n, ghost_mode, redistribute, ghost_mode_refined, option "ghost_mode_refined", [mesh.GhostMode.none, mesh.GhostMode.shared_vertex, mesh.GhostMode.shared_facet], ) -@pytest.mark.parametrize("redistribute", [True, False]) -def test_refine_interval_adaptive(n, ghost_mode, redistribute, ghost_mode_refined): +def test_refine_interval_adaptive(n, ghost_mode, ghost_mode_refined): msh = mesh.create_interval(MPI.COMM_WORLD, n, [0, 1], ghost_mode=ghost_mode) msh_refined, edges, vertices = mesh.refine( msh, np.arange(10, dtype=np.int32), - redistribute=redistribute, - ghost_mode=ghost_mode_refined, option=mesh.RefinementOption.parent_cell, ) + # TODO: add create_cell_partitioner(ghost_mode) when works # vertex count assert msh_refined.topology.index_map(0).size_global == n + 1 + 10 * MPI.COMM_WORLD.size diff --git a/python/test/unit/refinement/test_refinement.py b/python/test/unit/refinement/test_refinement.py index 6aaa4090dfb..0a8ee67b79a 100644 --- a/python/test/unit/refinement/test_refinement.py +++ b/python/test/unit/refinement/test_refinement.py @@ -11,6 +11,7 @@ from numpy import isclose import ufl +from dolfinx.cpp.mesh import create_cell_partitioner from dolfinx.fem import assemble_matrix, form, functionspace from dolfinx.mesh import ( CellType, @@ -32,7 +33,7 @@ def test_refine_create_unit_square(): """Refine mesh of unit square.""" mesh = create_unit_square(MPI.COMM_WORLD, 5, 7, ghost_mode=GhostMode.none) mesh.topology.create_entities(1) - mesh_refined, _, _ = refine(mesh, redistribute=False) + mesh_refined, _, _ = refine(mesh) assert mesh_refined.topology.index_map(0).size_global == 165 assert mesh_refined.topology.index_map(2).size_global == 280 @@ -46,7 +47,7 @@ def test_refine_create_unit_cube(ghost_mode, redistribute): """Refine mesh of unit cube.""" mesh = create_unit_cube(MPI.COMM_WORLD, 5, 7, 9, ghost_mode=ghost_mode) mesh.topology.create_entities(1) - mesh, _, _ = refine(mesh, redistribute=redistribute) + mesh, _, _ = refine(mesh, partitioner=create_cell_partitioner(ghost_mode)) assert mesh.topology.index_map(0).size_global == 3135 assert mesh.topology.index_map(3).size_global == 15120 @@ -58,8 +59,7 @@ def test_refine_create_form(): """Check that forms can be assembled on refined mesh""" mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3) mesh.topology.create_entities(1) - mesh, _, _ = refine(mesh, redistribute=True) - + mesh, _, _ = refine(mesh) V = functionspace(mesh, ("Lagrange", 1)) # Define variational problem @@ -71,20 +71,20 @@ def test_refine_create_form(): def test_sub_refine(): """Test that refinement of a subset of edges works""" - mesh = create_unit_square( + msh = create_unit_square( MPI.COMM_WORLD, 3, 4, diagonal=DiagonalType.left, ghost_mode=GhostMode.none ) - mesh.topology.create_entities(1) + msh.topology.create_entities(1) def left_corner_edge(x, tol=1e-7): return isclose(x[0], 0) & (x[1] < 1 / 4 + tol) - edges = locate_entities_boundary(mesh, 1, left_corner_edge) - if MPI.COMM_WORLD.size == 0: + edges = locate_entities_boundary(msh, 1, left_corner_edge) + if MPI.COMM_WORLD.size == 1: assert edges == 1 - mesh2, _, _ = refine(mesh, edges, redistribute=False) - assert mesh.topology.index_map(2).size_global + 3 == mesh2.topology.index_map(2).size_global + msh1, _, _ = refine(msh, edges) + assert msh.topology.index_map(2).size_global + 3 == msh1.topology.index_map(2).size_global def test_refine_from_cells(): @@ -100,12 +100,12 @@ def left_side(x, tol=1e-16): return x[0] <= 0.5 + tol cells = locate_entities(msh, msh.topology.dim, left_side) - if MPI.COMM_WORLD.size == 0: - assert cells.__len__() == Nx * Ny + if MPI.COMM_WORLD.size == 1: + assert len(cells) == Nx * Ny edges = compute_incident_entities(msh.topology, cells, 2, 1) - if MPI.COMM_WORLD.size == 0: - assert edges.__len__() == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny - mesh2, _, _ = refine(msh, edges, redistribute=True) + if MPI.COMM_WORLD.size == 1: + assert len(edges) == Nx // 2 * (2 * Ny + 1) + (Nx // 2 + 1) * Ny + mesh2, _, _ = refine(msh, edges) num_cells_global = mesh2.topology.index_map(2).size_global actual_cells = 3 * (Nx * Ny) + 3 * Ny + 2 * Nx * Ny @@ -116,61 +116,58 @@ def left_side(x, tol=1e-16): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda mesh: refine( - mesh, redistribute=False, option=RefinementOption.parent_cell_and_facet - ), - lambda mesh: refine( - mesh, - np.arange(mesh.topology.index_map(1).size_local), - redistribute=False, + lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet), + lambda msh: refine( + msh, + edges=np.arange(msh.topology.index_map(1).size_local), + partitioner=None, option=RefinementOption.parent_cell_and_facet, ), ], ) def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): if tdim == 3: - mesh = create_unit_cube( + msh = create_unit_cube( MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none ) else: - mesh = create_unit_square( - MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none - ) - mesh.topology.create_entities(tdim - 1) - mesh.topology.create_connectivity(tdim - 1, tdim) - mesh.topology.create_entities(1) - f_to_c = mesh.topology.connectivity(tdim - 1, tdim) + msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) + msh.topology.create_entities(tdim - 1) + msh.topology.create_connectivity(tdim - 1, tdim) + msh.topology.create_entities(1) + f_to_c = msh.topology.connectivity(tdim - 1, tdim) facet_indices = [] - for f in range(mesh.topology.index_map(tdim - 1).size_local): + for f in range(msh.topology.index_map(tdim - 1).size_local): if len(f_to_c.links(f)) == 1: facet_indices += [f] meshtag = meshtags( - mesh, + msh, tdim - 1, np.array(facet_indices, dtype=np.int32), np.arange(len(facet_indices), dtype=np.int32), ) - fine_mesh, parent_cell, parent_facet = refine_plaza_wrapper(mesh) - fine_mesh.topology.create_entities(tdim - 1) - new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell, parent_facet) + msh1, parent_cell, parent_facet = refine_plaza_wrapper(msh) + + msh1.topology.create_entities(tdim - 1) + new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) # New tags should be on facets with one cell (i.e. exterior) - fine_mesh.topology.create_connectivity(tdim - 1, tdim) - new_f_to_c = fine_mesh.topology.connectivity(tdim - 1, tdim) + msh1.topology.create_connectivity(tdim - 1, tdim) + new_f_to_c = msh1.topology.connectivity(tdim - 1, tdim) for f in new_meshtag.indices: assert len(new_f_to_c.links(f)) == 1 # Now mark all facets (including internal) - facet_indices = np.arange(mesh.topology.index_map(tdim - 1).size_local) + facet_indices = np.arange(msh.topology.index_map(tdim - 1).size_local) meshtag = meshtags( - mesh, + msh, tdim - 1, np.array(facet_indices, dtype=np.int32), np.arange(len(facet_indices), dtype=np.int32), ) - new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell, parent_facet) + new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell, parent_facet) assert len(new_meshtag.indices) == (tdim * 2 - 2) * len(meshtag.indices) @@ -178,44 +175,40 @@ def test_refine_facet_meshtag(tdim, refine_plaza_wrapper): @pytest.mark.parametrize( "refine_plaza_wrapper", [ - lambda mesh: refine( - mesh, redistribute=False, option=RefinementOption.parent_cell_and_facet - ), - lambda mesh: refine( - mesh, - np.arange(mesh.topology.index_map(1).size_local), - redistribute=False, + lambda msh: refine(msh, partitioner=None, option=RefinementOption.parent_cell_and_facet), + lambda msh: refine( + msh, + np.arange(msh.topology.index_map(1).size_local), + partitioner=None, option=RefinementOption.parent_cell_and_facet, ), ], ) def test_refine_cell_meshtag(tdim, refine_plaza_wrapper): if tdim == 3: - mesh = create_unit_cube( + msh = create_unit_cube( MPI.COMM_WORLD, 2, 3, 5, CellType.tetrahedron, ghost_mode=GhostMode.none ) else: - mesh = create_unit_square( - MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none - ) + msh = create_unit_square(MPI.COMM_WORLD, 2, 5, CellType.triangle, ghost_mode=GhostMode.none) - mesh.topology.create_entities(1) - cell_indices = np.arange(mesh.topology.index_map(tdim).size_local) + msh.topology.create_entities(1) + cell_indices = np.arange(msh.topology.index_map(tdim).size_local) meshtag = meshtags( - mesh, + msh, tdim, np.array(cell_indices, dtype=np.int32), np.arange(len(cell_indices), dtype=np.int32), ) - fine_mesh, parent_cell, _ = refine_plaza_wrapper(mesh) - new_meshtag = transfer_meshtag(meshtag, fine_mesh, parent_cell) + msh1, parent_cell, _ = refine_plaza_wrapper(msh) + new_meshtag = transfer_meshtag(meshtag, msh1, parent_cell) assert sum(new_meshtag.values) == (tdim * 4 - 4) * sum(meshtag.values) assert len(new_meshtag.indices) == (tdim * 4 - 4) * len(meshtag.indices) def test_refine_ufl_cargo(): - mesh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3) - mesh.topology.create_entities(1) - refined_mesh, _, _ = refine(mesh, redistribute=False) - assert refined_mesh.ufl_domain().ufl_cargo() != mesh.ufl_domain().ufl_cargo() + msh = create_unit_cube(MPI.COMM_WORLD, 4, 3, 3) + msh.topology.create_entities(1) + msh1, _, _ = refine(msh) + assert msh1.ufl_domain().ufl_cargo() != msh.ufl_domain().ufl_cargo()