Add intersects ellipsoid-point

Author: aprokop

src/geometry/algorithms/ArborX_Intersects.hpp

@@ -14,6 +14,7 @@
 #include "ArborX_Distance.hpp"
 #include "ArborX_Expand.hpp"
 #include <ArborX_GeometryTraits.hpp>
+#include <misc/ArborX_Vector.hpp>
 
 #include <Kokkos_Array.hpp>
 #include <Kokkos_MathematicalFunctions.hpp>
@@ -512,6 +513,43 @@ struct intersects<BoxTag, SegmentTag, Box, Segment>
   }
 };
 
+namespace
+{
+// Computes x^t R y
+template <int DIM, typename Coordinate>
+KOKKOS_INLINE_FUNCTION auto rmt_multiply(Details::Vector<DIM, Coordinate> x,
+                                         Coordinate const rmt[DIM][DIM],
+                                         Details::Vector<DIM, Coordinate> y)
+{
+  Coordinate r = 0;
+  for (int i = 0; i < DIM; ++i)
+    for (int j = 0; j < DIM; ++j)
+      r += x[i] * rmt[i][j] * y[j];
+  return r;
+}
+} // namespace
+
+template <typename Ellipsoid, typename Point>
+struct intersects<EllipsoidTag, PointTag, Ellipsoid, Point>
+{
+  KOKKOS_FUNCTION static constexpr bool apply(Ellipsoid const &ellipsoid,
+                                              Point const &point)
+  {
+    auto d = point - ellipsoid.centroid();
+    return rmt_multiply(d, ellipsoid.rmt(), d) <= 1;
+  }
+};
+
+template <typename Point, typename Ellipsoid>
+struct intersects<PointTag, EllipsoidTag, Point, Ellipsoid>
+{
+  KOKKOS_FUNCTION static constexpr bool apply(Point const &point,
+                                              Ellipsoid const &ellipsoid)
+  {
+    return Details::intersects(ellipsoid, point);
+  }
+};
+
 } // namespace Dispatch
 
 } // namespace ArborX::Details

test/tstDetailsAlgorithms.cpp

@@ -342,6 +342,21 @@ BOOST_AUTO_TEST_CASE(intersects)
   BOOST_TEST(!intersects(Segment2{{0.5, 1.6}, {2, 0}}, box2));
   BOOST_TEST(intersects(box2, Segment2{{-1, 2}, {2, -1}}));
   BOOST_TEST(!intersects(Segment2{{0.5, 1.6}, {2, 0}}, box2));
+
+  // ellipsoid
+  using Ellipse = ArborX::Experimental::Ellipsoid<2>;
+  constexpr Ellipse ellipse{{1.f, 1.f}, {{2.f, -1.5f}, {-1.5f, 2.f}}};
+  BOOST_TEST(intersects(ellipse, Point2{1.f, 1.f}));
+  BOOST_TEST(intersects(ellipse, Point2{0.f, 0.f}));
+  BOOST_TEST(!intersects(ellipse, Point2{-0.01f, -0.01f}));
+  BOOST_TEST(intersects(ellipse, Point2{0.f, 0.5f}));
+  BOOST_TEST(!intersects(ellipse, Point2{0.f, 0.6f}));
+  BOOST_TEST(intersects(ellipse, Point2{2.f, 2.f}));
+  BOOST_TEST(!intersects(ellipse, Point2{0.5f, 1.5f}));
+  BOOST_TEST(intersects(ellipse, Point2{1.f, 0.3f}));
+  BOOST_TEST(!intersects(ellipse, Point2{1.f, 0.29f}));
+  BOOST_TEST(intersects(ellipse, Point2{1.f, 1.69f}));
+  BOOST_TEST(intersects(ellipse, Point2{1.f, 1.70f}));
 }
 
 BOOST_AUTO_TEST_CASE(equals)