Add Tetrahedron geometry (arborx#1079)

Add Tetrahedron geometry

src/geometry/ArborX_GeometryTraits.hpp

@@ -61,14 +61,16 @@ DEFINE_GEOMETRY(box, BoxTag);
 DEFINE_GEOMETRY(sphere, SphereTag);
 DEFINE_GEOMETRY(triangle, TriangleTag);
 DEFINE_GEOMETRY(kdop, KDOPTag);
+DEFINE_GEOMETRY(tetrahedron, TetrahedronTag);
 DEFINE_GEOMETRY(ray, RayTag);
 #undef DEFINE_GEOMETRY
 
 template <typename Geometry>
 inline constexpr bool
     is_valid_geometry = (is_point_v<Geometry> || is_box_v<Geometry> ||
                          is_sphere_v<Geometry> || is_kdop_v<Geometry> ||
-                         is_triangle_v<Geometry> || is_ray_v<Geometry>);
+                         is_triangle_v<Geometry> ||
+                         is_tetrahedron_v<Geometry> || is_ray_v<Geometry>);
 
 template <typename Geometry>
 using DimensionNotSpecializedArchetypeAlias =

src/geometry/ArborX_Tetrahedron.hpp

@@ -0,0 +1,109 @@
+/****************************************************************************
+ * Copyright (c) 2017-2022 by the ArborX authors                            *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of the ArborX library. ArborX is                       *
+ * distributed under a BSD 3-clause license. For the licensing terms see    *
+ * the LICENSE file in the top-level directory.                             *
+ *                                                                          *
+ * SPDX-License-Identifier: BSD-3-Clause                                    *
+ ****************************************************************************/
+
+#ifndef ARBORX_TETRAHEDRON_HPP
+#define ARBORX_TETRAHEDRON_HPP
+
+#include <ArborX_DetailsAlgorithms.hpp>
+#include <ArborX_GeometryTraits.hpp>
+#include <ArborX_HyperPoint.hpp>
+
+#include <Kokkos_Array.hpp>
+
+namespace ArborX::ExperimentalHyperGeometry
+{
+// Need to add a protection that the points are not in the same plane
+template <class Coordinate = float>
+struct Tetrahedron
+{
+  ExperimentalHyperGeometry::Point<3, Coordinate> a;
+  ExperimentalHyperGeometry::Point<3, Coordinate> b;
+  ExperimentalHyperGeometry::Point<3, Coordinate> c;
+  ExperimentalHyperGeometry::Point<3, Coordinate> d;
+};
+
+template <class Coordinate>
+Tetrahedron(ExperimentalHyperGeometry::Point<3, Coordinate>,
+            ExperimentalHyperGeometry::Point<3, Coordinate>,
+            ExperimentalHyperGeometry::Point<3, Coordinate>,
+            ExperimentalHyperGeometry::Point<3, Coordinate>)
+    -> Tetrahedron<Coordinate>;
+
+} // namespace ArborX::ExperimentalHyperGeometry
+
+template <class Coordinate>
+struct ArborX::GeometryTraits::dimension<
+    ArborX::ExperimentalHyperGeometry::Tetrahedron<Coordinate>>
+{
+  static constexpr int value = 3;
+};
+template <class Coordinate>
+struct ArborX::GeometryTraits::tag<
+    ArborX::ExperimentalHyperGeometry::Tetrahedron<Coordinate>>
+{
+  using type = TetrahedronTag;
+};
+template <class Coordinate>
+struct ArborX::GeometryTraits::coordinate_type<
+    ArborX::ExperimentalHyperGeometry::Tetrahedron<Coordinate>>
+{
+  using type = Coordinate;
+};
+
+namespace ArborX::Details::Dispatch
+{
+using GeometryTraits::BoxTag;
+using GeometryTraits::PointTag;
+using GeometryTraits::TetrahedronTag;
+
+// expand a box to include a tetrahedron
+template <typename Box, typename Tetrahedron>
+struct expand<BoxTag, TetrahedronTag, Box, Tetrahedron>
+{
+  KOKKOS_FUNCTION static void apply(Box &box, Tetrahedron const &tet)
+  {
+    Details::expand(box, tet.a);
+    Details::expand(box, tet.b);
+    Details::expand(box, tet.c);
+    Details::expand(box, tet.d);
+  }
+};
+
+template <typename Point, typename Tetrahedron>
+struct intersects<PointTag, TetrahedronTag, Point, Tetrahedron>
+{
+  KOKKOS_FUNCTION static constexpr bool 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. This is a simple but possibly not very
+    // efficient algorithm.
+    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 false;
+    }
+    return true;
+  }
+};
+
+} // namespace ArborX::Details::Dispatch
+
+#endif

