diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 0ce2f7d76040..b974b13bd890 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -48,6 +48,61 @@ namespace Kratos ///@name Kratos Classes ///@{ +/** + * @brief This struct provides information related with distributed search + * @details The class contains the following: + * - PointCoordinates: The coordinates of the points to be created + * - Indexes: The indexes of the points to be created + * - Ranks: The ranks where the results will be defined + */ +struct DistributedSearchInformation +{ + /// Alias for the data type used to represent indices. + using IndexType = std::size_t; + + /// Alias for the data type used to represent sizes. + using SizeType = std::size_t; + + std::vector PointCoordinates; /// Vector to store point coordinates. + std::vector Indexes; /// Vector to store point indices. + std::vector SearchRanks; /// Ranks where search is invoked. + std::vector> Ranks; /// Vector of vectors to store ranks. + + /** + * @brief Reserve memory for the point data vectors. + * @param Size The size to reserve. + */ + void Reserve(const SizeType Size) + { + PointCoordinates.reserve(Size * 3); + Indexes.reserve(Size); + SearchRanks.reserve(Size); + Ranks.reserve(Size); + } + + /** + * @brief Shrink the capacity of the point data vectors to fit the data. + */ + void Shrink() + { + PointCoordinates.shrink_to_fit(); + Indexes.shrink_to_fit(); + SearchRanks.shrink_to_fit(); + Ranks.shrink_to_fit(); + } + + /** + * @brief Clear all the data in the point data vectors. + */ + void Clear() + { + PointCoordinates.clear(); + Indexes.clear(); + SearchRanks.clear(); + Ranks.clear(); + } +}; + /** * @class SearchUtilities * @ingroup KratosCore @@ -223,6 +278,41 @@ class SearchUtilities SynchronizePoints(itPointBegin, itPointEnd, rAllPointsCoordinates, rAllPointsIds, rDataCommunicator, number_of_points, total_number_of_points); } + /** + * @brief SynchronousPointSynchronization prepares synchronously the coordinates of the points for MPI search. + * @param itPointBegin Iterator to the beginning of the points range. + * @param itPointEnd Iterator to the end of the points range. + * @param rSearchInfo The class containing the result of the search. + * @param rBoundingBox The bounding box considered. + * @param ThresholdBoundingBox The threshold for computing is inside bounding box considered. + * @param rDataCommunicator The data communicator. + * @param IndexItIsJustCounter If the index considered it it just a counter. + * @tparam TPointIteratorType The type of the point iterator. + * @tparam TBoundingBoxType The type of the bounding box. + * @return The ids of all points. + */ + template + static std::vector SynchronousPointSynchronizationWithBoundingBox( + TPointIteratorType itPointBegin, + TPointIteratorType itPointEnd, + DistributedSearchInformation& rSearchInfo, + const TBoundingBoxType& rBoundingBox, + const double ThresholdBoundingBox, + const DataCommunicator& rDataCommunicator, + const bool IndexItIsJustCounter = false + ) + { + // First check that the points are the same in all processes + int number_of_points, total_number_of_points; + CalculateNumberOfPoints(itPointBegin, itPointEnd, number_of_points, total_number_of_points, rDataCommunicator); + + KRATOS_DEBUG_ERROR_IF(number_of_points < 0) << "The number of points is negative" << std::endl; + KRATOS_DEBUG_ERROR_IF(total_number_of_points < 0) << "The total number of points is negative" << std::endl; + + // We synchronize the points + return SynchronizePointsWithBoundingBox(itPointBegin, itPointEnd, rSearchInfo, rBoundingBox, ThresholdBoundingBox, rDataCommunicator, number_of_points, total_number_of_points, IndexItIsJustCounter); + } + /** * @brief SynchronousPointSynchronizationWithRecvSizes prepares synchronously the coordinates of the points for MPI search including the recv sizes * @details With recv sizes @@ -424,7 +514,7 @@ class SearchUtilities ///@} private: ///@name Private Operations - ///@{ + ///@{ /** * @details Synchronizes points between different processes. @@ -464,8 +554,12 @@ class SearchUtilities } // Initialize and resize vectors - rAllPointsCoordinates.resize(TotalNumberOfPoints * 3); - rAllPointsIds.resize(TotalNumberOfPoints); + if (rAllPointsCoordinates.size() != static_cast(TotalNumberOfPoints * 3)) { + rAllPointsCoordinates.resize(TotalNumberOfPoints * 3); + } + if (rAllPointsIds.size() != static_cast(TotalNumberOfPoints)) { + rAllPointsIds.resize(TotalNumberOfPoints); + } std::vector recv_sizes(world_size, 0); // Auxiliary variables @@ -543,6 +637,248 @@ class SearchUtilities return recv_sizes; } + /** + * @details Synchronizes points between different processes. + * @details Synchonously. + * @param itPointBegin Iterator pointing to the beginning of the range of points. + * @param itPointEnd Iterator pointing to the end of the range of points. + * @param rAllPointsCoordinates Vector to store the synchronized points' coordinates. + * @param rAllPointsIds The ids of all the points (just a counter for points, and ids for nodes). + * @param rDataCommunicator Object for data communication between processes. + * @param NumberOfPoints Local number of points to be synchronized. + * @param TotalNumberOfPoints Total number of points across all processes. + * @param IndexItIsJustCounter If the index considered it it just a counter. + * @return A vector containing the ranks of each point. + * @tparam TPointIteratorType The type of the point iterator. + */ + template + static std::vector SynchronizePointsWithRanks( + TPointIteratorType itPointBegin, + TPointIteratorType itPointEnd, + std::vector& rAllPointsCoordinates, + std::vector& rAllPointsIds, + const DataCommunicator& rDataCommunicator, + const int NumberOfPoints, + const int TotalNumberOfPoints, + const bool IndexItIsJustCounter = false + ) + { + // Define lambda to retrieve Id + const auto GetIdNode = [](std::vector& rIds, TPointIteratorType& ItPoint, const std::size_t Counter, const std::size_t InitialId) { + if constexpr (std::is_same::value || std::is_same::value) { + rIds[Counter] = ItPoint->Id(); + } else { + rIds[Counter] = InitialId + Counter; + } + }; + const auto GetIdJustCounter = [](std::vector& rIds, TPointIteratorType& ItPoint, const std::size_t Counter, const std::size_t InitialId) { + rIds[Counter] = Counter; + }; + const auto GetId = IndexItIsJustCounter ? GetIdJustCounter : GetIdNode; + + // Get the World Size and rank in MPI + const int world_size = rDataCommunicator.Size(); + const int rank = rDataCommunicator.Rank(); + + // Getting global number of points + std::vector points_per_partition(world_size); + std::vector send_points_per_partition(1, NumberOfPoints); + std::vector all_points_ranks(TotalNumberOfPoints); + rDataCommunicator.AllGather(send_points_per_partition, points_per_partition); + std::size_t initial_id = 0; + if constexpr (!std::is_same::value || std::is_same::value) { + initial_id = std::accumulate(points_per_partition.begin(), points_per_partition.begin() + rank, 0); + } + + // Initialize and resize vectors + if (rAllPointsCoordinates.size() != static_cast(TotalNumberOfPoints * 3)) { + rAllPointsCoordinates.resize(TotalNumberOfPoints * 3); + } + if (rAllPointsIds.size() != static_cast(TotalNumberOfPoints)) { + rAllPointsIds.resize(TotalNumberOfPoints); + } + std::vector recv_sizes(world_size, 0); + + // Auxiliary variables + std::size_t counter = 0; + unsigned int i_coord; + + // If distributed + if (rDataCommunicator.IsDistributed()) { + // Initialize local points coordinates + std::vector send_points_coordinates(NumberOfPoints * 3); + std::vector send_points_ids(NumberOfPoints); + std::vector send_points_ranks(NumberOfPoints); + for (auto it_point = itPointBegin ; it_point != itPointEnd ; ++it_point) { + // In case of considering nodes + if constexpr (std::is_same::value || std::is_same::value) { + if (it_point->FastGetSolutionStepValue(PARTITION_INDEX) != rank) { + continue; // Skip if not local + } + } + const auto& r_coordinates = it_point->Coordinates(); + GetId(send_points_ids, it_point, counter, initial_id); + send_points_ranks[counter] = rank; + for (i_coord = 0; i_coord < 3; ++i_coord) { + send_points_coordinates[3 * counter + i_coord] = r_coordinates[i_coord]; + } + ++counter; + } + + /* Sync coordinates */ + + // Generate vectors with sizes for AllGatherv + for (int i_rank = 0; i_rank < world_size; ++i_rank) { + recv_sizes[i_rank] = 3 * points_per_partition[i_rank]; + } + std::vector recv_offsets(world_size, 0); + for (int i_rank = 1; i_rank < world_size; ++i_rank) { + recv_offsets[i_rank] = recv_offsets[i_rank - 1] + recv_sizes[i_rank - 1]; + } + + // Invoke AllGatherv + rDataCommunicator.AllGatherv(send_points_coordinates, rAllPointsCoordinates, recv_sizes, recv_offsets); + + /* Sync Ids */ + + // Generate vectors with sizes for AllGatherv + for (int i_rank = 0; i_rank < world_size; ++i_rank) { + recv_sizes[i_rank] = points_per_partition[i_rank]; + } + for (int i_rank = 1; i_rank < world_size; ++i_rank) { + recv_offsets[i_rank] = recv_offsets[i_rank - 1] + recv_sizes[i_rank - 1]; + } + + // Invoke AllGatherv + rDataCommunicator.AllGatherv(send_points_ids, rAllPointsIds, recv_sizes, recv_offsets); + rDataCommunicator.AllGatherv(send_points_ranks, all_points_ranks, recv_sizes, recv_offsets); + } else { // Serial + // Assign values + for (auto it_point = itPointBegin ; it_point != itPointEnd ; ++it_point) { + const auto& r_coordinates = it_point->Coordinates(); + if constexpr (std::is_same::value || std::is_same::value) { + rAllPointsIds[counter] = it_point->Id(); + } else { + rAllPointsIds[counter] = initial_id + counter; + } + for (i_coord = 0; i_coord < 3; ++i_coord) { + rAllPointsCoordinates[3 * counter + i_coord] = r_coordinates[i_coord]; + } + ++counter; + } + } + + // Return the ranks + return all_points_ranks; + } + + /** + * @brief Synchronizes points between different processes. + * @details Synchonously. + * @param itPointBegin Iterator pointing to the beginning of the range of points. + * @param itPointEnd Iterator pointing to the end of the range of points. + * @param rSearchInfo The class containing the result of the search. + * @param rBoundingBox The bounding box considered. + * @param ThresholdBoundingBox The threshold for computing is inside bounding box considered. + * @param rDataCommunicator Object for data communication between processes. + * @param NumberOfPoints Local number of points to be synchronized. + * @param TotalNumberOfPoints Total number of points across all processes. + * @param IndexItIsJustCounter If the index considered it it just a counter. + * @return A vector containing the sizes of data for each process. + * @tparam TPointIteratorType The type of the point iterator. + * @tparam TBoundingBoxType The type of the bounding box. + * @return The ids of all points. + */ + template + static std::vector SynchronizePointsWithBoundingBox( + TPointIteratorType itPointBegin, + TPointIteratorType itPointEnd, + DistributedSearchInformation& rSearchInfo, + const TBoundingBoxType& rBoundingBox, + const double ThresholdBoundingBox, + const DataCommunicator& rDataCommunicator, + const int NumberOfPoints, + const int TotalNumberOfPoints, + const bool IndexItIsJustCounter = false + ) + { + // Initialize and resize vectors + rSearchInfo.Reserve(TotalNumberOfPoints); + std::vector all_points_coordinates(TotalNumberOfPoints * 3); + std::vector all_points_ids(TotalNumberOfPoints); + + // Sync all points first + std::vector all_points_ranks = SynchronizePointsWithRanks(itPointBegin, itPointEnd, all_points_coordinates, all_points_ids, rDataCommunicator, NumberOfPoints, TotalNumberOfPoints, IndexItIsJustCounter); + + // Some definitions + IndexType i_coord = 0; + + // If distributed + if (rDataCommunicator.IsDistributed()) { + // Prepare MPI data + const int rank = rDataCommunicator.Rank(); + const int world_size = rDataCommunicator.Size(); + std::vector ranks(1, rank); + std::vector inside_ranks(world_size); + + // Fill actual vectors + Point point_to_test; + for (IndexType i_point = 0; i_point < static_cast(TotalNumberOfPoints); ++i_point) { + point_to_test[0] = all_points_coordinates[3 * i_point + 0]; + point_to_test[1] = all_points_coordinates[3 * i_point + 1]; + point_to_test[2] = all_points_coordinates[3 * i_point + 2]; + const bool is_inside = PointIsInsideBoundingBox(rBoundingBox, point_to_test, ThresholdBoundingBox); + const int search_rank = all_points_ranks[i_point]; + const bool to_be_included = is_inside || search_rank == rank; + inside_ranks.resize(world_size); + if (to_be_included) { + for (i_coord = 0; i_coord < 3; ++i_coord) { + rSearchInfo.PointCoordinates.push_back(all_points_coordinates[3 * i_point + i_coord]); + } + rSearchInfo.Indexes.push_back(all_points_ids[i_point]); + ranks[0] = rank; + } else { + ranks[0] = -1; + } + rDataCommunicator.AllGather(ranks, inside_ranks); + // Use std::remove_if and vector::erase to remove elements less than 0 + inside_ranks.erase( + std::remove_if( + inside_ranks.begin(), + inside_ranks.end(), + [](const int rank) { return rank < 0; } + ), + inside_ranks.end() + ); + + // Now adding // NOTE: Should be ordered, so a priori is not required to reorder + if (to_be_included) { + rSearchInfo.Ranks.push_back(inside_ranks); + rSearchInfo.SearchRanks.push_back(search_rank); + } + } + } else { // Serial + // Assign values + const std::size_t total_number_of_points = itPointEnd - itPointBegin; + for (IndexType i_point = 0 ; i_point < total_number_of_points; ++i_point) { + auto it_point = itPointBegin + i_point; + if (PointIsInsideBoundingBox(rBoundingBox, *it_point, ThresholdBoundingBox)) { + const auto& r_coordinates = it_point->Coordinates(); + for (i_coord = 0; i_coord < 3; ++i_coord) { + rSearchInfo.PointCoordinates.push_back(r_coordinates[i_coord]); + } + rSearchInfo.Indexes.push_back(all_points_ids[i_point]); + rSearchInfo.Ranks.push_back({0}); + } + } + } + + // Shrink to actual size + rSearchInfo.Shrink(); + + return all_points_ids; + } + /** * @brief Synchronizes the radius of all points in a distributed system. * @param rRecvSizes a vector containing the number of points to be received from each rank