diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp index 67463bb8a043424bc173073cff3413ce884f5600..aa8ca779cfb965ffb3e85e4feb1b4b9feabb79ca 100644 --- a/src/mesh/ConnectivityComputer.cpp +++ b/src/mesh/ConnectivityComputer.cpp @@ -4,8 +4,6 @@ #include <utils/Exceptions.hpp> #include <iostream> -#include <map> -#include <sstream> template <typename ConnectivityType> PUGS_INLINE ConnectivityMatrix @@ -17,8 +15,6 @@ ConnectivityComputer::computeConnectivityMatrix(const ConnectivityType& connecti if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) { const ConnectivityMatrix& child_to_item_matrix = connectivity.getMatrix(child_item_type, item_type); - std::cout << "computing connectivity " << itemName(item_type) << " -> " << itemName(child_item_type) << '\n'; - item_to_child_item_matrix = this->_computeInverse(child_to_item_matrix); } else { std::stringstream error_msg; @@ -44,32 +40,46 @@ template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(cons ConnectivityMatrix ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_matrix) const { - std::map<unsigned int, std::vector<unsigned int>> child_to_item_vector_map; - const size_t& number_of_items = item_to_child_matrix.numberOfRows(); + const size_t& number_of_rows = item_to_child_matrix.numberOfRows(); + + if ((item_to_child_matrix.values().size() > 0)) { + const size_t& number_of_columns = max(item_to_child_matrix.values()); + + Array<uint32_t> transposed_next_free_column_index(number_of_columns + 1); + transposed_next_free_column_index.fill(0); - for (unsigned int j = 0; j < number_of_items; ++j) { - const auto& item_to_child = item_to_child_matrix[j]; - for (unsigned int r = 0; r < item_to_child.size(); ++r) { - child_to_item_vector_map[item_to_child[r]].push_back(j); + Array<uint32_t> transposed_rows_map(transposed_next_free_column_index.size() + 1); + transposed_rows_map.fill(0); + for (size_t i = 0; i < number_of_rows; ++i) { + for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) { + transposed_rows_map[item_to_child_matrix.values()[j] + 1]++; + } } - } + for (size_t i = 1; i < transposed_rows_map.size(); ++i) { + transposed_rows_map[i] += transposed_rows_map[i - 1]; + } + Array<uint32_t> transposed_column_indices(transposed_rows_map[transposed_rows_map.size() - 1]); - { - size_t i = 0; - for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) { - if (child_item_id != i) { - throw NotImplementedError("sparse item numerotation"); + for (size_t i = 0; i < number_of_rows; ++i) { + for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) { + size_t i_column_index = item_to_child_matrix.values()[j]; + uint32_t& shift = transposed_next_free_column_index[i_column_index]; + + transposed_column_indices[transposed_rows_map[i_column_index] + shift] = i; + Assert(shift < transposed_rows_map[i_column_index + 1]); + shift++; } - ++i; } - } - std::vector<std::vector<unsigned int>> child_to_items_vector(child_to_item_vector_map.size()); - for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) { - child_to_items_vector[child_item_id] = item_vector; - } + return ConnectivityMatrix{transposed_rows_map, transposed_column_indices}; + } else { + // empty connectivity + Array<uint32_t> transposed_rows_map{1}; + transposed_rows_map[0] = 0; + Array<uint32_t> transposed_column_indices; - return ConnectivityMatrix(child_to_items_vector); + return ConnectivityMatrix{transposed_rows_map, transposed_column_indices}; + } } template <typename ItemOfItem, typename ConnectivityType> diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp index e51ba189221faeef80fb7966310d9e5c82527d15..6b777288cd6b8e3793a9addff6629c661d45d82d 100644 --- a/src/mesh/ConnectivityMatrix.hpp +++ b/src/mesh/ConnectivityMatrix.hpp @@ -56,7 +56,22 @@ class ConnectivityMatrix } PUGS_INLINE - ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer) : m_is_built{true} + ConnectivityMatrix(const Array<const uint32_t>& row_map, + const Array<const uint32_t>& column_indices) noexcept(NO_ASSERT) + : m_row_map{row_map}, m_column_indices{column_indices}, m_is_built{true} + { + Assert(m_column_indices.size() == m_row_map[m_row_map.size() - 1], "incompatible row map and column indices"); + Assert(m_row_map.size() > 0, "invalid row map"); + Assert(m_row_map[0] == 0, "row map should start with 0"); +#ifndef NDEBUG + for (size_t i = 1; i < m_row_map.size(); ++i) { + Assert(m_row_map[i] > m_row_map[i - 1], "row map values must be strictly increasing"); + } +#endif // NDEBUG + } + + PUGS_INLINE + ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer) noexcept : m_is_built{true} { m_row_map = [&] { Array<uint32_t> row_map(initializer.size() + 1);