Merge pull request #1256 from aprokop/distance_point_tet

Author: aprokop

src/geometry/algorithms/ArborX_Distance.hpp

@@ -13,8 +13,10 @@
 
 #include "ArborX_Equals.hpp"
 #include <ArborX_GeometryTraits.hpp>
+#include <ArborX_Triangle.hpp>
 #include <misc/ArborX_Vector.hpp>
 
+#include <Kokkos_Array.hpp>
 #include <Kokkos_Clamp.hpp>
 #include <Kokkos_MathematicalFunctions.hpp>
 #include <Kokkos_MinMax.hpp>
@@ -198,6 +200,50 @@ struct distance<PointTag, TriangleTag, Point, Triangle>
   }
 };
 
+// distance point-tetrahedron
+template <typename Point, typename Tetrahedron>
+struct distance<PointTag, TetrahedronTag, Point, Tetrahedron>
+{
+  static constexpr int DIM = dimension_v<Point>;
+  using Coordinate = coordinate_type_t<Tetrahedron>;
+
+  KOKKOS_FUNCTION static auto apply(Point const &point, Tetrahedron const &tet)
+  {
+    static_assert(GeometryTraits::dimension_v<Point> == 3);
+
+    constexpr int N = 4;
+    Kokkos::Array<decltype(tet.a), N> v = {tet.a, tet.b, tet.c, tet.d};
+
+    // For every plane check that the vertex lies within the same halfspace as
+    // the other tetrahedron vertex (similar to the current intersects
+    // algorithm). If not, compute the distance to the corresponding triangle.
+    for (int j = 0; j < N; ++j)
+    {
+      auto normal = (v[(j + 1) % N] - v[j]).cross(v[(j + 2) % N] - v[j]);
+
+      bool same_half_space =
+          (normal.dot(v[(j + 3) % N] - v[j]) * normal.dot(point - v[j]) >= 0);
+      if (!same_half_space)
+        return Details::distance(
+            point, Triangle{v[j], v[(j + 1) % N], v[(j + 2) % N]});
+    }
+    return static_cast<Coordinate>(0);
+  }
+};
+
+// distance tetrahedron-point
+template <typename Tetrahedron, typename Point>
+struct distance<TetrahedronTag, PointTag, Tetrahedron, Point>
+{
+  static constexpr int DIM = dimension_v<Point>;
+  using Coordinate = coordinate_type_t<Tetrahedron>;
+
+  KOKKOS_FUNCTION static auto apply(Tetrahedron const &tet, Point const &p)
+  {
+    return Details::distance(p, tet);
+  }
+};
+
 // distance box-box
 template <typename Box1, typename Box2>
 struct distance<BoxTag, BoxTag, Box1, Box2>

test/tstDetailsAlgorithms.cpp

@@ -161,6 +161,39 @@ BOOST_AUTO_TEST_CASE(distance_point_triangle)
   BOOST_TEST(distance(Point3{1, -1, -1}, triangle3_2) == std::sqrt(3.f));
 }
 
+BOOST_AUTO_TEST_CASE(distance_point_tetrahedron)
+{
+  using ArborX::Details::distance;
+
+  using Point = ArborX::Point<3>;
+
+  constexpr ArborX::ExperimentalHyperGeometry::Tetrahedron tet{
+      Point{0, 0, 0}, Point{1, 0, 0}, Point{0, 1, 0}, Point{0, 0, 1}};
+
+  // vertices
+  BOOST_TEST(distance(Point{0, 0, 0}, tet) == 0);
+  BOOST_TEST(distance(Point{1, 0, 0}, tet) == 0);
+  BOOST_TEST(distance(Point{0, 1, 0}, tet) == 0);
+  BOOST_TEST(distance(Point{0, 0, 1}, tet) == 0);
+
+  // inside
+  BOOST_TEST(distance(Point{0.4f, 0.4f, 0.1f}, tet) == 0);
+  BOOST_TEST(distance(Point{0.8f, 0.05f, 0.05f}, tet) == 0);
+
+  // same plane as some side
+  BOOST_TEST(distance(Point{2, 0, 0}, tet) == 1);
+  BOOST_TEST(distance(Point{0.5f, -0.5f, 0}, tet) == 0.5f);
+  BOOST_TEST(distance(Point{-0.5f, 0.5f, 0}, tet) == 0.5f);
+  BOOST_TEST(distance(Point{0, 0, 2}, tet) == 1);
+  BOOST_TEST(distance(Point{0, -0.5f, 0}, tet) == 0.5f);
+
+  // outside
+  BOOST_TEST(distance(Point{-1, -1, -1}, tet) == std::sqrt(3.f));
+  BOOST_TEST(distance(Point{-1, -1, 2}, tet) == std::sqrt(3.f));
+  BOOST_TEST(distance(Point{1.5f, 1.5f, -1}, tet) == std::sqrt(3.f));
+  BOOST_TEST(distance(Point{1.5f, 1.5f, 0.5f}, tet) == 1.5f);
+}
+
 BOOST_AUTO_TEST_CASE(distance_box_box)
 {
   using ArborX::Details::distance;