test/tstDetailsAlgorithms.cpp

@@ -15,6 +15,7 @@
 #include <ArborX_HyperPoint.hpp>
 #include <ArborX_HyperSphere.hpp>
 #include <ArborX_HyperTriangle.hpp>
+#include <ArborX_Tetrahedron.hpp>
 
 #include <boost/mpl/list.hpp>
 
@@ -25,6 +26,7 @@ using Point = ArborX::ExperimentalHyperGeometry::Point<3>;
 using Box = ArborX::ExperimentalHyperGeometry::Box<3>;
 using Sphere = ArborX::ExperimentalHyperGeometry::Sphere<3>;
 using Triangle = ArborX::ExperimentalHyperGeometry::Triangle<3>;
+using Tetrahedron = ArborX::ExperimentalHyperGeometry::Tetrahedron<>;
 
 BOOST_AUTO_TEST_CASE(distance_point_point)
 {
@@ -276,6 +278,19 @@ BOOST_AUTO_TEST_CASE(intersects)
       ArborX::ExperimentalHyperGeometry::Triangle<3>{
           {{.1, .1, .1}}, {{.1, .9, .1}}, {{.9, .1, .1}}},
       unit_box));
+
+  // tetrahedron
+  constexpr Tetrahedron tet{{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+  BOOST_TEST(intersects(Point{0, 0, 0}, tet));
+  BOOST_TEST(intersects(Point{1, 0, 0}, tet));
+  BOOST_TEST(intersects(Point{0, 1, 0}, tet));
+  BOOST_TEST(intersects(Point{0, 0, 1}, tet));
+  BOOST_TEST(intersects(Point{0.2, 0.2, 0.1}, tet));
+  BOOST_TEST(!intersects(Point{-0.1, 0, 0}, tet));
+  BOOST_TEST(!intersects(Point{0, -0.1, 0}, tet));
+  BOOST_TEST(!intersects(Point{0, 0, 1.1}, tet));
+  BOOST_TEST(!intersects(Point{0.5, 0.5, 0.5}, tet));
+  BOOST_TEST(!intersects(Point{-0.5, 0.5, 0.5}, tet));
 }
 
 BOOST_AUTO_TEST_CASE(equals)
@@ -330,6 +345,13 @@ BOOST_AUTO_TEST_CASE(expand)
   BOOST_TEST(equals(box, Box{{{-24., -24., -24.}}, {{24., 24., 24.}}}));
   expand(box, Triangle{{{0, 0, 0}}, {{48, 0, 0}}, {{0, 48, 0}}});
   BOOST_TEST(equals(box, Box{{{-24., -24., -24.}}, {{48., 48., 24.}}}));
+
+  // expand box with tetrahedrons
+  box = Box{};
+  expand(box, Tetrahedron{{-1, -2, 3}, {1, 3, 2}, {0, 3, 7}, {-5, 4, 7}});
+  BOOST_TEST(equals(box, Box{{-5, -2, 2}, {1, 4, 7}}));
+  expand(box, Tetrahedron{{-3, -5, 2}, {2, 6, -1}, {3, 2, 3}, {5, 8, -3}});
+  BOOST_TEST(equals(box, Box{{-5, -5, -3}, {5, 8, 7}}));
 }
 
 BOOST_AUTO_TEST_CASE(centroid)