diff --git a/cpp/dolfinx/common/IndexMap.cpp b/cpp/dolfinx/common/IndexMap.cpp index ba2eb8a669a..d7913eb7294 100644 --- a/cpp/dolfinx/common/IndexMap.cpp +++ b/cpp/dolfinx/common/IndexMap.cpp @@ -1004,7 +1004,8 @@ std::vector IndexMap::global_indices() const //----------------------------------------------------------------------------- MPI_Comm IndexMap::comm() const { return _comm.comm(); } //---------------------------------------------------------------------------- -graph::AdjacencyList IndexMap::index_to_dest_ranks(int tag) const +graph::AdjacencyList> +IndexMap::index_to_dest_ranks(int tag) const { const std::int64_t offset = _local_range[0]; diff --git a/cpp/dolfinx/common/IndexMap.h b/cpp/dolfinx/common/IndexMap.h index 60b5c216832..bdae1681bdc 100644 --- a/cpp/dolfinx/common/IndexMap.h +++ b/cpp/dolfinx/common/IndexMap.h @@ -232,7 +232,7 @@ class IndexMap /// std::int64_t>,std::span,int) for an explanation of when /// `tag` is required. /// @return Shared indices. - graph::AdjacencyList index_to_dest_ranks( + graph::AdjacencyList> index_to_dest_ranks( int tag = static_cast(dolfinx::MPI::tag::consensus_nbx)) const; /// @brief Build a list of owned indices that are ghosted by another diff --git a/cpp/dolfinx/common/utils.h b/cpp/dolfinx/common/utils.h index 8c65d00a957..00be1fee2ec 100644 --- a/cpp/dolfinx/common/utils.h +++ b/cpp/dolfinx/common/utils.h @@ -121,7 +121,9 @@ std::size_t hash_global(MPI_Comm comm, const T& x) /// data is data volume (`weight`) and local/remote memory indicator /// (`local==1` is an edge to an shared memory process/rank, other /// wise the target node is a remote memory rank). -std::string comm_to_json( - const graph::AdjacencyList, - std::pair>& g); +std::string +comm_to_json(const graph::AdjacencyList< + std::vector>, + std::vector, + std::vector>>& g); } // namespace dolfinx::common diff --git a/cpp/dolfinx/fem/DofMap.cpp b/cpp/dolfinx/fem/DofMap.cpp index d3ad810eca9..0c2e46daf4a 100644 --- a/cpp/dolfinx/fem/DofMap.cpp +++ b/cpp/dolfinx/fem/DofMap.cpp @@ -100,7 +100,7 @@ fem::DofMap build_collapsed_dofmap(const DofMap& dofmap_view, [&old_to_new](auto idx_old) { return old_to_new[idx_old]; }); // Dimension sanity checks - assert((int)dofmap.size() + assert(dofmap.size() == (cells->num_nodes() * dofmap_view.element_dof_layout().num_dofs())); assert(dofmap.size() % dofmap_view.element_dof_layout().num_dofs() == 0); @@ -115,7 +115,7 @@ fem::DofMap build_collapsed_dofmap(const DofMap& dofmap_view, } // namespace //----------------------------------------------------------------------------- -graph::AdjacencyList fem::transpose_dofmap( +graph::AdjacencyList> fem::transpose_dofmap( md::mdspan> dofmap, std::int32_t num_cells) { @@ -201,14 +201,15 @@ DofMap DofMap::extract_sub_dofmap(std::span component) const this->index_map_bs(), std::move(dofmap), 1); } //----------------------------------------------------------------------------- -std::pair> DofMap::collapse( - MPI_Comm comm, const mesh::Topology& topology, - std::function(const graph::AdjacencyList&)>&& - reorder_fn) const +std::pair> +DofMap::collapse(MPI_Comm comm, const mesh::Topology& topology, + std::function( + const graph::AdjacencyList>&)>&& + reorder_fn) const { if (!reorder_fn) { - reorder_fn = [](const graph::AdjacencyList& g) + reorder_fn = [](const graph::AdjacencyList>& g) { return graph::reorder_gps(g); }; } // Create new dofmap @@ -249,7 +250,7 @@ std::pair> DofMap::collapse( auto cells = topology.connectivity(tdim, 0); assert(cells); const int bs = dofmap_new.bs(); - for (std::int32_t c = 0; c < cells->num_nodes(); ++c) + for (std::size_t c = 0; c < cells->num_nodes(); ++c) { std::span cell_dofs_view = this->cell_dofs(c); std::span cell_dofs = dofmap_new.cell_dofs(c); diff --git a/cpp/dolfinx/fem/DofMap.h b/cpp/dolfinx/fem/DofMap.h index 85efae6b41f..3b6115a1534 100644 --- a/cpp/dolfinx/fem/DofMap.h +++ b/cpp/dolfinx/fem/DofMap.h @@ -59,7 +59,7 @@ namespace dolfinx::fem /// typically used to exclude ghost cell contributions. /// @return Map from global (process-wise) index to positions in an /// unaassembled array. The links for each node are sorted. -graph::AdjacencyList transpose_dofmap( +graph::AdjacencyList> transpose_dofmap( md::mdspan> dofmap, std::int32_t num_cells); @@ -143,11 +143,11 @@ class DofMap /// @param[in] reorder_fn Graph re-ordering function to apply to the /// dof data /// @return The collapsed dofmap - std::pair> - collapse(MPI_Comm comm, const mesh::Topology& topology, - std::function( - const graph::AdjacencyList&)>&& reorder_fn - = nullptr) const; + std::pair> collapse( + MPI_Comm comm, const mesh::Topology& topology, + std::function( + const graph::AdjacencyList>&)>&& reorder_fn + = nullptr) const; /// @brief Get dofmap data /// @return The adjacency list with dof indices for each cell diff --git a/cpp/dolfinx/fem/dofmapbuilder.cpp b/cpp/dolfinx/fem/dofmapbuilder.cpp index 7b8d8d340fd..229dc1d2006 100644 --- a/cpp/dolfinx/fem/dofmapbuilder.cpp +++ b/cpp/dolfinx/fem/dofmapbuilder.cpp @@ -48,11 +48,11 @@ struct dofmap_t /// @param[in] reorder_fn The graph reordering function to apply /// @return Map from original_to_contiguous[i] to new index after /// reordering -std::vector -reorder_owned(const std::vector& dofmaps, std::int32_t owned_size, - const std::vector& original_to_contiguous, - const std::function( - const graph::AdjacencyList&)>& reorder_fn) +std::vector reorder_owned( + const std::vector& dofmaps, std::int32_t owned_size, + const std::vector& original_to_contiguous, + const std::function( + const graph::AdjacencyList>&)>& reorder_fn) { std::vector graph_data, graph_offsets; @@ -404,7 +404,7 @@ std::pair, std::int32_t> compute_reordering_map( const std::vector>& dof_entity, const std::vector>& index_maps, const std::function( - const graph::AdjacencyList&)>& reorder_fn) + const graph::AdjacencyList>&)>& reorder_fn) { common::Timer t0("Dofmap builder: compute dof reordering map"); @@ -642,7 +642,7 @@ fem::build_dofmap_data( MPI_Comm comm, const mesh::Topology& topology, const std::vector& element_dof_layouts, const std::function( - const graph::AdjacencyList&)>& reorder_fn) + const graph::AdjacencyList>&)>& reorder_fn) { common::Timer t0("Dofmap builder: build dofmap data"); diff --git a/cpp/dolfinx/fem/dofmapbuilder.h b/cpp/dolfinx/fem/dofmapbuilder.h index 77a12162799..b8efb61e24b 100644 --- a/cpp/dolfinx/fem/dofmapbuilder.h +++ b/cpp/dolfinx/fem/dofmapbuilder.h @@ -36,9 +36,10 @@ class ElementDofLayout; /// the dofmaps /// @return The index map, block size, and dofmaps for each element type std::tuple>> -build_dofmap_data(MPI_Comm comm, const mesh::Topology& topology, - const std::vector& element_dof_layouts, - const std::function( - const graph::AdjacencyList&)>& reorder_fn); +build_dofmap_data( + MPI_Comm comm, const mesh::Topology& topology, + const std::vector& element_dof_layouts, + const std::function( + const graph::AdjacencyList>&)>& reorder_fn); } // namespace dolfinx::fem diff --git a/cpp/dolfinx/fem/utils.cpp b/cpp/dolfinx/fem/utils.cpp index 1df0cec6f6a..31f2b6326db 100644 --- a/cpp/dolfinx/fem/utils.cpp +++ b/cpp/dolfinx/fem/utils.cpp @@ -34,7 +34,7 @@ fem::DofMap fem::create_dofmap( const std::function, std::uint32_t)>& permute_inv, const std::function( - const graph::AdjacencyList&)>& reorder_fn) + const graph::AdjacencyList>&)>& reorder_fn) { // Create required mesh entities const int D = topology.dim(); @@ -81,7 +81,7 @@ std::vector fem::create_dofmaps( const std::function, std::uint32_t)>& permute_inv, const std::function( - const graph::AdjacencyList&)>& reorder_fn) + const graph::AdjacencyList>&)>& reorder_fn) { std::int32_t D = topology.dim(); assert(layouts.size() == topology.entity_types(D).size()); @@ -189,8 +189,9 @@ fem::compute_integration_domains(fem::IntegralType integral_type, } auto get_connectivities = [tdim, &topology](int entity_dim) - -> std::pair>, - std::shared_ptr>> + -> std::pair< + std::shared_ptr>>, + std::shared_ptr>>> { auto e_to_c = topology.connectivity(entity_dim, tdim); if (!e_to_c) diff --git a/cpp/dolfinx/fem/utils.h b/cpp/dolfinx/fem/utils.h index 95cad88af93..499d1963c56 100644 --- a/cpp/dolfinx/fem/utils.h +++ b/cpp/dolfinx/fem/utils.h @@ -70,9 +70,9 @@ namespace impl /// @param[in] c_to_f Cell to facet connectivity /// @return Vector of (cell, local_facet) pairs template -std::array -get_cell_facet_pairs(std::int32_t f, std::span cells, - const graph::AdjacencyList& c_to_f) +std::array get_cell_facet_pairs( + std::int32_t f, std::span cells, + const graph::AdjacencyList>& c_to_f) { // Loop over cells sharing facet assert(cells.size() == num_cells); @@ -100,9 +100,9 @@ get_cell_facet_pairs(std::int32_t f, std::span cells, /// @param[in] c_to_e Cell to entity connectivity /// @return Vector of (cell, local_entity) pairs template -std::array -get_cell_entity_pairs(std::int32_t e, std::span cells, - const graph::AdjacencyList& c_to_e) +std::array get_cell_entity_pairs( + std::int32_t e, std::span cells, + const graph::AdjacencyList>& c_to_e) { static_assert(num_cells == 1); // Patch assembly not supported. @@ -355,13 +355,12 @@ ElementDofLayout create_element_dof_layout(const fem::FiniteElement& element, /// when transformation is not required. /// @param[in] reorder_fn Graph reordering function called on the dofmap /// @return A new dof map -DofMap -create_dofmap(MPI_Comm comm, const ElementDofLayout& layout, - mesh::Topology& topology, - const std::function, std::uint32_t)>& - permute_inv, - const std::function( - const graph::AdjacencyList&)>& reorder_fn); +DofMap create_dofmap( + MPI_Comm comm, const ElementDofLayout& layout, mesh::Topology& topology, + const std::function, std::uint32_t)>& + permute_inv, + const std::function( + const graph::AdjacencyList>&)>& reorder_fn); /// @brief Create a set of dofmaps on a given topology /// @param[in] comm MPI communicator @@ -379,7 +378,7 @@ std::vector create_dofmaps( const std::function, std::uint32_t)>& permute_inv, const std::function( - const graph::AdjacencyList&)>& reorder_fn); + const graph::AdjacencyList>&)>& reorder_fn); /// Get the name of each coefficient in a UFC form /// @param[in] ufcx_form The UFC form @@ -924,7 +923,8 @@ template FunctionSpace create_functionspace( std::shared_ptr> mesh, std::shared_ptr> e, - std::function(const graph::AdjacencyList&)> + std::function( + const graph::AdjacencyList>&)> reorder_fn = nullptr) { diff --git a/cpp/dolfinx/geometry/utils.h b/cpp/dolfinx/geometry/utils.h index 27f82271f81..3680020f7e3 100644 --- a/cpp/dolfinx/geometry/utils.h +++ b/cpp/dolfinx/geometry/utils.h @@ -476,7 +476,7 @@ std::vector compute_collisions(const BoundingBoxTree& tree0, /// @return For each point, the bounding box leaves that collide with /// the point. template -graph::AdjacencyList +graph::AdjacencyList> compute_collisions(const BoundingBoxTree& tree, std::span points) { if (tree.num_bboxes() > 0) @@ -626,9 +626,9 @@ compute_closest_entity(const BoundingBoxTree& tree, /// 3)`). Storage is row-major. /// @return For each point, the cells that collide with the point. template -graph::AdjacencyList compute_colliding_cells( +graph::AdjacencyList> compute_colliding_cells( const mesh::Mesh& mesh, - const graph::AdjacencyList& candidate_cells, + const graph::AdjacencyList>& candidate_cells, std::span points) { std::vector offsets = {0}; @@ -636,7 +636,7 @@ graph::AdjacencyList compute_colliding_cells( std::vector colliding_cells; constexpr T eps2 = 1e-12; const int tdim = mesh.topology()->dim(); - for (std::int32_t i = 0; i < candidate_cells.num_nodes(); i++) + for (std::size_t i = 0; i < candidate_cells.num_nodes(); ++i) { auto cells = candidate_cells.links(i); std::vector _point(3 * cells.size()); @@ -707,7 +707,8 @@ determine_point_ownership(const mesh::Mesh& mesh, std::span points, graph::AdjacencyList collisions = compute_collisions(global_bbtree, points); // Get unique list of outgoing ranks - std::vector out_ranks = collisions.array(); + std::vector out_ranks(collisions.array().begin(), + collisions.array().end()); std::ranges::sort(out_ranks); auto [unique_end, range_end] = std::ranges::unique(out_ranks); out_ranks.erase(unique_end, range_end); @@ -782,7 +783,7 @@ determine_point_ownership(const mesh::Mesh& mesh, std::span points, auto x_dofmap = geometry.dofmap(); // Compute candidate cells for collisions (and extrapolation) - const graph::AdjacencyList candidate_collisions + const graph::AdjacencyList> candidate_collisions = compute_collisions(bb, std::span(received_points.data(), received_points.size())); diff --git a/cpp/dolfinx/graph/AdjacencyList.h b/cpp/dolfinx/graph/AdjacencyList.h index 9963fc4a085..343b95290f4 100644 --- a/cpp/dolfinx/graph/AdjacencyList.h +++ b/cpp/dolfinx/graph/AdjacencyList.h @@ -1,4 +1,4 @@ -// Copyright (C) 2019-2025 Garth N. Wells +// Copyright (C) 2019-2026 Garth N. Wells // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -10,9 +10,10 @@ #include #include #include -#include +#include #include #include +#include #include #include @@ -31,73 +32,88 @@ namespace dolfinx::graph /// /// Node data can also be stored. /// -/// @tparam LinkData_t Graph link (edge) type. -/// @tparam NodeData_t Data type for graph node data. -template +/// @tparam EdgeContainer Edge container type. +/// @tparam OffsetContainer Offset container type. +/// @tparam NodeDataContainer Node data container type. +template , + typename NodeDataContainer = std::nullptr_t> + requires std::ranges::contiguous_range + and std::ranges::contiguous_range + and std::integral class AdjacencyList { -public: - /// @brief Adjacency list link (edge) type - using link_type = LinkData; - /// @brief Adjacency list node data type - using node_data_type = NodeData; - - /// @brief Construct trivial adjacency list where each of the n nodes - /// is connected to itself. - /// @param[in] n Number of nodes. - explicit AdjacencyList(const std::int32_t n) : _array(n), _offsets(n + 1) +private: + template + X create_iota(std::size_t n) { - std::iota(_array.begin(), _array.end(), 0); - std::iota(_offsets.begin(), _offsets.end(), 0); + X x(n); + std::iota(x.begin(), x.end(), 0); + return x; } +public: + /// @brief Adjacency list link (edge) type + using link_type = typename EdgeContainer::value_type; + /// @brief Adjacency list offset index type + using offset_type = typename OffsetContainer::value_type; + /// @brief Construct adjacency list from arrays of link (edge) data /// and offsets. - /// @param[in] data Adjacency lost data array. + /// + /// @param[in] data Adjacency list data array. /// @param[in] offsets Offsets into `data` for each node, where - /// `offsets[i]` is the first index in `data` for node `i`. The last - /// index in `offsets` is the equal to the length of `data`. array for - /// node `i`. - template - requires std::is_convertible_v, - std::vector> - and std::is_convertible_v, - std::vector> - AdjacencyList(U&& data, V&& offsets) - : _array(std::forward(data)), _offsets(std::forward(offsets)) + /// `offsets[i]` is the first index in `data` for node `i` and + /// `offsets[i+1] - offsets[i]` is the number of edges for node `i`. + /// The length of `offsets` is equal to the number of nodes plus 1. + template + requires std::is_convertible_v, EdgeContainer> + and std::is_convertible_v, + OffsetContainer> + AdjacencyList(W0&& data, W1&& offsets) + : _array(std::forward(data)), _offsets(std::forward(offsets)) { - _array.reserve(_offsets.back()); - assert(_offsets.back() == (std::int32_t)_array.size()); + assert(_offsets.back() <= (std::int32_t)_array.size()); } /// @brief Construct adjacency list from arrays of link (edge) data, /// offsets, and node data. - /// @param[in] data Adjacency lost data array. + /// + /// @param[in] data Adjacency list data array. /// @param[in] offsets Offsets into `data` for each node, where - /// `offsets[i]` is the first index in `data` for node `i`. The last - /// index in `offsets` is the equal to the length of `data`. + /// `offsets[i]` is the first index in `data` for node `i` and + /// `offsets[i+1] - offsets[i]` is the number of edges for node `i`. + /// The length of `offsets` is equal to the number of nodes plus 1. /// @param[in] node_data Node data array where `node_data[i]` is the /// data attached to node `i`. - template - requires std::is_convertible_v, - std::vector> - and std::is_convertible_v, - std::vector> - and std::is_convertible_v, - std::vector> - AdjacencyList(U&& data, V&& offsets, W&& node_data) - : _array(std::forward(data)), _offsets(std::forward(offsets)), - _node_data(std::forward(node_data)) + template + requires std::is_convertible_v, EdgeContainer> + and std::is_convertible_v, + OffsetContainer> + and std::is_convertible_v, + NodeDataContainer> + AdjacencyList(W0&& data, W1&& offsets, W2&& node_data) + : _array(std::forward(data)), _offsets(std::forward(offsets)), + _node_data(std::forward(node_data)) + { + assert(_node_data.size() == _offsets.size() - 1); + assert(_offsets.back() <= (std::int32_t)_array.size()); + } + + /// @brief Construct trivial adjacency list where each of the `n` + /// nodes is connected to itself. + /// + /// @param[in] n Number of nodes. + explicit AdjacencyList(std::int32_t n) + : AdjacencyList(create_iota>(n), + create_iota>(n + 1)) { - assert(_node_data.has_value() - and _node_data->size() == _offsets.size() - 1); - _array.reserve(_offsets.back()); - assert(_offsets.back() == (std::int32_t)_array.size()); } /// Set all connections for all entities (T is a '2D' container, e.g. /// a `std::vector<>`, /// `std::vector<>`, etc). + /// /// @param[in] data Adjacency list data, where `std::next(data, i)` /// points to the container of links (edges) for node `i`. template @@ -114,6 +130,30 @@ class AdjacencyList _array.insert(_array.end(), e.begin(), e.end()); } + /// @brief Construct an adjacency list with different template + /// parameters. + /// @param x Adjacency list to create the new adjacency list from. + template + AdjacencyList(const AdjacencyList& x) + : _array(x.array().begin(), x.array().end()), + _offsets(x.offsets().begin(), x.offsets().end()) + { + assert(_offsets.back() <= (std::int32_t)_array.size()); + } + + /// @brief Construct an adjacency list with different template + /// parameters. + /// @param x Adjacency list to create the new adjacency list from. + template + AdjacencyList(const AdjacencyList& x) + : _array(x.array().begin(), x.array().end()), + _offsets(x.offsets().begin(), x.offsets().end()), + _node_data(x.node_data().begin(), x.node_data().end()) + { + assert(_offsets.back() <= (std::int32_t)_array.size()); + assert(_node_data.size() == _offsets.size() - 1); + } + /// Copy constructor AdjacencyList(const AdjacencyList& list) = default; @@ -129,70 +169,132 @@ class AdjacencyList /// Move assignment operator AdjacencyList& operator=(AdjacencyList&& list) = default; - /// Equality operator - /// @return True is the adjacency lists are equal + /// @brief Equality operator + /// + /// @return True is the adjacency lists are equal. bool operator==(const AdjacencyList& list) const { - return this->_array == list._array and this->_offsets == list._offsets; + return this->_array == list._array and this->_offsets == list._offsets + and this->_node_data == list._node_data; } /// @brief Get the number of nodes. + /// /// @return The number of nodes in the adjacency list - std::int32_t num_nodes() const { return _offsets.size() - 1; } + std::size_t num_nodes() const { return _offsets.size() - 1; } /// @brief Number of connections for given node. - /// @param[in] node Node index. + /// @param[in] n Node index. /// @return The number of outgoing links (edges) from the node. - int num_links(std::size_t node) const + int num_links(std::size_t n) const { - assert((node + 1) < _offsets.size()); - return _offsets[node + 1] - _offsets[node]; + assert((n + 1) < _offsets.size()); + return *std::next(_offsets.begin(), n + 1) + - *std::next(_offsets.begin(), n); } /// @brief Get the links (edges) for given node. + /// /// @param[in] node Node index. /// @return Array of outgoing links for the node. The length will be /// `AdjacencyList::num_links(node)`. - std::span links(std::size_t node) + auto links(std::size_t node) { - return std::span(_array.data() + _offsets[node], - _offsets[node + 1] - _offsets[node]); + return std::span(_array.data() + *std::next(_offsets.begin(), node), + this->num_links(node)); } /// @brief Get the links (edges) for given node (const version). - /// @param[in] node Node index. + /// + /// @param[in] n Node index. /// @return Array of outgoing links for the node. The length will be /// `AdjacencyList:num_links(node)`. - std::span links(std::size_t node) const + std::span links(std::size_t n) const { - return std::span(_array.data() + _offsets[node], - _offsets[node + 1] - _offsets[node]); + return std::span(_array.data() + *std::next(_offsets.begin(), n), + this->num_links(n)); } - /// Return contiguous array of links for all nodes (const version). - const std::vector& array() const { return _array; } + /// @brief Create a sub-view of the adjacency list for a contiguous + /// range of nodes. + /// + /// The range of nodes is [n0, n1), i.e. it includes node n0 but not n1. + /// + /// @param n0 Start node index (inclusive). + /// @param n1 End node index (exclusive). + /// @return Sub-view of the adjacency list for nodes in the range [n0, n1). + auto sub(std::size_t n0, std::size_t n1) const + requires requires(NodeDataContainer) { + typename NodeDataContainer::value_type; + } + { + assert(n0 <= n1); + assert(n1 <= this->num_nodes()); + return AdjacencyList, + std::span, std::nullptr_t>( + std::span(_array), std::span(_offsets.data() + n0, n1 - n0 + 1), + std::span(_node_data.data() + n0, n1 - n0)); + } + + /// @brief Create a sub-view of the adjacency list for a contiguous + /// range of nodes. + /// + /// The range of nodes is [n0, n1), i.e. it includes node n0 but not n1. + /// + /// @param n0 Start node index (inclusive). + /// @param n1 End node index (exclusive). + /// @return Sub-view of the adjacency list for nodes in the range [n0, n1). + AdjacencyList, std::span, + std::nullptr_t> + sub(std::size_t n0, std::size_t n1) const + { + assert(n0 <= n1); + assert(n1 <= this->num_nodes()); + return AdjacencyList, + std::span, std::nullptr_t>( + std::span(_array), std::span(_offsets.data() + n0, n1 - n0 + 1)); + } - /// Return contiguous array of links for all nodes. - std::vector& array() { return _array; } + /// @brief Return container of links for all nodes (const version). + /// + /// @note The returned container may lager than the data pointed to by + /// AdjacencyList::offsets. + const EdgeContainer& array() const { return _array; } - /// Offset for each node in array() (const version). - const std::vector& offsets() const { return _offsets; } + /// @brief Return container of links for all nodes. + /// + /// @note The returned container may lager than the data pointed to by + /// AdjacencyList::offsets. + EdgeContainer& array() { return _array; } - /// Offset for each node in array(). - std::vector& offsets() { return _offsets; } + /// @brief Offset for each node in array() (const version). + const OffsetContainer& offsets() const { return _offsets; } + + /// @brief Offset for each node in array(). + OffsetContainer& offsets() { return _offsets; } /// Return node data (if present), where `node_data()[i]` is the data /// for node `i` (const version). - const std::optional>& node_data() const + const NodeDataContainer& node_data() const + requires requires(NodeDataContainer) { + typename NodeDataContainer::value_type; + } { return _node_data; } - /// Return node data (if present), where `node_data()[i]` is the data for node - /// `i`. - std::optional>& node_data() { return _node_data; } + /// Return node data (if present), where `node_data()[i]` is the data + /// for node `i`. + NodeDataContainer& node_data() + requires requires(NodeDataContainer) { + typename NodeDataContainer::value_type; + } + { + return _node_data; + } /// @brief Informal string representation (pretty-print). + /// /// @return String representation of the adjacency list. std::string str() const { @@ -212,23 +314,14 @@ class AdjacencyList private: // Connections (links/edges) for all entities stored as a contiguous // array - std::vector _array; + EdgeContainer _array; // Position of first connection for each entity (using local index) - std::vector _offsets; + OffsetContainer _offsets; // Node data, where _node_data[i] is the data associated with node `i` - std::optional> _node_data = std::nullopt; -}; - -/// @private Deduction -template -AdjacencyList(T, U) -> AdjacencyList; - -/// @private Deduction -template -AdjacencyList(T, U, W) - -> AdjacencyList; + NodeDataContainer _node_data; +}; // namespace dolfinx::graph /// @brief Construct a constant degree (valency) adjacency list. /// @@ -238,14 +331,12 @@ AdjacencyList(T, U, W) /// @param[in] data Adjacency array. /// @param[in] degree Number of (outgoing) links for each node. /// @return An adjacency list. -template - requires requires { - typename std::decay_t::value_type; - requires std::convertible_to< - U, std::vector::value_type>>; - } -AdjacencyList::value_type, V> -regular_adjacency_list(U&& data, int degree) +template > + requires std::ranges::contiguous_range + and std::ranges::contiguous_range +AdjacencyList +regular_adjacency_list(EdgeContainer&& data, int degree) { if (degree == 0 and !data.empty()) { @@ -259,12 +350,24 @@ regular_adjacency_list(U&& data, int degree) "Incompatible data size and degree for constant degree AdjacencyList"); } - std::int32_t num_nodes = degree == 0 ? data.size() : data.size() / degree; - std::vector offsets(num_nodes + 1, 0); - for (std::size_t i = 1; i < offsets.size(); ++i) - offsets[i] = offsets[i - 1] + degree; - return AdjacencyList::value_type, V>( - std::forward(data), std::move(offsets)); + std::size_t num_nodes = degree == 0 ? data.size() : data.size() / degree; + using int_type = typename std::decay_t::value_type; + auto offsets + = std::views::iota(int_type(0), int_type(num_nodes + 1)) + | std::views::transform([degree](auto i) { return degree * i; }); + return AdjacencyList(std::forward(data), + OffsetContainer(offsets.begin(), offsets.end())); } +/// @private Deduction +template +AdjacencyList(EdgeContainer, OffsetContainer) + -> AdjacencyList; + +/// @private Deduction +template +AdjacencyList(EdgeContainer, OffsetContainer, NodeDataContainer) + -> AdjacencyList; + } // namespace dolfinx::graph diff --git a/cpp/dolfinx/graph/ordering.cpp b/cpp/dolfinx/graph/ordering.cpp index 0c23ba8df35..f9dd4d5ca45 100644 --- a/cpp/dolfinx/graph/ordering.cpp +++ b/cpp/dolfinx/graph/ordering.cpp @@ -22,7 +22,7 @@ namespace // Compute the sets of connected components of the input "graph" which // contain the nodes in "indices". std::vector> -residual_graph_components(const graph::AdjacencyList& graph, +residual_graph_components(const graph::AdjacencyList>& graph, std::span indices) { if (indices.empty()) @@ -75,17 +75,18 @@ residual_graph_components(const graph::AdjacencyList& graph, } //----------------------------------------------------------------------------- // Get the (maximum) width of a level structure -int max_level_width(const graph::AdjacencyList& levels) +int max_level_width(const graph::AdjacencyList>& levels) { int wmax = 0; - for (int i = 0; i < levels.num_nodes(); ++i) + for (std::size_t i = 0; i < levels.num_nodes(); ++i) wmax = std::max(wmax, levels.num_links(i)); return wmax; } //----------------------------------------------------------------------------- // Create a level structure from graph, rooted at node s -graph::AdjacencyList -create_level_structure(const graph::AdjacencyList& graph, int s) +graph::AdjacencyList> +create_level_structure(const graph::AdjacencyList>& graph, + int s) { common::Timer t("GPS: create_level_structure"); @@ -125,9 +126,9 @@ create_level_structure(const graph::AdjacencyList& graph, int s) // Gibbs-Poole-Stockmeyer algorithm, finding a reordering for the given // graph, operating only on nodes which are yet unlabelled (indicated // with -1 in the vector rlabel). -std::vector -gps_reorder_unlabelled(const graph::AdjacencyList& graph, - std::span rlabel) +std::vector gps_reorder_unlabelled( + const graph::AdjacencyList>& graph, + std::span rlabel) { common::Timer timer("Gibbs-Poole-Stockmeyer ordering"); @@ -152,8 +153,8 @@ gps_reorder_unlabelled(const graph::AdjacencyList& graph, } // B. Generate a level structure Lv rooted at vertex v. - graph::AdjacencyList lv = create_level_structure(graph, v); - graph::AdjacencyList lu(0); + graph::AdjacencyList> lv = create_level_structure(graph, v); + graph::AdjacencyList> lu(0); bool done = false; int u = 0; std::vector S; @@ -172,7 +173,8 @@ gps_reorder_unlabelled(const graph::AdjacencyList& graph, // in order of increasing degree. for (int s : S) { - graph::AdjacencyList lstmp = create_level_structure(graph, s); + graph::AdjacencyList> lstmp + = create_level_structure(graph, s); if (lstmp.num_nodes() > lv.num_nodes()) { // Found a deeper level structure, so restart @@ -358,7 +360,7 @@ gps_reorder_unlabelled(const graph::AdjacencyList& graph, //----------------------------------------------------------------------------- std::vector -graph::reorder_gps(const graph::AdjacencyList& graph) +graph::reorder_gps(const graph::AdjacencyList>& graph) { const std::int32_t n = graph.num_nodes(); std::vector r(n, -1); diff --git a/cpp/dolfinx/graph/ordering.h b/cpp/dolfinx/graph/ordering.h index 42d11d66275..5d9e43dbb9a 100644 --- a/cpp/dolfinx/graph/ordering.h +++ b/cpp/dolfinx/graph/ordering.h @@ -22,6 +22,6 @@ namespace dolfinx::graph /// @return Reordering array `map`, where `map[i]` is the new index of /// node `i` std::vector -reorder_gps(const graph::AdjacencyList& graph); +reorder_gps(const graph::AdjacencyList>& graph); } // namespace dolfinx::graph diff --git a/cpp/dolfinx/graph/partition.cpp b/cpp/dolfinx/graph/partition.cpp index e77415d51ac..a1ebf20b173 100644 --- a/cpp/dolfinx/graph/partition.cpp +++ b/cpp/dolfinx/graph/partition.cpp @@ -17,10 +17,9 @@ using namespace dolfinx; //----------------------------------------------------------------------------- -graph::AdjacencyList -graph::partition_graph(MPI_Comm comm, int nparts, - const AdjacencyList& local_graph, - bool ghosting) +graph::AdjacencyList> graph::partition_graph( + MPI_Comm comm, int nparts, + const AdjacencyList>& local_graph, bool ghosting) { #if HAS_PARMETIS return graph::parmetis::partitioner()(comm, nparts, local_graph, ghosting); @@ -33,15 +32,15 @@ graph::partition_graph(MPI_Comm comm, int nparts, #endif } //----------------------------------------------------------------------------- -std::tuple, std::vector, +std::tuple>, std::vector, std::vector, std::vector> -graph::build::distribute(MPI_Comm comm, - const graph::AdjacencyList& list, - const graph::AdjacencyList& destinations) +graph::build::distribute( + MPI_Comm comm, const graph::AdjacencyList>& list, + const graph::AdjacencyList>& destinations) { common::Timer timer("Distribute AdjacencyList nodes to destination ranks"); - assert(list.num_nodes() == (int)destinations.num_nodes()); + assert(list.num_nodes() == destinations.num_nodes()); const int rank = dolfinx::MPI::rank(comm); // Get global offset for converting local index to global index for @@ -67,12 +66,12 @@ graph::build::distribute(MPI_Comm comm, // Build (dest, index, owning rank) list and sort std::vector> dest_to_index; dest_to_index.reserve(destinations.array().size()); - for (std::int32_t i = 0; i < destinations.num_nodes(); ++i) + for (std::size_t i = 0; i < destinations.num_nodes(); ++i) { auto di = destinations.links(i); std::ranges::transform(di, std::back_inserter(dest_to_index), [i, d0 = di.front()](auto d) -> std::array - { return {d, i, d0}; }); + { return {d, (std::int32_t)i, d0}; }); } std::ranges::sort(dest_to_index); @@ -219,23 +218,24 @@ graph::build::distribute(MPI_Comm comm, global_indices.shrink_to_fit(); ghost_index_owner.shrink_to_fit(); - return { - graph::AdjacencyList(std::move(data), std::move(offsets)), - std::move(src_ranks), std::move(global_indices), - std::move(ghost_index_owner)}; + return {graph::AdjacencyList>(std::move(data), + std::move(offsets)), + std::move(src_ranks), std::move(global_indices), + std::move(ghost_index_owner)}; } //----------------------------------------------------------------------------- std::tuple, std::vector, std::vector, std::vector> -graph::build::distribute(MPI_Comm comm, std::span list, - std::array shape, - const graph::AdjacencyList& destinations) +graph::build::distribute( + MPI_Comm comm, std::span list, + std::array shape, + const graph::AdjacencyList>& destinations) { common::Timer timer( "Distribute fixed-degree adjacency list to destination ranks"); assert(list.size() == shape[0] * shape[1]); - assert(destinations.num_nodes() == (std::int32_t)shape[0]); + assert(destinations.num_nodes() == shape[0]); int rank = dolfinx::MPI::rank(comm); std::int64_t num_owned = destinations.num_nodes(); @@ -251,12 +251,12 @@ graph::build::distribute(MPI_Comm comm, std::span list, // Build (dest, index, owning rank) list and sort std::vector> dest_to_index; dest_to_index.reserve(destinations.array().size()); - for (std::int32_t i = 0; i < destinations.num_nodes(); ++i) + for (std::size_t i = 0; i < destinations.num_nodes(); ++i) { auto di = destinations.links(i); std::ranges::transform(di, std::back_inserter(dest_to_index), [i, d0 = di.front()](auto d) -> std::array - { return {d, i, d0}; }); + { return {d, (std::int32_t)i, d0}; }); } std::ranges::sort(dest_to_index); diff --git a/cpp/dolfinx/graph/partition.h b/cpp/dolfinx/graph/partition.h index 95971a93538..3bf14de3549 100644 --- a/cpp/dolfinx/graph/partition.h +++ b/cpp/dolfinx/graph/partition.h @@ -28,8 +28,9 @@ namespace dolfinx::graph /// @param[in] ghosting Flag to enable ghosting of the output node /// distribution /// @return Destination rank for each input node -using partition_fn = std::function( - MPI_Comm, int, const AdjacencyList&, bool)>; +using partition_fn + = std::function>( + MPI_Comm, int, const AdjacencyList>&, bool)>; /// @brief Partition graph across processes using the default graph /// partitioner. @@ -41,9 +42,10 @@ using partition_fn = std::function( /// @param[in] ghosting Flag to enable ghosting of the output node /// distribution. /// @return Destination rank for each input node. -AdjacencyList +AdjacencyList> partition_graph(MPI_Comm comm, int nparts, - const AdjacencyList& local_graph, bool ghosting); + const AdjacencyList>& local_graph, + bool ghosting); /// Tools for distributed graphs /// @@ -64,10 +66,11 @@ namespace build /// 2. Source ranks for each node in the adjacency list /// 3. Original global index for each node in the adjacency list /// 4. Owning rank of ghost nodes. -std::tuple, std::vector, +std::tuple>, std::vector, std::vector, std::vector> -distribute(MPI_Comm comm, const graph::AdjacencyList& list, - const graph::AdjacencyList& destinations); +distribute(MPI_Comm comm, + const graph::AdjacencyList>& list, + const graph::AdjacencyList>& destinations); /// @brief Distribute fixed size nodes to destination ranks. /// @@ -90,7 +93,7 @@ std::tuple, std::vector, std::vector, std::vector> distribute(MPI_Comm comm, std::span list, std::array shape, - const graph::AdjacencyList& destinations); + const graph::AdjacencyList>& destinations); /// @brief Take a set of distributed input global indices, including /// ghosts, and determine the new global indices after remapping. diff --git a/cpp/dolfinx/graph/partitioners.cpp b/cpp/dolfinx/graph/partitioners.cpp index 05ae0fe8fd3..1cf94169a79 100644 --- a/cpp/dolfinx/graph/partitioners.cpp +++ b/cpp/dolfinx/graph/partitioners.cpp @@ -36,8 +36,9 @@ extern "C" using namespace dolfinx; template -graph::AdjacencyList dolfinx::graph::compute_destination_ranks( - MPI_Comm comm, const graph::AdjacencyList& graph, +graph::AdjacencyList> +dolfinx::graph::compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList>& graph, const std::vector& node_disp, const std::vector& part) { common::Timer timer("Extend graph destination ranks for halo"); @@ -45,13 +46,13 @@ graph::AdjacencyList dolfinx::graph::compute_destination_ranks( const int rank = dolfinx::MPI::rank(comm); const std::int64_t range0 = node_disp[rank]; const std::int64_t range1 = node_disp[rank + 1]; - assert(static_cast(range1 - range0) == graph.num_nodes()); + assert((range1 - range0) == (std::int64_t)graph.num_nodes()); // Wherever an owned 'node' goes, so must the nodes connected to it by // an edge ('node1'). Task is to let the owner of node1 know the extra // ranks that it needs to send node1 to. std::vector> node_to_dest; - for (int node0 = 0; node0 < graph.num_nodes(); ++node0) + for (std::size_t node0 = 0; node0 < graph.num_nodes(); ++node0) { // Wherever 'node' goes to, so must the attached 'node1' for (auto node1 : graph.links(node0)) @@ -177,10 +178,10 @@ graph::AdjacencyList dolfinx::graph::compute_destination_ranks( data[pos[x0]++] = x1; } - graph::AdjacencyList g(std::move(data), std::move(offsets)); + graph::AdjacencyList> g(std::move(data), std::move(offsets)); // Make sure the owning rank comes first for each node - for (std::int32_t i = 0; i < g.num_nodes(); ++i) + for (std::size_t i = 0; i < g.num_nodes(); ++i) { auto d = g.links(i); auto it = std::find(d.begin(), d.end(), part[i]); @@ -192,12 +193,14 @@ graph::AdjacencyList dolfinx::graph::compute_destination_ranks( } /// @cond -template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( - MPI_Comm comm, const graph::AdjacencyList& graph, +template graph::AdjacencyList> +dolfinx::graph::compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList>& graph, const std::vector& node_disp, const std::vector& part); -template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( - MPI_Comm comm, const graph::AdjacencyList& graph, +template graph::AdjacencyList> +dolfinx::graph::compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList>& graph, const std::vector& node_disp, const std::vector& part); /// @endcond @@ -205,9 +208,10 @@ template graph::AdjacencyList dolfinx::graph::compute_destination_ranks( //----------------------------------------------------------------------------- #ifdef HAS_PARMETIS template -std::vector adaptive_repartition(MPI_Comm comm, - const graph::AdjacencyList& adj_graph, - double weight) +std::vector +adaptive_repartition(MPI_Comm comm, + const graph::AdjacencyList>& adj_graph, + double weight) { common::Timer timer( "Compute graph partition (ParMETIS Adaptive Repartition)"); @@ -254,7 +258,8 @@ std::vector adaptive_repartition(MPI_Comm comm, } //----------------------------------------------------------------------------- template -std::vector refine(MPI_Comm comm, const graph::AdjacencyList& adj_graph) +std::vector refine(MPI_Comm comm, + const graph::AdjacencyList>& adj_graph) { common::Timer timer("Compute graph partition (ParMETIS Refine)"); @@ -310,9 +315,10 @@ std::vector refine(MPI_Comm comm, const graph::AdjacencyList& adj_graph) graph::partition_fn graph::scotch::partitioner(graph::scotch::strategy strategy, double imbalance, int seed) { - return [imbalance, strategy, seed](MPI_Comm comm, int nparts, - const AdjacencyList& graph, - bool ghosting) + return [imbalance, strategy, + seed](MPI_Comm comm, int nparts, + const AdjacencyList>& graph, + bool ghosting) { spdlog::info("Compute graph partition using PT-SCOTCH"); common::Timer timer("Compute graph partition (SCOTCH)"); @@ -463,7 +469,7 @@ graph::partition_fn graph::scotch::partitioner(graph::scotch::strategy strategy, // Create a map of local nodes to their additional destination // processes, due to ghosting std::map> local_node_to_dests; - for (std::int32_t node0 = 0; node0 < graph.num_nodes(); ++node0) + for (std::size_t node0 = 0; node0 < graph.num_nodes(); ++node0) { // Get all edges outward from node i const std::int32_t node0_rank = node_partition[node0]; @@ -479,7 +485,7 @@ graph::partition_fn graph::scotch::partitioner(graph::scotch::strategy strategy, timer5.stop(); offsets.reserve(graph.num_nodes() + 1); - for (std::int32_t i = 0; i < graph.num_nodes(); ++i) + for (std::size_t i = 0; i < graph.num_nodes(); ++i) { dests.push_back(node_partition[i]); if (auto it = local_node_to_dests.find(i); @@ -515,9 +521,10 @@ graph::partition_fn graph::scotch::partitioner(graph::scotch::strategy strategy, graph::partition_fn graph::parmetis::partitioner(double imbalance, std::array options) { - return [imbalance, options](MPI_Comm comm, idx_t nparts, - const graph::AdjacencyList& graph, - bool ghosting) + return [imbalance, + options](MPI_Comm comm, idx_t nparts, + const graph::AdjacencyList>& graph, + bool ghosting) { spdlog::info("Compute graph partition using ParMETIS"); common::Timer timer("Compute graph partition (ParMETIS)"); @@ -585,7 +592,7 @@ graph::partition_fn graph::parmetis::partitioner(double imbalance, if (ghosting and pcomm != MPI_COMM_NULL) { // FIXME: Is it implicit that the first entry is the owner? - graph::AdjacencyList dest + graph::AdjacencyList> dest = graph::compute_destination_ranks(pcomm, graph, node_disp, part); if (split_comm) MPI_Comm_free(&pcomm); @@ -612,7 +619,8 @@ graph::partition_fn graph::kahip::partitioner(int mode, int seed, { return [mode, seed, imbalance, suppress_output]( MPI_Comm comm, int nparts, - const graph::AdjacencyList& graph, bool ghosting) + const graph::AdjacencyList>& graph, + bool ghosting) { spdlog::info("Compute graph partition using (parallel) KaHIP"); diff --git a/cpp/dolfinx/graph/partitioners.h b/cpp/dolfinx/graph/partitioners.h index 444dc6d8ed6..2136e55e75c 100644 --- a/cpp/dolfinx/graph/partitioners.h +++ b/cpp/dolfinx/graph/partitioners.h @@ -24,8 +24,8 @@ namespace dolfinx::graph /// is the destination of the node with local index `i`. /// @return Destination ranks for each local node. template -graph::AdjacencyList compute_destination_ranks( - MPI_Comm comm, const graph::AdjacencyList& graph, +graph::AdjacencyList> compute_destination_ranks( + MPI_Comm comm, const graph::AdjacencyList>& graph, const std::vector& node_disp, const std::vector& part); namespace scotch diff --git a/cpp/dolfinx/graph/utils.cpp b/cpp/dolfinx/graph/utils.cpp index f015467a686..236117166ee 100644 --- a/cpp/dolfinx/graph/utils.cpp +++ b/cpp/dolfinx/graph/utils.cpp @@ -14,8 +14,9 @@ using namespace dolfinx; //----------------------------------------------------------------------------- -graph::AdjacencyList, - std::pair> +graph::AdjacencyList>, + std::vector, + std::vector>> graph::comm_graph(const common::IndexMap& map, int root) { MPI_Comm comm = map.comm(); @@ -101,17 +102,19 @@ graph::comm_graph(const common::IndexMap& map, int root) std::move(sizes_remote)); } //----------------------------------------------------------------------------- -std::string graph::comm_to_json( - const graph::AdjacencyList, - std::pair>& g) +std::string +graph::comm_to_json(const graph::AdjacencyList< + std::vector>, + std::vector, + std::vector>>& g) { const std::vector>& node_weights - = g.node_data().value(); + = g.node_data(); std::stringstream out; out << std::format("{{\"directed\": true, \"multigraph\": false, \"graph\": " "[], \"nodes\": ["); - for (std::int32_t n = 0; n < g.num_nodes(); ++n) + for (std::size_t n = 0; n < g.num_nodes(); ++n) { // Note: it is helpful to order map keys alphabetically out << std::format("{{\"num_ghosts\": {}, \"weight\": {}, \"id\": {}}}", @@ -121,7 +124,7 @@ std::string graph::comm_to_json( } out << "], "; out << "\"adjacency\": ["; - for (std::int32_t n = 0; n < g.num_nodes(); ++n) + for (std::size_t n = 0; n < g.num_nodes(); ++n) { out << "["; auto links = g.links(n); diff --git a/cpp/dolfinx/graph/utils.h b/cpp/dolfinx/graph/utils.h index ba42525935e..e441748cd77 100644 --- a/cpp/dolfinx/graph/utils.h +++ b/cpp/dolfinx/graph/utils.h @@ -44,8 +44,9 @@ namespace dolfinx::graph /// local/remote is memory indicator (`local==1` is an edge to a /// shared memory node). Node data is (number of owned indices, number /// of ghost indices). -AdjacencyList, - std::pair> +AdjacencyList>, + std::vector, + std::vector>> comm_graph(const common::IndexMap& map, int root = 0); /// @brief Build communication graph data as a JSON string. @@ -61,7 +62,8 @@ comm_graph(const common::IndexMap& map, int root = 0); /// data is data volume (`weight`) and local/remote memory indicator /// (`local==1` is an edge to an shared memory process/rank, other /// wise the target node is a remote memory rank). -std::string -comm_to_json(const AdjacencyList, - std::pair>& g); +std::string comm_to_json( + const AdjacencyList>, + std::vector, + std::vector>>& g); } // namespace dolfinx::graph diff --git a/cpp/dolfinx/io/XDMFFile.cpp b/cpp/dolfinx/io/XDMFFile.cpp index b4c0513e41b..e1bebfc73b9 100644 --- a/cpp/dolfinx/io/XDMFFile.cpp +++ b/cpp/dolfinx/io/XDMFFile.cpp @@ -398,7 +398,7 @@ XDMFFile::read_meshtags(const mesh::Mesh& mesh, const std::string& name, mesh::cell_entity_type(mesh.topology()->cell_type(), mesh::cell_dim(cell_type), 0), 0); - const graph::AdjacencyList entities_adj + const graph::AdjacencyList> entities_adj = graph::regular_adjacency_list(std::move(entities_values.first), num_vertices_per_entity); mesh::MeshTags meshtags = mesh::create_meshtags( diff --git a/cpp/dolfinx/io/utils.h b/cpp/dolfinx/io/utils.h index 51e3a386a86..f39d4e50850 100644 --- a/cpp/dolfinx/io/utils.h +++ b/cpp/dolfinx/io/utils.h @@ -435,7 +435,7 @@ std::pair, std::vector> distribute_entity_data( throw std::runtime_error("Missing cell-vertex connectivity."); std::map input_idx_to_vertex; - for (int c = 0; c < c_to_v->num_nodes(); ++c) + for (std::size_t c = 0; c < c_to_v->num_nodes(); ++c) { auto vertices = c_to_v->links(c); std::span xdofs(xdofmap.data_handle() + c * xdofmap.extent(1), diff --git a/cpp/dolfinx/mesh/Geometry.h b/cpp/dolfinx/mesh/Geometry.h index fca6a7db46c..08a29b50c3f 100644 --- a/cpp/dolfinx/mesh/Geometry.h +++ b/cpp/dolfinx/mesh/Geometry.h @@ -234,14 +234,15 @@ Geometry(std::shared_ptr, U&&, /// @return A mesh geometry. template Geometry> -create_geometry(const Topology& topology, - const std::vector>>& elements, - std::span nodes, - std::span xdofs, const U& x, int dim, - const std::function( - const graph::AdjacencyList&)>& reorder_fn - = nullptr) +create_geometry( + const Topology& topology, + const std::vector>>& elements, + std::span nodes, std::span xdofs, + const U& x, int dim, + const std::function( + const graph::AdjacencyList>&)>& reorder_fn + = nullptr) { spdlog::info("Create Geometry (multiple)"); diff --git a/cpp/dolfinx/mesh/MeshTags.h b/cpp/dolfinx/mesh/MeshTags.h index eb7c621fe97..2e40c73dfc5 100644 --- a/cpp/dolfinx/mesh/MeshTags.h +++ b/cpp/dolfinx/mesh/MeshTags.h @@ -135,9 +135,10 @@ class MeshTags /// @note Entities that do not exist on this rank are ignored. /// @warning `entities` must not contain duplicate entities. template -MeshTags create_meshtags(std::shared_ptr topology, int dim, - const graph::AdjacencyList& entities, - std::span values) +MeshTags +create_meshtags(std::shared_ptr topology, int dim, + const graph::AdjacencyList>& entities, + std::span values) { spdlog::info( "Building MeshTags object from tagged entities (defined by vertices)."); diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index af2ad8a363c..941289ce0d1 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -39,7 +39,7 @@ namespace /// for. /// @return Map from global index to sharing ranks for each index in /// indices. The owner rank is the first as the first in the of ranks. -graph::AdjacencyList +graph::AdjacencyList> determine_sharing_ranks(MPI_Comm comm, std::span indices) { common::Timer timer("Topology: determine shared index ownership"); @@ -367,7 +367,7 @@ std::array, 2> vertex_ownership_groups( /// 3. MPI rank of the owner std::vector exchange_indexing(MPI_Comm comm, std::span indices, - const graph::AdjacencyList& index_to_ranks, + const graph::AdjacencyList>& index_to_ranks, std::int64_t offset, std::span global_indices, std::span local_indices) @@ -376,7 +376,7 @@ exchange_indexing(MPI_Comm comm, std::span indices, // Build src and destination ranks std::vector src, dest; - for (std::int32_t i = 0; i < index_to_ranks.num_nodes(); ++i) + for (std::size_t i = 0; i < index_to_ranks.num_nodes(); ++i) { if (auto ranks = index_to_ranks.links(i); ranks.front() == mpi_rank) dest.insert(dest.end(), std::next(ranks.begin()), ranks.end()); @@ -399,7 +399,7 @@ exchange_indexing(MPI_Comm comm, std::span indices, // Pack send data. Use std::vector> since size will be // modest (equal to number of neighbour ranks) std::vector> send_buffer(dest.size()); - for (std::int32_t i = 0; i < index_to_ranks.num_nodes(); ++i) + for (std::size_t i = 0; i < index_to_ranks.num_nodes(); ++i) { // Get (global) ranks that share this vertex. Note that first rank // is the owner. @@ -743,7 +743,9 @@ Topology::Topology( std::vector cell_types, std::shared_ptr vertex_map, std::vector> cell_maps, - std::vector>> cells, + std::vector< + std::shared_ptr>>> + cells, const std::optional>>& original_index, int num_threads) : original_cell_index(original_index @@ -762,7 +764,7 @@ Topology::Topology( _index_maps.insert({{0, 0}, vertex_map}); _connectivity.insert( {{{0, 0}, {0, 0}}, - std::make_shared>( + std::make_shared>>( vertex_map->size_local() + vertex_map->num_ghosts())}); if (tdim > 0) { @@ -837,7 +839,7 @@ std::shared_ptr Topology::index_map(int dim) const return im.at(0); } //----------------------------------------------------------------------------- -std::shared_ptr> +std::shared_ptr>> Topology::connectivity(std::array d0, std::array d1) const { if (auto it = _connectivity.find({d0, d1}); it == _connectivity.end()) @@ -846,7 +848,7 @@ Topology::connectivity(std::array d0, std::array d1) const return it->second; } //----------------------------------------------------------------------------- -std::shared_ptr> +std::shared_ptr>> Topology::connectivity(int d0, int d1) const { if (this->entity_types(d0).size() > 1 or this->entity_types(d0).size() > 1) @@ -1066,7 +1068,7 @@ Topology mesh::create_topology( // For each vertex whose ownership needs determining, find the sharing // ranks. The first index in the list of ranks for a vertex is the // owner (as determined by determine_sharing_ranks). - const graph::AdjacencyList global_vertex_to_ranks + const graph::AdjacencyList> global_vertex_to_ranks = determine_sharing_ranks(comm, boundary_vertices); // Iterate over vertices that have 'unknown' ownership, and if flagged @@ -1281,13 +1283,16 @@ Topology mesh::create_topology( static_cast(dolfinx::MPI::tag::consensus_nbx) + cell_types.size()); // Set cell index map and connectivity - std::vector>> cells_c; + std::vector>>> + cells_c; cells_c.reserve(cell_types.size()); for (std::size_t i = 0; i < cell_types.size(); ++i) { - cells_c.push_back(std::make_shared>( - graph::regular_adjacency_list(std::move(_cells_local_idx[i]), - mesh::num_cell_vertices(cell_types[i])))); + cells_c.push_back( + std::make_shared>>( + graph::regular_adjacency_list( + std::move(_cells_local_idx[i]), + mesh::num_cell_vertices(cell_types[i])))); } // Save original cell index @@ -1386,8 +1391,9 @@ mesh::create_subtopology(const Topology& topology, int dim, sub_e_to_v_offsets.push_back(sub_e_to_v_vec.size()); } - auto sub_e_to_v = std::make_shared>( - std::move(sub_e_to_v_vec), std::move(sub_e_to_v_offsets)); + auto sub_e_to_v + = std::make_shared>>( + std::move(sub_e_to_v_vec), std::move(sub_e_to_v_offsets)); return {Topology({entity_type}, submap0, {submap}, {sub_e_to_v}), std::move(subentities), std::move(subvertices0)}; @@ -1495,7 +1501,7 @@ mesh::compute_mixed_cell_pairs(const Topology& topology, { if (fci) { - for (std::int32_t k = 0; k < fci->num_nodes(); ++k) + for (std::size_t k = 0; k < fci->num_nodes(); ++k) { if (fci->num_links(k) == 2) { @@ -1517,7 +1523,7 @@ mesh::compute_mixed_cell_pairs(const Topology& topology, if (fci and fcj) { assert(fci->num_nodes() == fcj->num_nodes()); - for (std::int32_t k = 0; k < fci->num_nodes(); ++k) + for (std::size_t k = 0; k < fci->num_nodes(); ++k) { if (fci->num_links(k) == 1 and fcj->num_links(k) == 1) { diff --git a/cpp/dolfinx/mesh/Topology.h b/cpp/dolfinx/mesh/Topology.h index 1b4beda3102..47eb35aee6b 100644 --- a/cpp/dolfinx/mesh/Topology.h +++ b/cpp/dolfinx/mesh/Topology.h @@ -58,15 +58,16 @@ class Topology /// @param[in] original_cell_index Original indices for each cell in /// `cells`. /// @param[in] num_threads Number of threads to use for entity creation. - Topology( - std::vector cell_types, - std::shared_ptr vertex_map, - std::vector> cell_maps, - std::vector>> cells, - const std::optional>>& - original_cell_index - = std::nullopt, - int num_threads = 1); + Topology(std::vector cell_types, + std::shared_ptr vertex_map, + std::vector> cell_maps, + std::vector< + std::shared_ptr>>> + cells, + const std::optional>>& + original_cell_index + = std::nullopt, + int num_threads = 1); /// Copy constructor Topology(const Topology& topology) = default; @@ -127,7 +128,7 @@ class Topology /// incident "entity type" within topological dimension). /// @return AdjacencyList of connectivity from entity type in `d0` to /// entity types in `d1`, or `nullptr` if not yet computed. - std::shared_ptr> + std::shared_ptr>> connectivity(std::array d0, std::array d1) const; /// @brief Return connectivity from entities of dimension `d0` to @@ -139,7 +140,7 @@ class Topology /// @return The adjacency list that for each entity of dimension `d0` /// gives the list of incident entities of dimension `d1`. Returns /// `nullptr` if connectivity has not been computed. - std::shared_ptr> + std::shared_ptr>> connectivity(int d0, int d1) const; /// @brief Returns the permutation information. @@ -232,7 +233,7 @@ class Topology // where dim0 and dim1 are topological dimensions and i0 and i1 // are the indices of cell types (following the order in _entity_types). std::map, std::array>, - std::shared_ptr>> + std::shared_ptr>>> _connectivity; // The facet permutations (local facet, cell)) diff --git a/cpp/dolfinx/mesh/cell_types.cpp b/cpp/dolfinx/mesh/cell_types.cpp index bcebe6d3346..d180aeb7721 100644 --- a/cpp/dolfinx/mesh/cell_types.cpp +++ b/cpp/dolfinx/mesh/cell_types.cpp @@ -61,21 +61,22 @@ mesh::CellType mesh::to_type(const std::string& cell) throw std::runtime_error("Unknown cell type (" + cell + ")"); } //----------------------------------------------------------------------------- -graph::AdjacencyList mesh::get_entity_vertices(CellType type, int dim) +graph::AdjacencyList> mesh::get_entity_vertices(CellType type, + int dim) { std::vector> topology = basix::cell::topology(cell_type_to_basix_type(type))[dim]; - return graph::AdjacencyList(topology); + return graph::AdjacencyList>(topology); } //----------------------------------------------------------------------------- -graph::AdjacencyList mesh::get_sub_entities(CellType type, int dim0, - int dim1) +graph::AdjacencyList> +mesh::get_sub_entities(CellType type, int dim0, int dim1) { // keep backward compatibility if (type == CellType::interval) - return graph::AdjacencyList(0); + return graph::AdjacencyList>(0); else if (type == CellType::point) - return graph::AdjacencyList(0); + return graph::AdjacencyList>(0); std::vector>> connectivity = basix::cell::sub_entity_connectivity( @@ -84,7 +85,7 @@ graph::AdjacencyList mesh::get_sub_entities(CellType type, int dim0, subset.reserve(connectivity.size()); for (auto& row : connectivity) subset.emplace_back(row[dim1]); - return graph::AdjacencyList(subset); + return graph::AdjacencyList>(subset); } //----------------------------------------------------------------------------- int mesh::cell_num_entities(CellType type, int dim) @@ -108,7 +109,8 @@ mesh::cell_entity_closure(CellType cell_type) for (int i = 0; i <= cell_dim; ++i) num_entities[i] = cell_num_entities(cell_type, i); - const graph::AdjacencyList edge_v = get_entity_vertices(cell_type, 1); + const graph::AdjacencyList> edge_v + = get_entity_vertices(cell_type, 1); const auto face_e = get_sub_entities(cell_type, 2, 1); std::map, std::vector>> entity_closure; diff --git a/cpp/dolfinx/mesh/cell_types.h b/cpp/dolfinx/mesh/cell_types.h index 98ca16fa101..855e1fe4a62 100644 --- a/cpp/dolfinx/mesh/cell_types.h +++ b/cpp/dolfinx/mesh/cell_types.h @@ -122,11 +122,13 @@ inline CellType cell_entity_type(CellType type, int d, int index) /// Return list of entities, where entities(e, k) is the local vertex /// index for the kth vertex of entity e of dimension dim -graph::AdjacencyList get_entity_vertices(CellType type, int dim); +graph::AdjacencyList> get_entity_vertices(CellType type, + int dim); /// Get entities of dimension dim1 and that make up entities of dimension /// dim0. -graph::AdjacencyList get_sub_entities(CellType type, int dim0, int dim1); +graph::AdjacencyList> get_sub_entities(CellType type, int dim0, + int dim1); /// @brief Number of entities of dimension. /// diff --git a/cpp/dolfinx/mesh/generation.h b/cpp/dolfinx/mesh/generation.h index 11437a882a5..ec0ff52b25b 100644 --- a/cpp/dolfinx/mesh/generation.h +++ b/cpp/dolfinx/mesh/generation.h @@ -402,12 +402,11 @@ Mesh build_tet(MPI_Comm comm, MPI_Comm subcomm, } template -mesh::Mesh -build_hex(MPI_Comm comm, MPI_Comm subcomm, std::array, 2> p, - std::array n, - const CellPartitionFunction& partitioner, - const std::function( - const graph::AdjacencyList&)>& reorder_fn) +mesh::Mesh build_hex( + MPI_Comm comm, MPI_Comm subcomm, std::array, 2> p, + std::array n, const CellPartitionFunction& partitioner, + const std::function( + const graph::AdjacencyList>&)>& reorder_fn) { common::Timer timer("Build BoxMesh (hexahedra)"); std::vector x; diff --git a/cpp/dolfinx/mesh/graphbuild.cpp b/cpp/dolfinx/mesh/graphbuild.cpp index a1180fa5cd6..a94b9ea1e88 100644 --- a/cpp/dolfinx/mesh/graphbuild.cpp +++ b/cpp/dolfinx/mesh/graphbuild.cpp @@ -48,11 +48,11 @@ namespace /// /// @return Global dual graph, including ghost edges (edges to /// off-procss cells) -graph::AdjacencyList compute_nonlocal_dual_graph( +graph::AdjacencyList> compute_nonlocal_dual_graph( const MPI_Comm comm, std::span facets, std::size_t local_max_vertices_per_facet, std::span cells, - const graph::AdjacencyList& local_dual_graph) + const graph::AdjacencyList>& local_dual_graph) { spdlog::info("Build nonlocal part of mesh dual graph"); common::Timer timer("Compute non-local part of mesh dual graph"); @@ -77,7 +77,8 @@ graph::AdjacencyList compute_nonlocal_dual_graph( return graph::AdjacencyList( std::vector(local_dual_graph.array().begin(), local_dual_graph.array().end()), - local_dual_graph.offsets()); + std::vector(local_dual_graph.offsets().begin(), + local_dual_graph.offsets().end())); } // Postoffice (PO) setup: @@ -468,7 +469,7 @@ graph::AdjacencyList compute_nonlocal_dual_graph( std::vector disp = offsets; // Copy local data and add cell offset - for (std::int32_t i = 0; i < local_dual_graph.num_nodes(); ++i) + for (std::size_t i = 0; i < local_dual_graph.num_nodes(); ++i) { auto e = local_dual_graph.links(i); disp[i] += e.size(); @@ -519,8 +520,8 @@ graph::AdjacencyList compute_nonlocal_dual_graph( //----------------------------------------------------------------------------- } // namespace //----------------------------------------------------------------------------- -std::tuple, std::vector, - std::size_t, std::vector> +std::tuple>, + std::vector, std::size_t, std::vector> mesh::build_local_dual_graph( std::span celltypes, const std::vector>& cells, @@ -536,8 +537,8 @@ mesh::build_local_dual_graph( ncells_local == 0) { // Empty mesh on this process - return {graph::AdjacencyList(0), std::vector(), - 0, std::vector()}; + return {graph::AdjacencyList>(0), + std::vector(), 0, std::vector()}; } if (cells.size() != celltypes.size()) @@ -572,12 +573,12 @@ mesh::build_local_dual_graph( cell_offsets.push_back(cell_offsets.back() + num_cells); facet_count += num_cell_facets * num_cells; - graph::AdjacencyList cell_facets + graph::AdjacencyList> cell_facets = mesh::get_entity_vertices(cell_type, tdim - 1); // Determine/update maximum number of vertices for facet std::ranges::for_each( - std::views::iota(0, cell_facets.num_nodes()), + std::views::iota(0, (int)cell_facets.num_nodes()), [&max = max_vertices_per_facet, &cell_facets](auto node) { max = std::max(max, cell_facets.num_links(node)); }); } @@ -603,14 +604,14 @@ mesh::build_local_dual_graph( int num_cell_vertices = mesh::cell_num_entities(cell_type, 0); std::int32_t num_cells = _cells.size() / num_cell_vertices; - graph::AdjacencyList cell_facets + graph::AdjacencyList> cell_facets = mesh::get_entity_vertices(cell_type, tdim - 1); for (std::int32_t c = 0; c < num_cells; ++c) { // Loop over cell facets std::span v = _cells.subspan(num_cell_vertices * c, num_cell_vertices); - for (int f = 0; f < cell_facets.num_nodes(); ++f) + for (std::size_t f = 0; f < cell_facets.num_nodes(); ++f) { std::span facet_vertices = cell_facets.links(f); std::ranges::transform(facet_vertices, std::back_inserter(facets), @@ -728,7 +729,7 @@ mesh::build_local_dual_graph( std::move(local_cells)}; } //----------------------------------------------------------------------------- -graph::AdjacencyList +graph::AdjacencyList> mesh::build_dual_graph(MPI_Comm comm, std::span celltypes, const std::vector>& cells, std::optional max_facet_to_cell_links) diff --git a/cpp/dolfinx/mesh/graphbuild.h b/cpp/dolfinx/mesh/graphbuild.h index 8fedf5747d4..10adb3dcb26 100644 --- a/cpp/dolfinx/mesh/graphbuild.h +++ b/cpp/dolfinx/mesh/graphbuild.h @@ -56,8 +56,8 @@ enum class CellType : std::int8_t; /// @note Facet (2) and cell (4) data will contain multiple entries for /// the same facet for branching meshes with `max_facet_to_cell_links>2` /// to account for all facet cell connectivies. -std::tuple, std::vector, - std::size_t, std::vector> +std::tuple>, + std::vector, std::size_t, std::vector> build_local_dual_graph(std::span celltypes, const std::vector>& cells, std::optional max_facet_to_cell_links); @@ -93,7 +93,7 @@ build_local_dual_graph(std::span celltypes, /// higher branching) across process boundaries to be picked up by the /// dual graph. If the joints do not live on the process boundary this /// is not a problem. -graph::AdjacencyList +graph::AdjacencyList> build_dual_graph(MPI_Comm comm, std::span celltypes, const std::vector>& cells, std::optional max_facet_to_cell_links); diff --git a/cpp/dolfinx/mesh/permutationcomputation.cpp b/cpp/dolfinx/mesh/permutationcomputation.cpp index 58fcf2c2b69..c81ca97089a 100644 --- a/cpp/dolfinx/mesh/permutationcomputation.cpp +++ b/cpp/dolfinx/mesh/permutationcomputation.cpp @@ -134,8 +134,12 @@ compute_triangle_quad_face_permutations(const mesh::Topology& topology, cell_face_types[i] = mesh::cell_facet_type(cell_type, i); // Connectivity for each face type - std::vector>> c_to_f; - std::vector>> f_to_v; + std::vector< + std::shared_ptr>>> + c_to_f; + std::vector< + std::shared_ptr>>> + f_to_v; // Create mapping for each face type to cell-local face index int tdim = topology.dim(); @@ -234,7 +238,7 @@ compute_edge_reflections(const mesh::Topology& topology) std::vector> edge_perm(num_cells, 0); std::vector cell_vertices, vertices; - for (int c = 0; c < c_to_v->num_nodes(); ++c) + for (std::size_t c = 0; c < c_to_v->num_nodes(); ++c) { cell_vertices.resize(c_to_v->num_links(c)); im->local_to_global(c_to_v->links(c), cell_vertices); diff --git a/cpp/dolfinx/mesh/topologycomputation.cpp b/cpp/dolfinx/mesh/topologycomputation.cpp index d2753f4f0e8..97dcb4abbc1 100644 --- a/cpp/dolfinx/mesh/topologycomputation.cpp +++ b/cpp/dolfinx/mesh/topologycomputation.cpp @@ -41,9 +41,8 @@ namespace /// /// This code is thread-safe. /// -/// @param[in] c0 Starting cell index. -/// @param[in] num_cells Number of cells to process. /// @param[in,out] entity_list +/// @param[in,out] entity_list_sorted /// @param[in] cells Cell-to-vertex connectivity. /// @param[in] e_vertices Entity-to-vertices, where /// `e_vertices.links(e)[i]` is the `i`th local (to the cell) vertex @@ -53,11 +52,14 @@ namespace /// @param[in] vertex_index_map Index map for the vertices. void build_entity_list(std::span entity_list, std::span entity_list_sorted, - std::int32_t c0, std::int32_t num_cells, - const graph::AdjacencyList& cells, - const graph::AdjacencyList& e_vertices, + graph::AdjacencyList, + std::span> + cells, + graph::AdjacencyList, + std::span> + e_vertices, mesh::CellType entity_type, - const std::vector& cell_type_entities, + std::span cell_type_entities, const common::IndexMap& vertex_index_map) { int num_vertices_per_entity = num_cell_vertices(entity_type); @@ -68,10 +70,10 @@ void build_entity_list(std::span entity_list, std::vector perm(num_vertices_per_entity); // Iterate over cells - for (std::int32_t c = 0; c < num_cells; ++c) + for (std::size_t c = 0; c < cells.num_nodes(); ++c) { // Get vertices from each cell - auto vertices = cells.links(c + c0); + auto vertices = cells.links(c); // Iterate over cell entities of given type for (int e = 0; e < num_entities_per_cell; ++e) @@ -81,20 +83,20 @@ void build_entity_list(std::span entity_list, // Get entity vertices. Padded with -1 if fewer than // max_vertices_per_entity // - // NOTE: Entity orientation is determined by vertex - // ordering. The orientation of an entity with respect to - // the cell may differ from its global mesh orientation. - // Hence, we reorder the vertices so that each entity's - // orientation agrees with their global orientation. + // NOTE: Entity orientation is determined by vertex ordering. The + // orientation of an entity with respect to the cell may differ + // from its global mesh orientation. Hence, we reorder the + // vertices so that each entity's orientation agrees with their + // global orientation. // - // FIXME: This might be better below when the entity to - // vertex connectivity is computed + // FIXME: This might be better below when the entity to vertex + // connectivity is computed assert(ev.size() == entity_vertices.size()); for (std::size_t j = 0; j < ev.size(); ++j) entity_vertices[j] = vertices[ev[j]]; - // Orient the entities. Simply sort according to global - // vertex index for simplices. + // Orient the entities. Simply sort according to global vertex + // index for simplices. assert(entity_vertices.size() == global_vertices.size()); vertex_index_map.local_to_global(entity_vertices, global_vertices); @@ -102,8 +104,8 @@ void build_entity_list(std::span entity_list, std::ranges::sort(perm, [&global_vertices](std::size_t i0, std::size_t i1) { return global_vertices[i0] < global_vertices[i1]; }); - // For quadrilaterals, the vertex opposite the lowest - // vertex should be last + // For quadrilaterals, the vertex opposite the lowest vertex + // should be last if (entity_type == mesh::CellType::quadrilateral) { std::size_t min_vertex_idx = perm[0]; @@ -127,16 +129,17 @@ void build_entity_list(std::span entity_list, } } -/// @brief Create an adjacency list from array of pairs, where -/// the first value in the pair is the node and the second value -/// is the edge. +/// @brief Create an adjacency list from array of pairs, where the first +/// value in the pair is the node and the second value is the edge. +/// /// @param[in] data List if pairs -/// @param[in] size The number of edges in the graph. For -/// example, this can be used to build an adjacency list that -/// includes 'owned' nodes only. +/// @param[in] size The number of edges in the graph. For example, this +/// can be used to build an adjacency list that includes 'owned' nodes +/// only. /// @pre The `data` array must be sorted. template -graph::AdjacencyList create_adj_list(U& data, std::int32_t size) +graph::AdjacencyList> create_adj_list(U& data, + std::int32_t size) { std::ranges::sort(data); auto [unique_end, range_end] = std::ranges::unique(data); @@ -169,9 +172,8 @@ graph::AdjacencyList create_adj_list(U& data, std::int32_t size) template int get_ownership(const U& processes, const V& vertices) { - // Use a deterministic random number generator, seeded with - // global vertex indices ensuring all processes get the same - // answer + // Use a deterministic random number generator, seeded with global + // vertex indices ensuring all processes get the same answer std::mt19937 gen; std::seed_seq seq(vertices.begin(), vertices.end()); gen.seed(seq); @@ -181,16 +183,17 @@ int get_ownership(const U& processes, const V& vertices) return owner; } //----------------------------------------------------------------------------- -/// Communicate with sharing processes to find out which -/// entities are ghosts and return a map (vector) to move these -/// local indices to the end of the local range. Also returns -/// the index map, and shared entities, i.e. the set of all -/// processes which share each shared entity. + +/// Communicate with sharing processes to find out which / entities are +/// ghosts and return a map (vector) to move these / local indices to +/// the end of the local range. Also returns / the index map, and shared +/// entities, i.e. the set of all / processes which share each shared +/// entity. /// @param[in] comm MPI Communicator /// @param[in] cell_map Index map for cell distribution /// @param[in] vertex_map Index map for vertex distribution -/// @param[in] entity_list List of entities, each entity -/// represented by its local vertex indices +/// @param[in] entity_list List of entities, each entity / represented +/// by its local vertex indices /// @param[in] num_vertices_per_e Number of vertices per entity /// @param[in] num_entities_per_cell Number of entities per cell /// @param[in] entity_index Initial numbering for each row in @@ -225,7 +228,8 @@ get_local_indexing(MPI_Comm comm, const common::IndexMap& vertex_map, // Create a symmetric neighbor_comm from vertex_ranks // Get sharing ranks for each vertex - graph::AdjacencyList vertex_ranks = vertex_map.index_to_dest_ranks(); + graph::AdjacencyList> vertex_ranks + = vertex_map.index_to_dest_ranks(); // Create unique list of ranks that share vertices (owners of) std::vector ranks(vertex_ranks.array().begin(), @@ -421,10 +425,10 @@ get_local_indexing(MPI_Comm comm, const common::IndexMap& vertex_map, } } - const graph::AdjacencyList shared_entities + const graph::AdjacencyList> shared_entities = create_adj_list(shared_entities_data, entity_count); - const graph::AdjacencyList shared_entities_v + const graph::AdjacencyList> shared_entities_v = create_adj_list(shared_entity_to_global_vertices_data, entity_count); //--------- @@ -556,15 +560,16 @@ get_local_indexing(MPI_Comm comm, const common::IndexMap& vertex_map, /// @return Returns the (cell-entity connectivity, entity-vertex /// connectivity, index map for the entity distribution across /// processes, shared entities) -std::tuple>>, - graph::AdjacencyList, common::IndexMap, +std::tuple>>>, + graph::AdjacencyList>, common::IndexMap, std::vector> compute_entities_by_key_matching( MPI_Comm comm, - std::vector>, - std::reference_wrapper>> + std::vector>>, + std::reference_wrapper>> cell_lists, const common::IndexMap& vertex_index_map, mesh::CellType entity_type, int dim, int num_threads) @@ -591,7 +596,7 @@ compute_entities_by_key_matching( cell_type_entities[k].push_back(e); } - const graph::AdjacencyList& cells + const graph::AdjacencyList>& cells = std::get<1>(cell_lists[k]); std::size_t num_cells = cells.num_nodes(); cell_type_offsets.push_back(cell_type_offsets.back() @@ -615,7 +620,7 @@ compute_entities_by_key_matching( common::Timer t_thread("Threaded part"); - const graph::AdjacencyList& cells + const graph::AdjacencyList>& cells = std::get<1>(cell_lists[k]); int num_entities_per_cell = cell_type_entities[k].size(); std::vector threads(num_threads); @@ -623,7 +628,6 @@ compute_entities_by_key_matching( { auto [c0, c1] = dolfinx::MPI::local_range(i, cells.num_nodes(), num_threads); - std::size_t offset = cell_type_offsets[k] * num_vertices_per_entity + c0 * num_vertices_per_entity * num_entities_per_cell; @@ -631,9 +635,9 @@ compute_entities_by_key_matching( = (c1 - c0) * num_vertices_per_entity * num_entities_per_cell; threads[i] = std::jthread( build_entity_list, std::span(entity_list.data() + offset, count), - std::span(entity_list_sorted.data() + offset, count), c0, c1 - c0, - std::cref(cells), std::cref(e_vertices), entity_type, - std::cref(cell_type_entities[k]), std::cref(vertex_index_map)); + std::span(entity_list_sorted.data() + offset, count), + cells.sub(c0, c1), e_vertices, entity_type, cell_type_entities[k], + std::cref(vertex_index_map)); } } @@ -687,7 +691,7 @@ compute_entities_by_key_matching( // which only appear in ghost cells tagged. const common::IndexMap& cell_map = std::get<2>(cell_lists[k]); assert(cell_map.size_local() + cell_map.num_ghosts() - == std::get<1>(cell_lists[k]).get().num_nodes()); + == (std::int32_t)std::get<1>(cell_lists[k]).get().num_nodes()); const std::int32_t ghost_offset = cell_map.size_local(); int num_entities_per_cell = cell_type_entities[k].size(); std::size_t offset = cell_type_offsets[k]; @@ -716,15 +720,15 @@ compute_entities_by_key_matching( num_vertices_per_entity, ev.links(local_index[i]).begin()); } - std::vector>> ce( - cell_lists.size()); + std::vector>>> + ce(cell_lists.size()); for (std::size_t k = 0; k < cell_lists.size(); ++k) { if (!cell_type_entities[k].empty()) { std::vector tmp(std::next(local_index.begin(), cell_type_offsets[k]), std::next(local_index.begin(), cell_type_offsets[k + 1])); - ce[k] = std::make_shared>( + ce[k] = std::make_shared>>( graph::regular_adjacency_list(std::move(tmp), cell_type_entities[k].size())); } @@ -745,18 +749,16 @@ compute_entities_by_key_matching( /// dimension d0. /// @return The connectivity from entities of dimension d0 to /// entities of dimension d1. -graph::AdjacencyList -compute_from_transpose(const graph::AdjacencyList& c_d1_d0, - const int num_entities_d0) +graph::AdjacencyList> compute_from_transpose( + const graph::AdjacencyList>& c_d1_d0, + const int num_entities_d0) { // Compute number of connections for each e0 std::vector num_connections(num_entities_d0, 0); - for (int e1 = 0; e1 < c_d1_d0.num_nodes(); ++e1) - { + for (std::size_t e1 = 0; e1 < c_d1_d0.num_nodes(); ++e1) for (std::int32_t e0 : c_d1_d0.links(e1)) num_connections[e0]++; - } // Compute offsets std::vector offsets(num_connections.size() + 1, 0); @@ -765,7 +767,7 @@ compute_from_transpose(const graph::AdjacencyList& c_d1_d0, std::vector counter(num_connections.size(), 0); std::vector connections(offsets[offsets.size() - 1]); - for (int e1 = 0; e1 < c_d1_d0.num_nodes(); ++e1) + for (std::size_t e1 = 0; e1 < c_d1_d0.num_nodes(); ++e1) for (std::int32_t e0 : c_d1_d0.links(e1)) connections[offsets[e0] + counter[e0]++] = e1; @@ -781,16 +783,16 @@ compute_from_transpose(const graph::AdjacencyList& c_d1_d0, /// @param[in] cell_type_d0 The cell type for entities of /// dimension d0 /// @return The d0 -> d1 connectivity -graph::AdjacencyList -compute_from_map(const graph::AdjacencyList& c_d0_0, - const graph::AdjacencyList& c_d1_0) +graph::AdjacencyList> +compute_from_map(const graph::AdjacencyList>& c_d0_0, + const graph::AdjacencyList>& c_d1_0) { // Make a map from the sorted edge vertices to the edge index boost::unordered_map, std::int32_t> edge_to_index; edge_to_index.reserve(c_d1_0.num_nodes()); std::array key; - for (int e = 0; e < c_d1_0.num_nodes(); ++e) + for (std::size_t e = 0; e < c_d1_0.num_nodes(); ++e) { std::span v = c_d1_0.links(e); assert(v.size() == key.size()); @@ -802,14 +804,15 @@ compute_from_map(const graph::AdjacencyList& c_d0_0, // vertices so AdjacencyList will have same offset pattern std::vector connections; connections.reserve(c_d0_0.array().size()); - std::vector offsets(c_d0_0.offsets()); + std::vector offsets(c_d0_0.offsets().begin(), + c_d0_0.offsets().end()); // Search for edges of facet in map, and recover index - const graph::AdjacencyList tri_vertices_ref + const graph::AdjacencyList> tri_vertices_ref = get_entity_vertices(mesh::CellType::triangle, 1); - const graph::AdjacencyList quad_vertices_ref + const graph::AdjacencyList> quad_vertices_ref = get_entity_vertices(mesh::CellType::quadrilateral, 1); - for (int e = 0; e < c_d0_0.num_nodes(); ++e) + for (std::size_t e = 0; e < c_d0_0.num_nodes(); ++e) { auto e0 = c_d0_0.links(e); auto vref = (e0.size() == 3) ? &tri_vertices_ref : &quad_vertices_ref; @@ -832,8 +835,9 @@ compute_from_map(const graph::AdjacencyList& c_d0_0, } // namespace //----------------------------------------------------------------------------- -std::tuple>>, - std::shared_ptr>, +std::tuple>>>, + std::shared_ptr>>, std::shared_ptr, std::vector> mesh::compute_entities(const Topology& topology, int dim, CellType entity_type, int num_threads) @@ -843,8 +847,10 @@ mesh::compute_entities(const Topology& topology, int dim, CellType entity_type, // Vertices must always exist if (dim == 0) { - return {std::vector>>(), - nullptr, nullptr, std::vector()}; + return { + std::vector< + std::shared_ptr>>>(), + nullptr, nullptr, std::vector()}; } { @@ -853,9 +859,9 @@ mesh::compute_entities(const Topology& topology, int dim, CellType entity_type, int index = std::distance(topology.entity_types(dim).begin(), idx); if (topology.connectivity({dim, index}, {0, 0})) { - return { - std::vector>>(), - nullptr, nullptr, std::vector()}; + return {std::vector>>>(), + nullptr, nullptr, std::vector()}; } } @@ -863,10 +869,10 @@ mesh::compute_entities(const Topology& topology, int dim, CellType entity_type, // Lists of all cells by cell type std::vector cell_types = topology.entity_types(tdim); - std::vector>, - std::reference_wrapper>> + std::vector>>, + std::reference_wrapper>> cell_lists; auto cell_index_maps = topology.index_maps(tdim); @@ -888,12 +894,13 @@ mesh::compute_entities(const Topology& topology, int dim, CellType entity_type, topology.comm(), cell_lists, *vertex_map, entity_type, dim, num_threads); return {d0, - std::make_shared>(std::move(d1)), + std::make_shared>>( + std::move(d1)), std::make_shared(std::move(im)), std::move(interprocess_entities)}; } //----------------------------------------------------------------------------- -std::array>, 2> +std::array>>, 2> mesh::compute_connectivity(const Topology& topology, std::array d0, std::array d1) { @@ -924,7 +931,7 @@ mesh::compute_connectivity(const Topology& topology, std::array d0, } // Get entities if they exist - std::shared_ptr> c_d0_0 + std::shared_ptr>> c_d0_0 = topology.connectivity(d0, {0, 0}); if (d0[0] > 0 and !topology.connectivity(d0, {0, 0})) { @@ -932,7 +939,7 @@ mesh::compute_connectivity(const Topology& topology, std::array d0, + std::to_string(d0[0]) + "."); } - std::shared_ptr> c_d1_0 + std::shared_ptr>> c_d1_0 = topology.connectivity(d1, {0, 0}); if (d1[0] > 0 and !topology.connectivity(d1, {0, 0})) { @@ -947,7 +954,7 @@ mesh::compute_connectivity(const Topology& topology, std::array d0, // Decide how to compute the connectivity if (d0 == d1) { - return {std::make_shared>( + return {std::make_shared>>( c_d0_0->num_nodes()), nullptr}; } @@ -959,13 +966,15 @@ mesh::compute_connectivity(const Topology& topology, std::array d0, { // Only possible case is edge->facet assert(d0[0] == 1 and d1[0] == 2); - auto c_d1_d0 = std::make_shared>( - compute_from_map(*c_d1_0, *c_d0_0)); + auto c_d1_d0 + = std::make_shared>>( + compute_from_map(*c_d1_0, *c_d0_0)); spdlog::info("Computing mesh connectivity {}-{} from transpose.", d0[0], d1[0]); - auto c_d0_d1 = std::make_shared>( - compute_from_transpose(*c_d1_d0, c_d0_0->num_nodes())); + auto c_d0_d1 + = std::make_shared>>( + compute_from_transpose(*c_d1_d0, c_d0_0->num_nodes())); return {c_d0_d1, c_d1_d0}; } else @@ -975,9 +984,10 @@ mesh::compute_connectivity(const Topology& topology, std::array d0, spdlog::info("Computing mesh connectivity {}-{} from transpose.", std::to_string(d0[0]), std::to_string(d1[0])); - auto c_d0_d1 = std::make_shared>( - compute_from_transpose(*topology.connectivity(d1, d0), - c_d0_0->num_nodes())); + auto c_d0_d1 + = std::make_shared>>( + compute_from_transpose(*topology.connectivity(d1, d0), + c_d0_0->num_nodes())); return {c_d0_d1, nullptr}; } } @@ -988,8 +998,9 @@ mesh::compute_connectivity(const Topology& topology, std::array d0, // Only possible case is facet->edge assert(d0[0] == 2 and d1[0] == 1); - auto c_d0_d1 = std::make_shared>( - compute_from_map(*c_d0_0, *c_d1_0)); + auto c_d0_d1 + = std::make_shared>>( + compute_from_map(*c_d0_0, *c_d1_0)); return {c_d0_d1, nullptr}; } else diff --git a/cpp/dolfinx/mesh/topologycomputation.h b/cpp/dolfinx/mesh/topologycomputation.h index be0550c31c4..584a826ab82 100644 --- a/cpp/dolfinx/mesh/topologycomputation.h +++ b/cpp/dolfinx/mesh/topologycomputation.h @@ -43,8 +43,9 @@ class Topology; /// owned cells of each process. If entities of type `entity_type` /// already exists, then {nullptr, nullptr, nullptr, std::vector()} is /// returned. -std::tuple>>, - std::shared_ptr>, +std::tuple>>>, + std::shared_ptr>>, std::shared_ptr, std::vector> compute_entities(const Topology& topology, int dim, CellType entity_type, int num_threads = 1); @@ -61,7 +62,7 @@ compute_entities(const Topology& topology, int dim, CellType entity_type, /// (d0, d1) is computed and the computation of (d1, d0) was required as /// part of computing (d0, d1), the (d1, d0) is returned as the second /// entry. The second entry is otherwise nullptr. -std::array>, 2> +std::array>>, 2> compute_connectivity(const Topology& topology, std::array d0, std::array d1); diff --git a/cpp/dolfinx/mesh/utils.cpp b/cpp/dolfinx/mesh/utils.cpp index 9a3e3f1db7c..1f768734ec4 100644 --- a/cpp/dolfinx/mesh/utils.cpp +++ b/cpp/dolfinx/mesh/utils.cpp @@ -104,7 +104,7 @@ mesh::CellPartitionFunction mesh::create_cell_partitioner( return [partfn, ghost_mode, max_facet_to_cell_links]( MPI_Comm comm, int nparts, const std::vector& cell_types, const std::vector>& cells) - -> graph::AdjacencyList + -> graph::AdjacencyList> { spdlog::info("Compute partition of cells across ranks"); diff --git a/cpp/dolfinx/mesh/utils.h b/cpp/dolfinx/mesh/utils.h index cfde69e8ba6..52f13d116df 100644 --- a/cpp/dolfinx/mesh/utils.h +++ b/cpp/dolfinx/mesh/utils.h @@ -206,16 +206,17 @@ std::vector exterior_facet_indices(const Topology& topology); /// points, should not be included. /// @return Destination ranks for each cell on this process. /// @note Cells can have multiple destination ranks, when ghosted. -using CellPartitionFunction = std::function( - MPI_Comm comm, int nparts, const std::vector& cell_types, - const std::vector>& cells)>; +using CellPartitionFunction + = std::function>( + MPI_Comm comm, int nparts, const std::vector& cell_types, + const std::vector>& cells)>; /// @brief Function that reorders (locally) cells that /// are owned by this process. It takes the local mesh dual graph as an /// argument and returns a list whose `i`th entry is the new index of /// cell `i`. using CellReorderFunction = std::function( - const graph::AdjacencyList&)>; + const graph::AdjacencyList>&)>; /// @brief Creates the default boundary vertices routine for a given reorder /// function. @@ -635,7 +636,7 @@ compute_vertex_coords(const mesh::Mesh& mesh) auto x_dofmap = mesh.geometry().dofmap(cell_type_idx); auto c_to_v = topology->connectivity({tdim, cell_type_idx}, {0, 0}); assert(c_to_v); - for (int c = 0; c < c_to_v->num_nodes(); ++c) + for (std::size_t c = 0; c < c_to_v->num_nodes(); ++c) { auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent); auto vertices = c_to_v->links(c); @@ -1076,7 +1077,7 @@ Mesh> create_mesh( { spdlog::info("Using partitioner with cell data ({} cell types)", num_cell_types); - graph::AdjacencyList dest(0); + graph::AdjacencyList> dest(0); if (commt != MPI_COMM_NULL) { int size = dolfinx::MPI::size(comm); @@ -1107,7 +1108,7 @@ Mesh> create_mesh( std::int32_t offset_0 = offsets_i.front(); std::ranges::for_each(offsets_i, [&offset_0](std::int32_t& j) { j -= offset_0; }); - graph::AdjacencyList dest_i(data_i, offsets_i); + graph::AdjacencyList> dest_i(data_i, offsets_i); cell_offset += num_cells; // Distribute cells (topology, includes higher-order 'nodes') to diff --git a/cpp/dolfinx/refinement/interval.h b/cpp/dolfinx/refinement/interval.h index df7fcdb9993..c4208599a8b 100644 --- a/cpp/dolfinx/refinement/interval.h +++ b/cpp/dolfinx/refinement/interval.h @@ -33,7 +33,7 @@ namespace dolfinx::refinement::interval /// @return New mesh data: cell topology, vertex coordinates and parent /// cell indices. template -std::tuple, std::vector, +std::tuple>, std::vector, std::array, std::optional>, std::optional>> compute_refinement_data(const mesh::Mesh& mesh, @@ -55,11 +55,12 @@ compute_refinement_data(const mesh::Mesh& mesh, // TODO: creation of sharing ranks in external function? Also same // code in use for plaza // Get sharing ranks for each cell - graph::AdjacencyList cell_ranks = map_c->index_to_dest_ranks(); + graph::AdjacencyList> cell_ranks + = map_c->index_to_dest_ranks(); // Create unique list of ranks that share cells (owners of ghosts plus // ranks that ghost owned indices) - std::vector ranks = cell_ranks.array(); + std::vector ranks(cell_ranks.array().begin(), cell_ranks.array().end()); std::ranges::sort(ranks); auto to_remove = std::ranges::unique(ranks); ranks.erase(to_remove.begin(), to_remove.end()); diff --git a/cpp/dolfinx/refinement/plaza.cpp b/cpp/dolfinx/refinement/plaza.cpp index 8189fe321a4..5b11977beed 100644 --- a/cpp/dolfinx/refinement/plaza.cpp +++ b/cpp/dolfinx/refinement/plaza.cpp @@ -183,11 +183,10 @@ get_tetrahedra(std::span indices, } // namespace //----------------------------------------------------------------------------- -void plaza::impl::enforce_rules(MPI_Comm comm, - const graph::AdjacencyList& shared_edges, - std::span marked_edges, - const mesh::Topology& topology, - std::span long_edge) +void plaza::impl::enforce_rules( + MPI_Comm comm, const graph::AdjacencyList>& shared_edges, + std::span marked_edges, const mesh::Topology& topology, + std::span long_edge) { common::Timer t0("PLAZA: Enforce rules"); diff --git a/cpp/dolfinx/refinement/plaza.h b/cpp/dolfinx/refinement/plaza.h index bcb469a896c..181dbb8aa16 100644 --- a/cpp/dolfinx/refinement/plaza.h +++ b/cpp/dolfinx/refinement/plaza.h @@ -128,7 +128,8 @@ get_simplices(std::span indices, /// Propagate edge markers according to rules (longest edge of each /// face must be marked, if any edge of face is marked) -void enforce_rules(MPI_Comm comm, const graph::AdjacencyList& shared_edges, +void enforce_rules(MPI_Comm comm, + const graph::AdjacencyList>& shared_edges, std::span marked_edges, const mesh::Topology& topology, std::span long_edge); @@ -217,7 +218,7 @@ face_long_edge(const mesh::Mesh& mesh) assert(f_to_e); const std::vector global_indices = mesh.topology()->index_map(0)->global_indices(); - for (int f = 0; f < f_to_v->num_nodes(); ++f) + for (std::size_t f = 0; f < f_to_v->num_nodes(); ++f) { auto face_edges = f_to_e->links(f); @@ -282,12 +283,12 @@ face_long_edge(const mesh::Mesh& mesh) /// Shape of the new geometry_shape, (4) Map from new cells to parent cells /// and (5) map from refined facets to parent facets. template -std::tuple, std::vector, +std::tuple>, std::vector, std::array, std::optional>, std::optional>> compute_refinement(MPI_Comm neighbor_comm, std::span marked_edges, - const graph::AdjacencyList& shared_edges, + const graph::AdjacencyList>& shared_edges, const mesh::Mesh& mesh, std::span long_edge, std::span edge_ratio_ok, Option option) @@ -456,7 +457,7 @@ compute_refinement(MPI_Comm neighbor_comm, /// @return New mesh data: cell topology, vertex coordinates and parent /// cell index, and stored parent facet indices (if requested). template -std::tuple, std::vector, +std::tuple>, std::vector, std::array, std::optional>, std::optional>> compute_refinement_data(const mesh::Mesh& mesh, @@ -478,7 +479,8 @@ compute_refinement_data(const mesh::Mesh& mesh, throw std::runtime_error("Edges must be initialised"); // Get sharing ranks for each edge - graph::AdjacencyList edge_ranks = map_e->index_to_dest_ranks(); + graph::AdjacencyList> edge_ranks + = map_e->index_to_dest_ranks(); // Create unique list of ranks that share edges (owners of ghosts plus // ranks that ghost owned indices) diff --git a/cpp/dolfinx/refinement/refine.h b/cpp/dolfinx/refinement/refine.h index ca7f4c65d69..bc84bf8585a 100644 --- a/cpp/dolfinx/refinement/refine.h +++ b/cpp/dolfinx/refinement/refine.h @@ -50,7 +50,7 @@ create_identity_partitioner(const mesh::Mesh& parent_mesh, parent_cell](MPI_Comm comm, int /*nparts*/, std::vector cell_types, std::vector> cells) - -> graph::AdjacencyList + -> graph::AdjacencyList> { auto parent_cell_im = parent_mesh.topology()->index_map(parent_mesh.topology()->dim()); diff --git a/cpp/dolfinx/refinement/utils.cpp b/cpp/dolfinx/refinement/utils.cpp index af8bae95a1e..c2b8fcab970 100644 --- a/cpp/dolfinx/refinement/utils.cpp +++ b/cpp/dolfinx/refinement/utils.cpp @@ -236,8 +236,8 @@ std::array, 2> refinement::transfer_facet_meshtag( offset_child.front() = 0; std::partial_sum(count_child.begin(), count_child.end(), std::next(offset_child.begin())); - graph::AdjacencyList p_to_c_facet(std::move(child_facet), - std::move(offset_child)); + graph::AdjacencyList> p_to_c_facet( + std::move(child_facet), std::move(offset_child)); // Copy facet meshtag from parent to child std::vector facet_indices, tag_values; diff --git a/cpp/dolfinx/refinement/utils.h b/cpp/dolfinx/refinement/utils.h index 177fae24422..4f869b47f33 100644 --- a/cpp/dolfinx/refinement/utils.h +++ b/cpp/dolfinx/refinement/utils.h @@ -143,10 +143,10 @@ void update_logical_edgefunction( /// coordinates of the new vertices (row-major storage) and (2) the shape of the /// new coordinates. template -std::tuple, std::vector, +std::tuple>, std::vector, std::array> create_new_vertices(MPI_Comm comm, - const graph::AdjacencyList& shared_edges, + const graph::AdjacencyList>& shared_edges, const mesh::Mesh& mesh, std::span marked_edges) { @@ -273,7 +273,8 @@ create_new_vertices(MPI_Comm comm, offset[m + 1] = 1; std::partial_sum(offset.begin(), offset.end(), offset.begin()); - graph::AdjacencyList local_edge_to_new_vertex_adj(data, offset); + graph::AdjacencyList> local_edge_to_new_vertex_adj( + data, offset); return {std::move(local_edge_to_new_vertex_adj), std::move(new_vertex_coords), xshape}; } diff --git a/cpp/test/graph.cpp b/cpp/test/graph.cpp index 6aea18a3436..928ff7817fe 100644 --- a/cpp/test/graph.cpp +++ b/cpp/test/graph.cpp @@ -45,7 +45,7 @@ void test_adjacency_list_create() std::vector node_data{-1, 5, -20}; graph::AdjacencyList g1(edges, offsets, node_data); - CHECK(std::ranges::equal(g1.node_data().value(), node_data)); + CHECK(std::ranges::equal(g1.node_data(), node_data)); } void test_comm_graphs() diff --git a/cpp/test/mesh/generation.cpp b/cpp/test/mesh/generation.cpp index 5f588192970..6c961ffbd7e 100644 --- a/cpp/test/mesh/generation.cpp +++ b/cpp/test/mesh/generation.cpp @@ -25,13 +25,13 @@ namespace { template void CHECK_adjacency_list_equal( - const dolfinx::graph::AdjacencyList& adj_list, + const dolfinx::graph::AdjacencyList>& adj_list, const std::vector>& expected_list) { REQUIRE(static_cast(adj_list.num_nodes()) == expected_list.size()); - for (T i = 0; i < adj_list.num_nodes(); i++) + for (std::size_t i = 0; i < adj_list.num_nodes(); ++i) { CHECK_THAT(adj_list.links(i), Catch::Matchers::RangeEquals(expected_list[i])); @@ -114,7 +114,7 @@ TEMPLATE_TEST_CASE("Interval mesh (parallel)", "[mesh][interval]", float, else SKIP("Test only supports <= 3 processes"); - return graph::AdjacencyList(data); + return graph::AdjacencyList>(data); }; mesh::Mesh mesh = mesh::create_interval( diff --git a/cpp/test/mesh/refinement/interval.cpp b/cpp/test/mesh/refinement/interval.cpp index a3297d49de8..4f0b73e07e0 100644 --- a/cpp/test/mesh/refinement/interval.cpp +++ b/cpp/test/mesh/refinement/interval.cpp @@ -27,12 +27,11 @@ namespace { template void CHECK_adjacency_list_equal( - const dolfinx::graph::AdjacencyList& adj_list, + const dolfinx::graph::AdjacencyList>& adj_list, const std::vector>& expected_list) { - REQUIRE(static_cast(adj_list.num_nodes()) - == expected_list.size()); - for (T i = 0; i < adj_list.num_nodes(); i++) + REQUIRE(adj_list.num_nodes() == expected_list.size()); + for (std::size_t i = 0; i < adj_list.num_nodes(); ++i) { CHECK_THAT(adj_list.links(i), Catch::Matchers::RangeEquals(expected_list[i])); @@ -177,9 +176,9 @@ TEMPLATE_TEST_CASE("Interval Refinement (parallel)", = [](MPI_Comm /* comm */, int /* nparts */, const std::vector& /* cell_types */, const std::vector>& /* cells */) - -> graph::AdjacencyList + -> graph::AdjacencyList> { - return graph::AdjacencyList( + return graph::AdjacencyList>( dolfinx::MPI::size(MPI_COMM_WORLD)); }; diff --git a/cpp/test/mesh/refinement/rectangle.cpp b/cpp/test/mesh/refinement/rectangle.cpp index 8c0534fa951..ffd7800baaf 100644 --- a/cpp/test/mesh/refinement/rectangle.cpp +++ b/cpp/test/mesh/refinement/rectangle.cpp @@ -29,7 +29,7 @@ namespace { template void CHECK_adjacency_list_equal( - const dolfinx::graph::AdjacencyList& adj_list, + const dolfinx::graph::AdjacencyList>& adj_list, const std::vector>& expected_list) { REQUIRE(static_cast(adj_list.num_nodes()) diff --git a/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h b/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h index 3ebe2c01040..ac1c9bfad09 100644 --- a/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h +++ b/python/dolfinx/wrappers/dolfinx_wrappers/mesh.h @@ -18,7 +18,8 @@ template auto create_partitioner_cpp(Functor&& p) { return [p](MPI_Comm comm, int nparts, - const dolfinx::graph::AdjacencyList& local_graph, + const dolfinx::graph::AdjacencyList>& + local_graph, bool ghosting) { return p(dolfinx_wrappers::MPICommWrapper(comm), nparts, local_graph, @@ -43,13 +44,13 @@ auto create_cell_partitioner_py(Functor&& p) } using PythonCellPartitionFunction - = std::function( + = std::function>( dolfinx_wrappers::MPICommWrapper, int, const std::vector&, std::vector>)>; using CppCellPartitionFunction - = std::function( + = std::function>( MPI_Comm, int, const std::vector& q, const std::vector>&)>; diff --git a/python/dolfinx/wrappers/fem.cpp b/python/dolfinx/wrappers/fem.cpp index 188481af4e0..8d429662d49 100644 --- a/python/dolfinx/wrappers/fem.cpp +++ b/python/dolfinx/wrappers/fem.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -1270,8 +1271,8 @@ void fem(nb::module_& m) assert(topology.entity_types(topology.dim()).size() == 1); auto [map, bs, dofmap] = dolfinx::fem::build_dofmap_data( comm.get(), topology, {layout}, - [](const dolfinx::graph::AdjacencyList& g) - { return dolfinx::graph::reorder_gps(g); }); + [](const dolfinx::graph::AdjacencyList>& + g) { return dolfinx::graph::reorder_gps(g); }); return std::tuple(std::move(map), bs, std::move(dofmap)); }, nb::arg("comm"), nb::arg("topology"), nb::arg("layout"), @@ -1325,10 +1326,14 @@ void fem(nb::module_& m) const dolfinx::fem::ElementDofLayout& element, std::shared_ptr index_map, int index_map_bs, - const dolfinx::graph::AdjacencyList& dofmap, int bs) + nb::ndarray, nb::c_contig> dofmap, + int bs) { - new (self) dolfinx::fem::DofMap(element, index_map, index_map_bs, - dofmap.array(), bs); + new (self) dolfinx::fem::DofMap( + element, index_map, index_map_bs, + std::vector(dofmap.data(), + dofmap.data() + dofmap.size()), + bs); }, nb::arg("element_dof_layout"), nb::arg("index_map"), nb::arg("index_map_bs"), nb::arg("dofmap"), nb::arg("bs")) diff --git a/python/dolfinx/wrappers/geometry.cpp b/python/dolfinx/wrappers/geometry.cpp index 4d04f54c6a9..31f06335509 100644 --- a/python/dolfinx/wrappers/geometry.cpp +++ b/python/dolfinx/wrappers/geometry.cpp @@ -153,7 +153,7 @@ void declare_bbtree(nb::module_& m, const std::string& type) m.def( "compute_colliding_cells", [](const dolfinx::mesh::Mesh& mesh, - const dolfinx::graph::AdjacencyList& candidate_cells, + const dolfinx::graph::AdjacencyList>& candidate_cells, nb::ndarray, nb::c_contig> points) { return dolfinx::geometry::compute_colliding_cells( @@ -163,7 +163,7 @@ void declare_bbtree(nb::module_& m, const std::string& type) m.def( "compute_colliding_cells", [](const dolfinx::mesh::Mesh& mesh, - const dolfinx::graph::AdjacencyList& candidate_cells, + const dolfinx::graph::AdjacencyList>& candidate_cells, nb::ndarray, nb::c_contig> points) { return dolfinx::geometry::compute_colliding_cells( diff --git a/python/dolfinx/wrappers/graph.cpp b/python/dolfinx/wrappers/graph.cpp index b9269f260b3..2f7ed8666d5 100644 --- a/python/dolfinx/wrappers/graph.cpp +++ b/python/dolfinx/wrappers/graph.cpp @@ -1,4 +1,4 @@ -// Copyright (C) 2017-2025 Chris N. Richardson and Garth N. Wells +// Copyright (C) 2017-2026 Chris N. Richardson and Garth N. Wells // // This file is part of DOLFINx (https://www.fenicsproject.org) // @@ -34,7 +34,8 @@ template auto create_partitioner_py(Functor&& p_cpp) { return [p_cpp](dolfinx_wrappers::MPICommWrapper comm, int nparts, - const dolfinx::graph::AdjacencyList& local_graph, + const dolfinx::graph::AdjacencyList>& + local_graph, bool ghosting) { return p_cpp(comm.get(), nparts, local_graph, ghosting); }; } @@ -43,45 +44,54 @@ template void declare_adjacency_list_init(nb::module_& m, std::string type) { std::string pyclass_name = std::string("AdjacencyList_") + type; - nb::class_>(m, pyclass_name.c_str(), - "Adjacency List") + nb::class_, + std::vector, U>>( + m, pyclass_name.c_str(), "Adjacency List") .def( "__init__", - [](dolfinx::graph::AdjacencyList* a, + [](dolfinx::graph::AdjacencyList, + std::vector, U>* a, nb::ndarray, nb::c_contig> adj) { std::vector data(adj.data(), adj.data() + adj.size()); - new (a) dolfinx::graph::AdjacencyList( - dolfinx::graph::regular_adjacency_list(std::move(data), 1)); + new (a) dolfinx::graph::AdjacencyList, + std::vector, U>( + dolfinx::graph::regular_adjacency_list(std::move(data), 1)); }, nb::arg("adj").noconvert()) .def( "__init__", - [](dolfinx::graph::AdjacencyList* a, + [](dolfinx::graph::AdjacencyList, + std::vector, U>* a, nb::ndarray, nb::c_contig> adj) { std::vector data(adj.data(), adj.data() + adj.size()); - new (a) dolfinx::graph::AdjacencyList( - dolfinx::graph::regular_adjacency_list(std::move(data), - adj.shape(1))); + new (a) dolfinx::graph::AdjacencyList, + std::vector, U>( + dolfinx::graph::regular_adjacency_list(std::move(data), + adj.shape(1))); }, nb::arg("adj").noconvert()) .def( "__init__", - [](dolfinx::graph::AdjacencyList* a, + [](dolfinx::graph::AdjacencyList, + std::vector, U>* a, nb::ndarray, nb::c_contig> array, nb::ndarray, nb::c_contig> displ) { std::vector data(array.data(), array.data() + array.size()); std::vector offsets(displ.data(), displ.data() + displ.size()); - new (a) dolfinx::graph::AdjacencyList(std::move(data), - std::move(offsets)); + new (a) dolfinx::graph::AdjacencyList, + std::vector, U>( + std::move(data), std::move(offsets)); }, nb::arg("data").noconvert(), nb::arg("offsets")) .def( "links", - [](const dolfinx::graph::AdjacencyList& self, int i) + [](const dolfinx::graph::AdjacencyList< + std::vector, std::vector, U>& self, + int i) { std::span link = self.links(i); return nb::ndarray(link.data(), {link.size()}); @@ -90,7 +100,8 @@ void declare_adjacency_list_init(nb::module_& m, std::string type) "Links (edges) of a node") .def_prop_ro( "array", - [](const dolfinx::graph::AdjacencyList& self) + [](const dolfinx::graph::AdjacencyList< + std::vector, std::vector, U>& self) { return nb::ndarray(self.array().data(), {self.array().size()}); @@ -98,37 +109,57 @@ void declare_adjacency_list_init(nb::module_& m, std::string type) nb::rv_policy::reference_internal) .def_prop_ro( "offsets", - [](const dolfinx::graph::AdjacencyList& self) + [](const dolfinx::graph::AdjacencyList< + std::vector, std::vector, U>& self) { return nb::ndarray( self.offsets().data(), {self.offsets().size()}); }, nb::rv_policy::reference_internal) - .def_prop_ro("num_nodes", &dolfinx::graph::AdjacencyList::num_nodes) - .def("__eq__", &dolfinx::graph::AdjacencyList::operator==, + .def_prop_ro("num_nodes", + &dolfinx::graph::AdjacencyList< + std::vector, std::vector, U>::num_nodes) + .def("__eq__", + &dolfinx::graph::AdjacencyList< + std::vector, std::vector, U>::operator==, nb::is_operator()) - .def("__repr__", &dolfinx::graph::AdjacencyList::str) - .def("__len__", &dolfinx::graph::AdjacencyList::num_nodes); + .def("__repr__", + &dolfinx::graph::AdjacencyList, + std::vector, U>::str) + .def("__len__", + &dolfinx::graph::AdjacencyList< + std::vector, std::vector, U>::num_nodes); } template void declare_adjacency_list(nb::module_& m, std::string type) { std::string pyclass_name = std::string("AdjacencyList_") + type; - nb::class_>(m, pyclass_name.c_str(), - "Adjacency List") + nb::class_, std::vector, std::vector>>( + m, pyclass_name.c_str(), "Adjacency List") .def_prop_ro( "offsets", - [](const dolfinx::graph::AdjacencyList& self) + [](const dolfinx::graph::AdjacencyList< + std::vector, std::vector, U>& self) { return nb::ndarray( self.offsets().data(), {self.offsets().size()}); }, nb::rv_policy::reference_internal) - .def_prop_ro("num_nodes", &dolfinx::graph::AdjacencyList::num_nodes) - .def("__eq__", &dolfinx::graph::AdjacencyList::operator==, + .def_prop_ro("num_nodes", + &dolfinx::graph::AdjacencyList, + std::vector, + std::vector>::num_nodes) + .def("__eq__", + &dolfinx::graph::AdjacencyList, + std::vector, + std::vector>::operator==, nb::is_operator()) - .def("__len__", &dolfinx::graph::AdjacencyList::num_nodes); + .def("__len__", + &dolfinx::graph::AdjacencyList, + std::vector, + std::vector>::num_nodes); } } // namespace @@ -143,9 +174,12 @@ void graph(nb::module_& m) m, "int_sizet_int8__int32_int32"); using partition_fn - = std::function( + = std::function, + std::vector>( MPICommWrapper, int, - const dolfinx::graph::AdjacencyList&, bool)>; + const dolfinx::graph::AdjacencyList, + std::vector>&, + bool)>; m.def( "partitioner", []() -> partition_fn { return create_partitioner_py(dolfinx::graph::partition_graph); }, @@ -198,12 +232,13 @@ void graph(nb::module_& m) m.def( "comm_graph_data", [](dolfinx::graph::AdjacencyList< - std::tuple, - std::pair>& g) + std::vector>, + std::vector, + std::vector>>& g) { std::vector>> adj; - for (std::int32_t n = 0; n < g.num_nodes(); ++n) + for (std::size_t n = 0; n < g.num_nodes(); ++n) { for (auto [e, w, local] : g.links(n)) { @@ -217,7 +252,7 @@ void graph(nb::module_& m) std::pair>> nodes; std::ranges::transform( - g.node_data().value(), std::ranges::views::iota(0), + g.node_data(), std::ranges::views::iota(0), std::back_inserter(nodes), [](auto data, auto n) { @@ -234,8 +269,9 @@ void graph(nb::module_& m) m.def( "comm_to_json", [](dolfinx::graph::AdjacencyList< - std::tuple, - std::pair>& g) + std::vector>, + std::vector, + std::vector>>& g) { return dolfinx::graph::comm_to_json(g); }, "Build a JSON string representation of a parallel communication " "graph that can use used by build a NetworkX graph."); diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index 4def5d96769..a52d1b8d9a4 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -114,7 +114,8 @@ void declare_meshtags(nb::module_& m, const std::string& type) m.def("create_meshtags", [](std::shared_ptr topology, int dim, - const dolfinx::graph::AdjacencyList& entities, + const dolfinx::graph::AdjacencyList>& + entities, nb::ndarray, nb::c_contig> values) { return dolfinx::mesh::create_meshtags( @@ -549,7 +550,7 @@ void mesh(nb::module_& m) m.def( "build_dual_graph", [](const MPICommWrapper comm, dolfinx::mesh::CellType cell_type, - const dolfinx::graph::AdjacencyList& cells, + const dolfinx::graph::AdjacencyList>& cells, std::optional max_facet_to_cell_links) { std::vector c = {cell_type}; @@ -642,7 +643,9 @@ void mesh(nb::module_& m) [](dolfinx::mesh::Topology* t, dolfinx::mesh::CellType cell_type, std::shared_ptr vertex_map, std::shared_ptr cell_map, - std::shared_ptr> cells, + std::shared_ptr< + dolfinx::graph::AdjacencyList>> + cells, std::optional< nb::ndarray, nb::c_contig>> original_index) @@ -805,10 +808,12 @@ void mesh(nb::module_& m) "Create default cell partitioner."); m.def( "create_cell_partitioner", - [](const std::function( - MPICommWrapper comm, int nparts, - const dolfinx::graph::AdjacencyList& local_graph, - bool ghosting)>& part, + [](const std::function< + dolfinx::graph::AdjacencyList>( + MPICommWrapper comm, int nparts, + const dolfinx::graph::AdjacencyList>& + local_graph, + bool ghosting)>& part, dolfinx::mesh::GhostMode mode, std::optional max_facet_to_cell_links) -> part::impl::PythonCellPartitionFunction