From c35879c72df4872e31ad8008276b31948c994c01 Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Thu, 23 Nov 2023 18:30:15 +0100 Subject: [PATCH 01/13] Adding `SynchronizePointsWithBoundingBox` --- kratos/utilities/search_utilities.h | 120 +++++++++++++++++++++++++++- 1 file changed, 117 insertions(+), 3 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 0ce2f7d76040..0dff7e3c1736 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -223,6 +223,38 @@ 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 rAllPointsCoordinates vector where the computed coordinates will be stored + * @param rAllPointsIds The ids of all the points (just a counter for points, and ids for nodes) + * @param rBoundingBox The bounding box considered + * @param rDataCommunicator The data communicator + * @tparam TPointIteratorType The type of the point iterator + * @tparam TBoundingBoxType The type of the bounding box + */ + template + static std::vector SynchronousPointSynchronizationWithBoundingBox( + TPointIteratorType itPointBegin, + TPointIteratorType itPointEnd, + std::vector& rAllPointsCoordinates, + std::vector& rAllPointsIds, + const TBoundingBoxType& rBoundingBox, + const DataCommunicator& rDataCommunicator + ) + { + // 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, rAllPointsCoordinates, rAllPointsIds, rBoundingBox, rDataCommunicator, number_of_points, total_number_of_points); + } + /** * @brief SynchronousPointSynchronizationWithRecvSizes prepares synchronously the coordinates of the points for MPI search including the recv sizes * @details With recv sizes @@ -424,7 +456,7 @@ class SearchUtilities ///@} private: ///@name Private Operations - ///@{ + ///@{ /** * @details Synchronizes points between different processes. @@ -464,8 +496,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 +579,84 @@ 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 rBoundingBox The 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 + * @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 + */ + template + static std::vector SynchronizePointsWithBoundingBox( + TPointIteratorType itPointBegin, + TPointIteratorType itPointEnd, + std::vector& rAllPointsCoordinates, + std::vector& rAllPointsIds, + const TBoundingBoxType& rBoundingBox, + const DataCommunicator& rDataCommunicator, + const int NumberOfPoints, + const int TotalNumberOfPoints + ) + { + // Initialize and resize vectors + rAllPointsCoordinates.reserve(TotalNumberOfPoints * 3); + std::vector all_points_coordinates(TotalNumberOfPoints * 3); + rAllPointsIds.reserve(TotalNumberOfPoints); + std::vector all_points_ids(TotalNumberOfPoints); + + // Sync all points first + SynchronizePoints(itPointBegin, itPointEnd, all_points_coordinates, all_points_ids, rDataCommunicator, NumberOfPoints, TotalNumberOfPoints); + + // Some definitions + IndexType i_coord = 0; + array_1d coordinates; + + // If distributed + if (rDataCommunicator.IsDistributed()) { + // 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]; + if (PointIsInsideBoundingBox(rBoundingBox, point_to_test)) { + for (i_coord = 0; i_coord < 3; ++i_coord) { + rAllPointsCoordinates.push_back(all_points_coordinates[3 * i_point + i_coord]); + } + rAllPointsIds.push_back(all_points_ids[i_point]); + } + } + } 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)) { + noalias(coordinates) = it_point->Coordinates(); + rAllPointsIds.push_back(all_points_ids[i_point]); + for (i_coord = 0; i_coord < 3; ++i_coord) { + rAllPointsCoordinates.push_back(coordinates[i_coord]); + } + } + } + } + + // Shrink to actual size + rAllPointsIds.shrink_to_fit(); + rAllPointsCoordinates.shrink_to_fit(); + + // Return the recv sizes + 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 From e443fcfbbe0e0df92bc192b4d0eb2ed2339639cc Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Thu, 23 Nov 2023 19:34:22 +0100 Subject: [PATCH 02/13] Update --- kratos/utilities/search_utilities.h | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 0dff7e3c1736..87ccb380b685 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -230,6 +230,7 @@ class SearchUtilities * @param rAllPointsCoordinates vector where the computed coordinates will be stored * @param rAllPointsIds The ids of all the points (just a counter for points, and ids for nodes) * @param rBoundingBox The bounding box considered + * @param ThresholdBoundingBox The threshold for computing is inside bounding box considered * @param rDataCommunicator The data communicator * @tparam TPointIteratorType The type of the point iterator * @tparam TBoundingBoxType The type of the bounding box @@ -241,6 +242,7 @@ class SearchUtilities std::vector& rAllPointsCoordinates, std::vector& rAllPointsIds, const TBoundingBoxType& rBoundingBox, + const double ThresholdBoundingBox, const DataCommunicator& rDataCommunicator ) { @@ -252,7 +254,7 @@ class SearchUtilities 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, rAllPointsCoordinates, rAllPointsIds, rBoundingBox, rDataCommunicator, number_of_points, total_number_of_points); + return SynchronizePointsWithBoundingBox(itPointBegin, itPointEnd, rAllPointsCoordinates, rAllPointsIds, rBoundingBox, ThresholdBoundingBox, rDataCommunicator, number_of_points, total_number_of_points); } /** @@ -587,6 +589,7 @@ class SearchUtilities * @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 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 @@ -601,6 +604,7 @@ class SearchUtilities std::vector& rAllPointsCoordinates, std::vector& rAllPointsIds, const TBoundingBoxType& rBoundingBox, + const double ThresholdBoundingBox, const DataCommunicator& rDataCommunicator, const int NumberOfPoints, const int TotalNumberOfPoints @@ -627,7 +631,7 @@ class SearchUtilities 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]; - if (PointIsInsideBoundingBox(rBoundingBox, point_to_test)) { + if (PointIsInsideBoundingBox(rBoundingBox, point_to_test, ThresholdBoundingBox)) { for (i_coord = 0; i_coord < 3; ++i_coord) { rAllPointsCoordinates.push_back(all_points_coordinates[3 * i_point + i_coord]); } @@ -639,7 +643,7 @@ class SearchUtilities 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)) { + if (PointIsInsideBoundingBox(rBoundingBox, *it_point, ThresholdBoundingBox)) { noalias(coordinates) = it_point->Coordinates(); rAllPointsIds.push_back(all_points_ids[i_point]); for (i_coord = 0; i_coord < 3; ++i_coord) { From 1ecd720ff4c908c52bf61dea3be914bd4bcba537 Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Tue, 5 Dec 2023 15:04:21 +0100 Subject: [PATCH 03/13] [Core] Refactor search utilities in order to filter by partition --- kratos/utilities/search_utilities.h | 119 ++++++++++++++++++++++------ 1 file changed, 95 insertions(+), 24 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 87ccb380b685..adff24f13f8a 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -48,6 +48,57 @@ 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. // NOTE: This could be removed in the future + 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); + 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(); + Ranks.shrink_to_fit(); + } + + /** + * @brief Clear all the data in the point data vectors. + */ + void Clear() + { + PointCoordinates.clear(); + Indexes.clear(); + Ranks.clear(); + } +}; + /** * @class SearchUtilities * @ingroup KratosCore @@ -227,8 +278,7 @@ class SearchUtilities * @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 rAllPointsCoordinates vector where the computed coordinates will be stored - * @param rAllPointsIds The ids of all the points (just a counter for points, and ids for nodes) + * @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 @@ -236,11 +286,10 @@ class SearchUtilities * @tparam TBoundingBoxType The type of the bounding box */ template - static std::vector SynchronousPointSynchronizationWithBoundingBox( + static void SynchronousPointSynchronizationWithBoundingBox( TPointIteratorType itPointBegin, TPointIteratorType itPointEnd, - std::vector& rAllPointsCoordinates, - std::vector& rAllPointsIds, + DistributedSearchInformation& rSearchInfo, const TBoundingBoxType& rBoundingBox, const double ThresholdBoundingBox, const DataCommunicator& rDataCommunicator @@ -254,7 +303,7 @@ class SearchUtilities 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, rAllPointsCoordinates, rAllPointsIds, rBoundingBox, ThresholdBoundingBox, rDataCommunicator, number_of_points, total_number_of_points); + SynchronizePointsWithBoundingBox(itPointBegin, itPointEnd, rSearchInfo, rBoundingBox, ThresholdBoundingBox, rDataCommunicator, number_of_points, total_number_of_points); } /** @@ -582,12 +631,11 @@ class SearchUtilities } /** - * @details Synchronizes points between different processes. + * @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 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 @@ -596,13 +644,13 @@ class SearchUtilities * @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 distributed search information */ template - static std::vector SynchronizePointsWithBoundingBox( + static void SynchronizePointsWithBoundingBox( TPointIteratorType itPointBegin, TPointIteratorType itPointEnd, - std::vector& rAllPointsCoordinates, - std::vector& rAllPointsIds, + DistributedSearchInformation& rSearchInfo, const TBoundingBoxType& rBoundingBox, const double ThresholdBoundingBox, const DataCommunicator& rDataCommunicator, @@ -611,9 +659,8 @@ class SearchUtilities ) { // Initialize and resize vectors - rAllPointsCoordinates.reserve(TotalNumberOfPoints * 3); + rSearchInfo.Reserve(TotalNumberOfPoints); std::vector all_points_coordinates(TotalNumberOfPoints * 3); - rAllPointsIds.reserve(TotalNumberOfPoints); std::vector all_points_ids(TotalNumberOfPoints); // Sync all points first @@ -625,17 +672,43 @@ class SearchUtilities // 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]; - if (PointIsInsideBoundingBox(rBoundingBox, point_to_test, ThresholdBoundingBox)) { + const bool is_inside = PointIsInsideBoundingBox(rBoundingBox, point_to_test, ThresholdBoundingBox); + inside_ranks.resize(world_size); + if (is_inside) { for (i_coord = 0; i_coord < 3; ++i_coord) { - rAllPointsCoordinates.push_back(all_points_coordinates[3 * i_point + i_coord]); + rSearchInfo.PointCoordinates.push_back(all_points_coordinates[3 * i_point + i_coord]); } - rAllPointsIds.push_back(all_points_ids[i_point]); + 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 (is_inside) { + rSearchInfo.Ranks.push_back(inside_ranks); } } } else { // Serial @@ -645,20 +718,18 @@ class SearchUtilities auto it_point = itPointBegin + i_point; if (PointIsInsideBoundingBox(rBoundingBox, *it_point, ThresholdBoundingBox)) { noalias(coordinates) = it_point->Coordinates(); - rAllPointsIds.push_back(all_points_ids[i_point]); for (i_coord = 0; i_coord < 3; ++i_coord) { - rAllPointsCoordinates.push_back(coordinates[i_coord]); + rSearchInfo.PointCoordinates.push_back(coordinates[i_coord]); } + rSearchInfo.Indexes.push_back(all_points_ids[i_point]); + std::vector ranks(1, 0); + rSearchInfo.Ranks.push_back(ranks); } } } // Shrink to actual size - rAllPointsIds.shrink_to_fit(); - rAllPointsCoordinates.shrink_to_fit(); - - // Return the recv sizes - return all_points_ids; + rSearchInfo.Shrink(); } /** From 0175fbad869c7a337cb6e391c3ee565b01f63e4f Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Wed, 6 Dec 2023 11:08:42 +0100 Subject: [PATCH 04/13] Missing return --- kratos/utilities/search_utilities.h | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index adff24f13f8a..ed4732ae488f 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -284,9 +284,10 @@ class SearchUtilities * @param rDataCommunicator The data communicator * @tparam TPointIteratorType The type of the point iterator * @tparam TBoundingBoxType The type of the bounding box + * @return The ids of all points */ template - static void SynchronousPointSynchronizationWithBoundingBox( + static std::vector SynchronousPointSynchronizationWithBoundingBox( TPointIteratorType itPointBegin, TPointIteratorType itPointEnd, DistributedSearchInformation& rSearchInfo, @@ -303,7 +304,7 @@ class SearchUtilities KRATOS_DEBUG_ERROR_IF(total_number_of_points < 0) << "The total number of points is negative" << std::endl; // We synchronize the points - SynchronizePointsWithBoundingBox(itPointBegin, itPointEnd, rSearchInfo, rBoundingBox, ThresholdBoundingBox, rDataCommunicator, number_of_points, total_number_of_points); + return SynchronizePointsWithBoundingBox(itPointBegin, itPointEnd, rSearchInfo, rBoundingBox, ThresholdBoundingBox, rDataCommunicator, number_of_points, total_number_of_points); } /** @@ -644,10 +645,10 @@ class SearchUtilities * @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 distributed search information + * @return The ids of all points */ template - static void SynchronizePointsWithBoundingBox( + static std::vector SynchronizePointsWithBoundingBox( TPointIteratorType itPointBegin, TPointIteratorType itPointEnd, DistributedSearchInformation& rSearchInfo, @@ -730,6 +731,8 @@ class SearchUtilities // Shrink to actual size rSearchInfo.Shrink(); + + return all_points_ids; } /** From 7283a9dd748aedf93b9ac022d25f36e911fd44a9 Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Thu, 7 Dec 2023 14:52:01 +0100 Subject: [PATCH 05/13] Fix empty cases --- kratos/utilities/search_utilities.h | 132 +++++++++++++++++++++++++++- 1 file changed, 129 insertions(+), 3 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index ed4732ae488f..99b3b875c6f6 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -631,6 +631,131 @@ 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 + * @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 + ) + { + // 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; + array_1d coordinates; + 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 + } + } + noalias(coordinates) = it_point->Coordinates(); + if constexpr (std::is_same::value || std::is_same::value) { + send_points_ids[counter] = it_point->Id(); + } else { + send_points_ids[counter] = initial_id + counter; + } + send_points_ranks[counter] = rank; + for (i_coord = 0; i_coord < 3; ++i_coord) { + send_points_coordinates[3 * counter + i_coord] = 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) { + noalias(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] = coordinates[i_coord]; + } + ++counter; + } + } + + // Return the ranks + return all_points_ranks; + } + /** * @details Synchronizes points between different processes. * @details Synchonously @@ -665,7 +790,7 @@ class SearchUtilities std::vector all_points_ids(TotalNumberOfPoints); // Sync all points first - SynchronizePoints(itPointBegin, itPointEnd, all_points_coordinates, all_points_ids, rDataCommunicator, NumberOfPoints, TotalNumberOfPoints); + std::vector all_points_ranks = SynchronizePointsWithRanks(itPointBegin, itPointEnd, all_points_coordinates, all_points_ids, rDataCommunicator, NumberOfPoints, TotalNumberOfPoints); // Some definitions IndexType i_coord = 0; @@ -686,8 +811,9 @@ class SearchUtilities 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 bool to_be_included = is_inside || all_points_ranks[i_point] == rank; inside_ranks.resize(world_size); - if (is_inside) { + 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]); } @@ -708,7 +834,7 @@ class SearchUtilities ); // Now adding // NOTE: Should be ordered, so a priori is not required to reorder - if (is_inside) { + if (to_be_included) { rSearchInfo.Ranks.push_back(inside_ranks); } } From 867183fb0166fb442bd52d7b0326395e962b4563 Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Fri, 15 Dec 2023 10:51:15 +0100 Subject: [PATCH 06/13] Adding more info --- kratos/utilities/search_utilities.h | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 99b3b875c6f6..a636605f2d6f 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -64,7 +64,8 @@ struct DistributedSearchInformation using SizeType = std::size_t; std::vector PointCoordinates; /// Vector to store point coordinates. - std::vector Indexes; /// Vector to store point indices. // NOTE: This could be removed in the future + 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. /** @@ -75,6 +76,7 @@ struct DistributedSearchInformation { PointCoordinates.reserve(Size * 3); Indexes.reserve(Size); + SearchRanks.reserve(Size); Ranks.reserve(Size); } @@ -85,6 +87,7 @@ struct DistributedSearchInformation { PointCoordinates.shrink_to_fit(); Indexes.shrink_to_fit(); + SearchRanks.shrink_to_fit(); Ranks.shrink_to_fit(); } @@ -95,6 +98,7 @@ struct DistributedSearchInformation { PointCoordinates.clear(); Indexes.clear(); + SearchRanks.clear(); Ranks.clear(); } }; @@ -811,7 +815,8 @@ class SearchUtilities 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 bool to_be_included = is_inside || all_points_ranks[i_point] == rank; + 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) { @@ -836,6 +841,7 @@ class SearchUtilities // 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 From 1a8b3605855ea80a5131a425441dfa339f9bcbc5 Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Fri, 15 Dec 2023 16:16:01 +0100 Subject: [PATCH 07/13] More generic --- kratos/utilities/search_utilities.h | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index a636605f2d6f..3568b00fd144 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -656,9 +656,23 @@ class SearchUtilities std::vector& rAllPointsIds, const DataCommunicator& rDataCommunicator, const int NumberOfPoints, - const int TotalNumberOfPoints + 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(); @@ -701,11 +715,7 @@ class SearchUtilities } } noalias(coordinates) = it_point->Coordinates(); - if constexpr (std::is_same::value || std::is_same::value) { - send_points_ids[counter] = it_point->Id(); - } else { - send_points_ids[counter] = initial_id + counter; - } + 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] = coordinates[i_coord]; From e55357326acadfa5b966a3bf8c0752c8af46ebd2 Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Fri, 15 Dec 2023 16:21:09 +0100 Subject: [PATCH 08/13] Update --- kratos/utilities/search_utilities.h | 57 +++++++++++++++-------------- 1 file changed, 30 insertions(+), 27 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 3568b00fd144..ed7ebce18853 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -636,17 +636,18 @@ class SearchUtilities } /** - * @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 - * @return A vector containing the ranks of each point - * @tparam TPointIteratorType The type of the point iterator + * @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( @@ -771,20 +772,21 @@ class SearchUtilities } /** - * @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 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 - * @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 + * @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( @@ -795,7 +797,8 @@ class SearchUtilities const double ThresholdBoundingBox, const DataCommunicator& rDataCommunicator, const int NumberOfPoints, - const int TotalNumberOfPoints + const int TotalNumberOfPoints, + const bool IndexItIsJustCounter = false ) { // Initialize and resize vectors @@ -804,7 +807,7 @@ class SearchUtilities 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); + std::vector all_points_ranks = SynchronizePointsWithRanks(itPointBegin, itPointEnd, all_points_coordinates, all_points_ids, rDataCommunicator, NumberOfPoints, TotalNumberOfPoints, IndexItIsJustCounter); // Some definitions IndexType i_coord = 0; From 3e9a095a4b3bc407422fa9f25dabf7387a389c27 Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Fri, 15 Dec 2023 16:53:07 +0100 Subject: [PATCH 09/13] Missing --- kratos/utilities/search_utilities.h | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index ed7ebce18853..2afc3cfa9632 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -280,15 +280,16 @@ class SearchUtilities /** * @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 - * @tparam TPointIteratorType The type of the point iterator - * @tparam TBoundingBoxType The type of the bounding box - * @return The ids of all points + * @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( @@ -297,7 +298,8 @@ class SearchUtilities DistributedSearchInformation& rSearchInfo, const TBoundingBoxType& rBoundingBox, const double ThresholdBoundingBox, - const DataCommunicator& rDataCommunicator + const DataCommunicator& rDataCommunicator, + const bool IndexItIsJustCounter = false ) { // First check that the points are the same in all processes @@ -308,7 +310,7 @@ class SearchUtilities 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); + return SynchronizePointsWithBoundingBox(itPointBegin, itPointEnd, rSearchInfo, rBoundingBox, ThresholdBoundingBox, rDataCommunicator, number_of_points, total_number_of_points, IndexItIsJustCounter); } /** From 71bfa0cd1006b9a43b947d0e5bbe2079782a0719 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Wed, 20 Dec 2023 15:18:37 +0100 Subject: [PATCH 10/13] Update kratos/utilities/search_utilities.h Co-authored-by: Carlos Roig --- kratos/utilities/search_utilities.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 2afc3cfa9632..89cf74ee3682 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -691,7 +691,7 @@ class SearchUtilities } // Initialize and resize vectors - if (rAllPointsCoordinates.size() != static_cast(TotalNumberOfPoints * 3)) { + if (rAllPointsCoordinates.size() != static_cast(TotalNumberOfPoints * 3)) { rAllPointsCoordinates.resize(TotalNumberOfPoints * 3); } if (rAllPointsIds.size() != static_cast(TotalNumberOfPoints)) { From d63bfa87768d8a3c805eef79b0029ccc355495f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Wed, 20 Dec 2023 15:18:43 +0100 Subject: [PATCH 11/13] Update kratos/utilities/search_utilities.h Co-authored-by: Carlos Roig --- kratos/utilities/search_utilities.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 89cf74ee3682..fe0152544c4b 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -694,7 +694,7 @@ class SearchUtilities if (rAllPointsCoordinates.size() != static_cast(TotalNumberOfPoints * 3)) { rAllPointsCoordinates.resize(TotalNumberOfPoints * 3); } - if (rAllPointsIds.size() != static_cast(TotalNumberOfPoints)) { + if (rAllPointsIds.size() != static_cast(TotalNumberOfPoints)) { rAllPointsIds.resize(TotalNumberOfPoints); } std::vector recv_sizes(world_size, 0); From e2884c3fc39bb872ca7bcd433b5f310364dea98e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vicente=20Mataix=20Ferr=C3=A1ndiz?= Date: Wed, 20 Dec 2023 16:12:50 +0100 Subject: [PATCH 12/13] Recommendation Co-authored-by: Carlos Roig --- kratos/utilities/search_utilities.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index fe0152544c4b..2bf740f06d93 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -870,8 +870,7 @@ class SearchUtilities rSearchInfo.PointCoordinates.push_back(coordinates[i_coord]); } rSearchInfo.Indexes.push_back(all_points_ids[i_point]); - std::vector ranks(1, 0); - rSearchInfo.Ranks.push_back(ranks); + rSearchInfo.Ranks.push_back({0}); } } } From b8faa2596bef7d41207342806a0deaec218fd1ad Mon Sep 17 00:00:00 2001 From: Vicente Mataix Ferrandiz Date: Wed, 20 Dec 2023 16:20:01 +0100 Subject: [PATCH 13/13] Suggestion by @roigcarlo --- kratos/utilities/search_utilities.h | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/kratos/utilities/search_utilities.h b/kratos/utilities/search_utilities.h index 2bf740f06d93..b974b13bd890 100644 --- a/kratos/utilities/search_utilities.h +++ b/kratos/utilities/search_utilities.h @@ -701,7 +701,6 @@ class SearchUtilities // Auxiliary variables std::size_t counter = 0; - array_1d coordinates; unsigned int i_coord; // If distributed @@ -717,11 +716,11 @@ class SearchUtilities continue; // Skip if not local } } - noalias(coordinates) = it_point->Coordinates(); + 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] = coordinates[i_coord]; + send_points_coordinates[3 * counter + i_coord] = r_coordinates[i_coord]; } ++counter; } @@ -756,14 +755,14 @@ class SearchUtilities } else { // Serial // Assign values for (auto it_point = itPointBegin ; it_point != itPointEnd ; ++it_point) { - noalias(coordinates) = it_point->Coordinates(); + 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] = coordinates[i_coord]; + rAllPointsCoordinates[3 * counter + i_coord] = r_coordinates[i_coord]; } ++counter; } @@ -813,7 +812,6 @@ class SearchUtilities // Some definitions IndexType i_coord = 0; - array_1d coordinates; // If distributed if (rDataCommunicator.IsDistributed()) { @@ -865,9 +863,9 @@ class SearchUtilities 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)) { - noalias(coordinates) = it_point->Coordinates(); + const auto& r_coordinates = it_point->Coordinates(); for (i_coord = 0; i_coord < 3; ++i_coord) { - rSearchInfo.PointCoordinates.push_back(coordinates[i_coord]); + rSearchInfo.PointCoordinates.push_back(r_coordinates[i_coord]); } rSearchInfo.Indexes.push_back(all_points_ids[i_point]); rSearchInfo.Ranks.push_back({0});