[Core] Adding SearchInBoundingBox into GeometricalObjectsBins

kratos/python/add_search_strategies_to_python.cpp

@@ -575,6 +575,34 @@ void AddSearchStrategiesToPython(pybind11::module& m)
         }
         return list_results;
     })
+    .def("SearchInBoundingBox", [&](GeometricalObjectsBins& self, const Point& rPoint, const array_1d<double, 3>& rMinPoint, const array_1d<double, 3>& rMaxPoint) {
+        // Perform the search
+        std::vector<ResultTypeGeometricalObject> results;
+        self.SearchInBoundingBox(rPoint, rMinPoint, rMaxPoint, results);
+
+        // Copy the results to the python list
+        py::list list_results;
+        for (auto& r_result : results) {
+            list_results.append(r_result);
+        }
+        return list_results;
+    })
+    .def("SearchInBoundingBox", [&](GeometricalObjectsBins& self, const NodesContainerType& rNodes, const array_1d<double, 3>& rMinPoint, const array_1d<double, 3>& rMaxPoint) {
+        // Perform the search
+        std::vector<std::vector<ResultTypeGeometricalObject>> results;
+        self.SearchInBoundingBox(rNodes.begin(), rNodes.end(), rMinPoint, rMaxPoint, results);
+
+        // Copy the results to the python list
+        py::list list_results;
+        for (auto& r_result : results) {
+            py::list sub_list_results;
+            for (auto& r_sub_result : r_result) {
+                sub_list_results.append(r_sub_result);
+            }
+            list_results.append(sub_list_results);
+        }
+        return list_results;
+    })
     .def("SearchNearest", [&](GeometricalObjectsBins& self, const Point& rPoint) {
         // Perform the search
         return self.SearchNearest(rPoint);

kratos/spatial_containers/geometrical_objects_bins.cpp

@@ -158,6 +158,91 @@ GeometricalObjectsBins::ResultType GeometricalObjectsBins::SearchNearestInRadius
 /***********************************************************************************/
 /***********************************************************************************/
 
+void GeometricalObjectsBins::SearchInBoundingBox(
+    const PointType& rPoint,
+    const array_1d<double, 3>& rMinPoint,
+    const array_1d<double, 3>& rMaxPoint,
+    std::vector<ResultType>& rResults
+    )
+{
+    // Clear the results
+    rResults.clear();
+
+    // Define bounding box
+    BoundingBox<PointType> bounding_box;
+    auto& r_min_point = bounding_box.GetMinPoint();
+    auto& r_max_point = bounding_box.GetMaxPoint();
+    noalias(r_min_point.Coordinates()) = rMinPoint;
+    noalias(r_max_point.Coordinates()) = rMaxPoint;
+
+    // Initialize the processed objects and collected results
+    std::unordered_set<GeometricalObject*> processed_objects;
+    std::unordered_set<GeometricalObject*> candidates;
+    std::vector<GeometricalObject*> collected_results;
+
+    // Initialize the position bounds
+    array_1d<std::size_t, Dimension> min_position;
+    array_1d<std::size_t, Dimension> max_position;
+
+    // Calculate the position bounds
+    for(unsigned int i = 0; i < Dimension; i++ ) {
+        min_position[i] = CalculatePosition(rMinPoint[i], i);
+        max_position[i] = CalculatePosition(rMaxPoint[i], i) + 1;
+    }
+
+    // Loop over the cells and gather candidates
+    for(std::size_t k = min_position[2]; k < max_position[2]; k++) {
+        for(std::size_t j = min_position[1]; j < max_position[1]; j++) {
+            for(std::size_t i = min_position[0]; i < max_position[0]; i++) {
+                auto& r_cell = GetCell(i, j, k);
+
+                // Check if the cell's bounding box intersects with the search bounding box
+                if (IsCellInsideBoundingBox(i, j, k, rPoint, bounding_box)) {
+                    for(auto p_geometrical_object : r_cell) {
+                        // Insert the object into the processed set to avoid duplicates
+                        if (processed_objects.insert(p_geometrical_object).second) {
+                            // Collect the object
+                            collected_results.push_back(p_geometrical_object);
+                        }
+                    }
+                } else {
+                    for(auto p_geometrical_object : r_cell) {
+                        candidates.insert(p_geometrical_object);
+                    }
+                }
+            }
+        }
+    }
+
+    // Loop over all candidates and check for intersection
+    for(auto& p_geometrical_object : candidates) {
+        // Skip already processed objects
+        if (processed_objects.find(p_geometrical_object) != processed_objects.end()) {
+            continue;
+        }
+
+        const auto& r_geometry = p_geometrical_object->GetGeometry();
+        if(r_geometry.HasIntersection(r_min_point, r_max_point)) {
+            // Collect the object
+            collected_results.push_back(p_geometrical_object);
+        }
+    }
+
+    // Reserve space in rResults to prevent reallocations
+    rResults.reserve(collected_results.size());
+
+    // Insert all collected results into rResults
+    for (auto& p_geometrical_object : collected_results) {
+        rResults.push_back(ResultType(p_geometrical_object));
+    }
+
+    // Note: We avoid shrink_to_fit for potential performance issues
+    // rResults.shrink_to_fit();
+}
+
+/***********************************************************************************/
+/***********************************************************************************/
+
 GeometricalObjectsBins::ResultType GeometricalObjectsBins::SearchNearest(const PointType& rPoint)
 {
     ResultType current_result;
@@ -190,6 +275,41 @@ GeometricalObjectsBins::ResultType GeometricalObjectsBins::SearchIsInside(const
 /***********************************************************************************/
 /***********************************************************************************/
 
+bool GeometricalObjectsBins::IsCellInsideBoundingBox(
+    const std::size_t I,
+    const std::size_t J,
+    const std::size_t K,
+    const Point& rPoint,
+    const BoundingBox<PointType>& rBoundingBox
+    )
+{
+    // Check the indices
+    KRATOS_DEBUG_ERROR_IF(I > mNumberOfCells[0]) << "Index " << I << " is larger than number of cells in x direction : " << mNumberOfCells[0] << std::endl;
+    KRATOS_DEBUG_ERROR_IF(J > mNumberOfCells[1]) << "Index " << J << " is larger than number of cells in y direction : " << mNumberOfCells[1] << std::endl;
+    KRATOS_DEBUG_ERROR_IF(K > mNumberOfCells[2]) << "Index " << K << " is larger than number of cells in z direction : " << mNumberOfCells[2] << std::endl;
+
+    // Get the bounding box points min point
+    const array_1d<double, 3>& r_min_point = mBoundingBox.GetMinPoint();
+
+    // Calculate the minimum point in the bounding box
+    array_1d<double, 3> min_point;
+    min_point[0] = r_min_point[0] + I * mCellSizes[0];
+    min_point[1] = r_min_point[1] + J * mCellSizes[1];
+    min_point[2] = r_min_point[2] + K * mCellSizes[2];
+
+    // Calculate the maximum point in the bounding box
+    array_1d<double, 3> max_point;
+    max_point[0] = r_min_point[0] + (I + 1) * mCellSizes[0];
+    max_point[1] = r_min_point[1] + (J + 1) * mCellSizes[1];
+    max_point[2] = r_min_point[2] + (K + 1) * mCellSizes[2];
+
+    // We do the checks of both points
+    return rBoundingBox.IsInside(min_point) && rBoundingBox.IsInside(max_point);
+}
+
+/***********************************************************************************/
+/***********************************************************************************/
+
 bool GeometricalObjectsBins::PointIsInsideBoundingBox(const array_1d<double, 3>& rCoords)
 {
     // Get the bounding box points

kratos/spatial_containers/geometrical_objects_bins.h

@@ -222,6 +222,47 @@ class KRATOS_API(KRATOS_CORE) GeometricalObjectsBins
         return results;
     }
 
+    /**
+     * @brief This method takes a point and finds all of the objects in the given bounding box to it.
+     * @details The result contains the object and also its distance to the point.
+     * @param rPoint The point to be checked
+     * @param rMinPoint The minimum point of the bounding box
+     * @param rMaxPoint The maximum point of the bounding box
+     * @param rResults The results of the search
+     */
+    void SearchInBoundingBox(
+        const PointType& rPoint,
+        const array_1d<double, 3>& rMinPoint,
+        const array_1d<double, 3>& rMaxPoint,
+        std::vector<ResultType>& rResults
+        );
+
+    /**
+     * @brief This method takes a point and finds all of the objects in the given bounding box to it (iterative version).
+     * @details The result contains the object and also its distance to the point.
+     * @param itPointBegin The first point iterator
+     * @param itPointEnd The last point iterator
+     * @param rMinPoint The minimum point of the bounding box
+     * @param rMaxPoint The maximum point of the bounding box
+     * @param rResults The results of the search
+     * @tparam TPointIteratorType The type of the point iterator
+     */
+    template<typename TPointIteratorType>
+    void SearchInBoundingBox(
+        TPointIteratorType itPointBegin,
+        TPointIteratorType itPointEnd,
+        const array_1d<double, 3>& rMinPoint,
+        const array_1d<double, 3>& rMaxPoint,
+        std::vector<std::vector<ResultType>>& rResults
+        )
+    {
+        const std::size_t number_of_points = std::distance(itPointBegin, itPointEnd);
+        rResults.resize(number_of_points);
+        for (auto it_point = itPointBegin ; it_point != itPointEnd ; it_point++){
+            SearchInBoundingBox(*it_point, rMinPoint, rMaxPoint, rResults[it_point - itPointBegin]);
+        }
+    }
+
     /**
      * @brief This method takes a point and finds the nearest object to it.
      * @details If there are more than one object in the same minimum distance only one is returned
@@ -386,6 +427,24 @@ class KRATOS_API(KRATOS_CORE) GeometricalObjectsBins
     ///@name Protected Operations
     ///@{
 
+    /**
+     * @brief Checks if a 3D bounding box is fully inside a given bounding box from a point.
+     * @details This function determines if a 3D bounding box, defined by the indices of the cell studied, is fully within a specified bounding box from a given point. It calculates the farthest point in the bounding box from the reference point and checks if the distance from this point to the reference point is within the specified bounding box.
+     * @param I The index in x direction
+     * @param J The index in y direction
+     * @param K The index in z direction
+     * @param rPoint The reference point from which the distance is measured.
+     * @param rBoundingBox The bounding box within which the bounding box should be.
+     * @return True if the bounding box is fully inside the bounding box, false otherwise.
+     */
+    bool IsCellInsideBoundingBox(
+        const std::size_t I,
+        const std::size_t J,
+        const std::size_t K,
+        const Point& rPoint,
+        const BoundingBox<PointType>& rBoundingBox
+        );
+
     /**
      * @brief This method checks if a point is inside any bounding box of the global bounding boxes
      * @param rCoords The coordinates of the point

kratos/tests/cpp_tests/spatial_containers/test_geometrical_objects_bins.cpp

@@ -170,6 +170,77 @@ KRATOS_TEST_CASE_IN_SUITE(GeometricalObjectsBinsSearchNearestInRadius, KratosCor
     KRATOS_EXPECT_TRUE(id == 3);
 }
 
+/** Checks bins search in bounding box
+*/
+KRATOS_TEST_CASE_IN_SUITE(GeometricalObjectsBinsSearchInBoundingBox, KratosCoreFastSuite)
+{
+    Model current_model;
+
+    // Generate the cube skin
+    ModelPart& r_skin_part = CppTestsUtilities::CreateCubeSkinModelPart(current_model);
+
+    GeometricalObjectsBins bins(r_skin_part.ElementsBegin(), r_skin_part.ElementsEnd());
+
+    std::vector<GeometricalObjectsBins::ResultType> results;
+    Point center_point{0.00,0.00,0.00};
+
+    const double limit_1 = 0.29;
+    Point point_min{-limit_1,-limit_1,-limit_1};
+    Point point_max{limit_1,limit_1,limit_1};
+    bins.SearchInBoundingBox(center_point, point_min, point_max, results);
+    KRATOS_EXPECT_EQ(results.size(), 0);
+
+    const double limit_2 = 0.3;
+    point_min.X() = -limit_2;
+    point_min.Y() = -limit_2;
+    point_min.Z() = -limit_2;
+    point_max.X() = limit_2;
+    point_max.Y() = limit_2;
+    point_max.Z() = limit_2;
+    bins.SearchInBoundingBox(center_point, point_min, point_max, results);
+    KRATOS_EXPECT_EQ(results.size(), 4);
+
+    const double limit_3 = 0.4;
+    point_min.X() = -limit_3;
+    point_min.Y() = -limit_3;
+    point_min.Z() = -limit_3;
+    point_max.X() = limit_3;
+    point_max.Y() = limit_3;
+    point_max.Z() = limit_3;
+    bins.SearchInBoundingBox(center_point, point_min, point_max, results);
+    KRATOS_EXPECT_EQ(results.size(), 4);
+
+    const double limit_4 = 0.6;
+    point_min.X() = -limit_4;
+    point_min.Y() = -limit_4;
+    point_min.Z() = -limit_4;
+    point_max.X() = limit_4;
+    point_max.Y() = limit_4;
+    point_max.Z() = limit_4;
+    bins.SearchInBoundingBox(center_point, point_min, point_max, results);
+    KRATOS_EXPECT_EQ(results.size(), 8);
+
+    const double limit_5 = 0.7;
+    point_min.X() = -limit_5;
+    point_min.Y() = -limit_5;
+    point_min.Z() = -limit_5;
+    point_max.X() = limit_5;
+    point_max.Y() = limit_5;
+    point_max.Z() = limit_5;
+    bins.SearchInBoundingBox(center_point, point_min, point_max, results);
+    KRATOS_EXPECT_EQ(results.size(), 8);
+
+    const double limit_6 = 0.9;
+    point_min.X() = -limit_6;
+    point_min.Y() = -limit_6;
+    point_min.Z() = -limit_6;
+    point_max.X() = limit_6;
+    point_max.Y() = limit_6;
+    point_max.Z() = limit_6;
+    bins.SearchInBoundingBox(center_point, point_min, point_max, results);
+    KRATOS_EXPECT_EQ(results.size(), 12);
+}
+
 /** Checks bins search nearest
 */
 KRATOS_TEST_CASE_IN_SUITE(GeometricalObjectsBinsSearchNearest, KratosCoreFastSuite) 

kratos/tests/test_geometrical_object_bins.py

@@ -106,6 +106,29 @@ def test_GeometricalObjectsBins_SearchNearestInRadius(self):
             self.assertEqual(len(results), 1)
             self.assertEqual(results[0].Get().Id, 1)
 
+    def test_GeometricalObjectsBins_SearchInBoundingBox(self):
+        # Define radius
+        point_min = KM.Point(-0.25, -0.25, -0.10)
+        point_max = KM.Point(0.25, 0.25, 0.40)
+
+        # Reference solution
+        cond_id_ref = [34, 12, 117, 89, 14, 116, 120, 103, 88, 2, 98, 95, 112, 51, 110, 77, 92]
+
+        # One node search
+        if not KM.IsDistributedRun():
+            results = self.search.SearchInBoundingBox(self.node, point_min, point_max)
+            self.assertEqual(len(results), 17)
+            for result in results:
+                self.assertTrue(result.Get().Id in cond_id_ref)
+
+        # Nodes array search
+        results = self.search.SearchInBoundingBox(self.sub_model_part.Nodes, point_min, point_max)
+        if self.is_first_partition():
+            self.assertEqual(len(results), 1)
+            self.assertEqual(len(results[0]), 17)
+            for result in results[0]:
+                self.assertTrue(result.Get().Id in cond_id_ref)
+
     def test_GeometricalObjectsBins_SearchNearest(self):
         # One node search
         if not KM.IsDistributedRun():