diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp index 65a8804950abad22fe2205569217cd8ff162b043..44c07b62877451940e975fbc929b6da1530a36f4 100644 --- a/src/mesh/Connectivity.cpp +++ b/src/mesh/Connectivity.cpp @@ -46,14 +46,6 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor) m_cell_number = WeakCellValue<const int>(*this, descriptor.cellNumberVector()); m_node_number = WeakNodeValue<const int>(*this, descriptor.nodeNumberVector()); - { - WeakCellValue<int> cell_global_index(*this); - int first_index = 0; - parallel_for( - this->numberOfCells(), PUGS_LAMBDA(CellId j) { cell_global_index[j] = first_index + j; }); - m_cell_global_index = cell_global_index; - } - m_cell_owner = WeakCellValue<const int>(*this, descriptor.cellOwnerVector()); { @@ -310,6 +302,80 @@ Connectivity<Dimension>::_write(std::ostream& os) const return os; } +template <size_t Dimension> +CRSGraph +Connectivity<Dimension>::ownCellToCellGraph() const +{ + std::vector<std::vector<int>> cell_cells(this->numberOfCells()); + const auto& face_to_cell_matrix = this->faceToCellMatrix(); + + for (FaceId l = 0; l < this->numberOfFaces(); ++l) { + const auto& face_to_cell = face_to_cell_matrix[l]; + if (face_to_cell.size() > 1) { + const CellId cell_0 = face_to_cell[0]; + const CellId cell_1 = face_to_cell[1]; + + if (m_cell_is_owned[cell_0]) { + cell_cells[cell_0].push_back(cell_1); + } + if (m_cell_is_owned[cell_1]) { + cell_cells[cell_1].push_back(cell_0); + } + } + } + + std::vector<std::vector<int>> owned_cell_cells; + owned_cell_cells.reserve(this->numberOfCells()); + + size_t number_of_owned_cells = 0; + for (CellId cell_id = 0; cell_id < this->numberOfCells(); ++cell_id) { + if (m_cell_is_owned[cell_id]) { + ++number_of_owned_cells; + } + } + + for (size_t i_cell = 0; i_cell < this->numberOfCells(); ++i_cell) { + if (cell_cells[i_cell].size() > 0) { + owned_cell_cells.emplace_back(std::move(cell_cells[i_cell])); + } + } + + Array number_of_owned_cells_by_rank = parallel::allGather(number_of_owned_cells); + + // index is defined locally. First we compute the first index on the + // processor + size_t global_index = 0; + for (size_t i = 0; i < parallel::rank(); ++i) { + global_index += number_of_owned_cells_by_rank[i]; + } + + CellValue<int> global_cell_index{*this}; + for (CellId cell_id = 0; cell_id < this->numberOfCells(); ++cell_id) { + if (m_cell_is_owned[cell_id]) { + global_cell_index[cell_id] = global_index++; + } + } + synchronize(global_cell_index); + + Array<int> entries(owned_cell_cells.size() + 1); + entries[0] = 0; + for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) { + entries[i_cell + 1] = entries[i_cell] + owned_cell_cells[i_cell].size(); + } + Array<int> neighbors(entries[owned_cell_cells.size()]); + { + size_t i = 0; + for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) { + for (CellId cell_id : owned_cell_cells[i_cell]) { + neighbors[i] = global_cell_index[cell_id]; + ++i; + } + } + } + + return CRSGraph(entries, neighbors); +} + template std::ostream& Connectivity<1>::_write(std::ostream&) const; template std::ostream& Connectivity<2>::_write(std::ostream&) const; template std::ostream& Connectivity<3>::_write(std::ostream&) const; @@ -333,3 +399,7 @@ template Connectivity<3>::Connectivity(); template void Connectivity<1>::_buildFrom(const ConnectivityDescriptor&); template void Connectivity<2>::_buildFrom(const ConnectivityDescriptor&); template void Connectivity<3>::_buildFrom(const ConnectivityDescriptor&); + +template CRSGraph Connectivity<1>::ownCellToCellGraph() const; +template CRSGraph Connectivity<2>::ownCellToCellGraph() const; +template CRSGraph Connectivity<3>::ownCellToCellGraph() const; diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp index faee4c28c34a879117c6418f2fe13e0f753a1e22..6d89a820fc25f2ed662f0e4f277a7f44fb136830 100644 --- a/src/mesh/Connectivity.hpp +++ b/src/mesh/Connectivity.hpp @@ -59,8 +59,6 @@ class Connectivity final : public IConnectivity ConnectivityMatrix m_item_to_item_matrix[Dimension + 1][Dimension + 1]; WeakCellValue<const CellType> m_cell_type; - WeakCellValue<const int> m_cell_global_index; - WeakCellValue<const int> m_cell_number; WeakFaceValue<const int> m_face_number; WeakEdgeValue<const int> m_edge_number; @@ -655,54 +653,7 @@ class Connectivity final : public IConnectivity } } - PUGS_INLINE - CRSGraph - cellToCellGraph() const - { - std::vector<std::vector<int>> cell_cells(this->numberOfCells()); - const auto& face_to_cell_matrix = this->faceToCellMatrix(); - - for (FaceId l = 0; l < this->numberOfFaces(); ++l) { - const auto& face_to_cell = face_to_cell_matrix[l]; - if (face_to_cell.size() > 1) { - const CellId cell_0 = face_to_cell[0]; - const CellId cell_1 = face_to_cell[1]; - - if (m_cell_is_owned[cell_0]) { - cell_cells[cell_0].push_back(cell_1); - } - if (m_cell_is_owned[cell_1]) { - cell_cells[cell_1].push_back(cell_0); - } - } - } - - std::vector<std::vector<int>> owned_cell_cells; - owned_cell_cells.reserve(this->numberOfCells()); - - for (size_t i_cell = 0; i_cell < this->numberOfCells(); ++i_cell) { - if (cell_cells[i_cell].size() > 0) { - owned_cell_cells.emplace_back(std::move(cell_cells[i_cell])); - } - } - - Array<int> entries(owned_cell_cells.size() + 1); - entries[0] = 0; - for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) { - entries[i_cell + 1] = entries[i_cell] + owned_cell_cells[i_cell].size(); - } - Array<int> neighbors(entries[owned_cell_cells.size()]); - { - size_t i = 0; - for (size_t i_cell = 0; i_cell < owned_cell_cells.size(); ++i_cell) { - for (CellId cell_id : owned_cell_cells[i_cell]) { - neighbors[i] = m_cell_global_index[cell_id]; - ++i; - } - } - } - return CRSGraph(entries, neighbors); - } + CRSGraph ownCellToCellGraph() const; PUGS_INLINE size_t diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp index 65133fda261190c008b7169c4d32090488379237..8c205a87183f698dee24b4335ea6c5356cc3e677 100644 --- a/src/mesh/ConnectivityDispatcher.cpp +++ b/src/mesh/ConnectivityDispatcher.cpp @@ -13,7 +13,7 @@ void ConnectivityDispatcher<Dimension>::_buildNewOwner() { if constexpr (item_type == ItemType::cell) { - CRSGraph connectivity_graph = m_connectivity.cellToCellGraph(); + CRSGraph connectivity_graph = m_connectivity.ownCellToCellGraph(); Partitioner P; Array new_owner_array = P.partition(connectivity_graph); CellValue cell_is_owned = m_connectivity.cellIsOwned();