Skip to content

Commit 23e3069

Browse files
aprokopaprokop
committed
Add barycentric coordinates algorithm (triangle and tetrahedron only)
1 parent c59461f commit 23e3069

4 files changed

Lines changed: 237 additions & 0 deletions

File tree

‎src/geometry/algorithms/ArborX_BarycentricCoordinates.hpp‎

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
/****************************************************************************
2+
* Copyright (c) 2025, ArborX authors *
3+
* All rights reserved. *
4+
* *
5+
* This file is part of the ArborX library. ArborX is *
6+
* distributed under a BSD 3-clause license. For the licensing terms see *
7+
* the LICENSE file in the top-level directory. *
8+
* *
9+
* SPDX-License-Identifier: BSD-3-Clause *
10+
****************************************************************************/
11+
#ifndef ARBORX_GEOMETRY_BARYCENTRIC_COORDINATES_HPP
12+
#define ARBORX_GEOMETRY_BARYCENTRIC_COORDINATES_HPP
13+
14+
#include <ArborX_GeometryTraits.hpp>
15+
#include <misc/ArborX_Vector.hpp>
16+
17+
#include <Kokkos_Array.hpp>
18+
#include <Kokkos_Macros.hpp>
19+
20+
namespace ArborX::Details
21+
{
22+
namespace Dispatch
23+
{
24+
template <typename Tag1, typename Tag2, typename Geometry1, typename Geometry2>
25+
struct barycentric;
26+
}
27+
28+
template <typename Geometry1, typename Geometry2>
29+
KOKKOS_INLINE_FUNCTION auto barycentricCoordinates(Geometry1 const &geometry1,
30+
Geometry2 const &geometry2)
31+
{
32+
static_assert(GeometryTraits::dimension_v<Geometry1> ==
33+
GeometryTraits::dimension_v<Geometry2>);
34+
return Dispatch::barycentric<GeometryTraits::tag_t<Geometry1>,
35+
GeometryTraits::tag_t<Geometry2>, Geometry1,
36+
Geometry2>::apply(geometry1, geometry2);
37+
}
38+
39+
namespace Dispatch
40+
{
41+
42+
using namespace GeometryTraits;
43+
44+
template <typename Triangle, typename Point>
45+
struct barycentric<TriangleTag, PointTag, Triangle, Point>
46+
{
47+
KOKKOS_FUNCTION static constexpr auto apply(Triangle const &triangle,
48+
Point const &point)
49+
{
50+
using Coordinate = GeometryTraits::coordinate_type_t<Point>;
51+
52+
// ArborX should probably provide the interpolation function
53+
auto const a = triangle.a;
54+
auto const b = triangle.b;
55+
auto const c = triangle.c;
56+
57+
// Find coefficients alpha and beta such that
58+
// x = a + alpha * (b - a) + beta * (c - a)
59+
// = (1 - alpha - beta) * a + alpha * b + beta * c
60+
// recognizing the linear system
61+
// ((b - a) (c - a)) (alpha beta)^T = (x - a)
62+
Coordinate u[2] = {b[0] - a[0], b[1] - a[1]};
63+
Coordinate v[2] = {c[0] - a[0], c[1] - a[1]};
64+
Coordinate const det = v[1] * u[0] - v[0] * u[1];
65+
if (det == 0)
66+
Kokkos::abort("Degenerate triangles are not supported!");
67+
auto const inv_det = 1 / det;
68+
69+
Coordinate alpha[2] = {v[1] * inv_det, -v[0] * inv_det};
70+
Coordinate beta[2] = {-u[1] * inv_det, u[0] * inv_det};
71+
72+
auto alpha_coeff =
73+
alpha[0] * (point[0] - a[0]) + alpha[1] * (point[1] - a[1]);
74+
auto beta_coeff = beta[0] * (point[0] - a[0]) + beta[1] * (point[1] - a[1]);
75+
76+
return Kokkos::Array<Coordinate, 3>{1 - alpha_coeff - beta_coeff,
77+
alpha_coeff, beta_coeff};
78+
}
79+
};
80+
81+
template <typename Tetrahedron, typename Point>
82+
struct barycentric<TetrahedronTag, PointTag, Tetrahedron, Point>
83+
{
84+
template <typename Coordinate>
85+
using Vector = ::ArborX::Details::Vector<3, Coordinate>;
86+
87+
template <typename Coordinate>
88+
KOKKOS_FUNCTION static constexpr auto
89+
triple_product(Vector<Coordinate> const &u, Vector<Coordinate> const &v,
90+
Vector<Coordinate> const &w)
91+
{
92+
return u.dot(v.cross(w));
93+
}
94+
95+
KOKKOS_FUNCTION static constexpr auto apply(Tetrahedron const &tet,
96+
Point const &point)
97+
{
98+
using Coordinate = GeometryTraits::coordinate_type_t<Point>;
99+
100+
auto ap = point - tet.a;
101+
auto bp = point - tet.b;
102+
103+
auto ab = tet.b - tet.a;
104+
auto ac = tet.c - tet.a;
105+
auto ad = tet.d - tet.a;
106+
107+
auto bc = tet.c - tet.b;
108+
auto bd = tet.d - tet.b;
109+
110+
auto denom = 1 / triple_product(ab, ac, ad);
111+
return Kokkos::Array<Coordinate, 4>{
112+
triple_product(bp, bd, bc) * denom, triple_product(ap, ac, ad) * denom,
113+
triple_product(ap, ad, ab) * denom, triple_product(ap, ab, ac) * denom};
114+
}
115+
};
116+
117+
} // namespace Dispatch
118+
119+
} // namespace ArborX::Details
120+
121+
#endif

‎test/ArborX_EnableArrayComparison.hpp‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
/****************************************************************************
2+
* Copyright (c) 2025, ArborX authors *
3+
* All rights reserved. *
4+
* *
5+
* This file is part of the ArborX library. ArborX is *
6+
* distributed under a BSD 3-clause license. For the licensing terms see *
7+
* the LICENSE file in the top-level directory. *
8+
* *
9+
* SPDX-License-Identifier: BSD-3-Clause *
10+
****************************************************************************/
11+
12+
#ifndef ARBORX_ENABLE_ARRAY_COMPARISON_HPP
13+
#define ARBORX_ENABLE_ARRAY_COMPARISON_HPP
14+
15+
#include <Kokkos_Array.hpp>
16+
17+
#include "BoostTest_CUDA_clang_workarounds.hpp"
18+
#include <boost/test/unit_test.hpp>
19+
#include <boost/test/utils/is_forward_iterable.hpp>
20+
21+
namespace boost::unit_test
22+
{
23+
24+
template <typename T, size_t N>
25+
struct is_forward_iterable<Kokkos::Array<T, N>> : public boost::mpl::true_
26+
{};
27+
28+
template <typename T, size_t N>
29+
struct bt_iterator_traits<Kokkos::Array<T, N>, true>
30+
{
31+
using array_type = Kokkos::Array<T, N>;
32+
using value_type = typename array_type::value_type;
33+
using const_iterator = typename array_type::const_pointer;
34+
static const_iterator begin(array_type const &v) { return v.data(); }
35+
static const_iterator end(array_type const &v) { return v.data() + v.size(); }
36+
static std::size_t size(array_type const &v) { return v.size(); }
37+
};
38+
39+
} // namespace boost::unit_test
40+
41+
#endif

‎test/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,7 @@ target_include_directories(ArborX_Test_KokkosExt.exe PRIVATE ${CMAKE_CURRENT_BIN
7171
add_test(NAME ArborX_Test_KokkosExt COMMAND ArborX_Test_KokkosExt.exe)
7272

7373
add_executable(ArborX_Test_Geometry.exe
74+
tstGeometryBarycentricCoordinates.cpp
7475
tstGeometryCentroid.cpp
7576
tstGeometryDistance.cpp
7677
tstGeometryExpand.cpp

‎test/tstGeometryBarycentricCoordinates.cpp‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
/****************************************************************************
2+
* Copyright (c) 2025, ArborX authors *
3+
* All rights reserved. *
4+
* *
5+
* This file is part of the ArborX library. ArborX is *
6+
* distributed under a BSD 3-clause license. For the licensing terms see *
7+
* the LICENSE file in the top-level directory. *
8+
* *
9+
* SPDX-License-Identifier: BSD-3-Clause *
10+
****************************************************************************/
11+
12+
#include "ArborX_EnableArrayComparison.hpp"
13+
#include <ArborX_Box.hpp>
14+
#include <ArborX_GeometryTraits.hpp>
15+
#include <ArborX_Point.hpp>
16+
#include <ArborX_Tetrahedron.hpp>
17+
#include <ArborX_Triangle.hpp>
18+
#include <algorithms/ArborX_BarycentricCoordinates.hpp>
19+
20+
#include "BoostTest_CUDA_clang_workarounds.hpp"
21+
#include <boost/test/unit_test.hpp>
22+
23+
namespace tt = boost::test_tools;
24+
25+
template <int N>
26+
using Array = Kokkos::Array<float, N>;
27+
28+
BOOST_AUTO_TEST_CASE(barycentric_triangle)
29+
{
30+
using ArborX::Point;
31+
using ArborX::Triangle;
32+
using ArborX::Details::barycentricCoordinates;
33+
34+
Triangle<2> tri{{-1, -1}, {1, -1}, {-1, 1}};
35+
// clang-format off
36+
// vertices
37+
BOOST_TEST(barycentricCoordinates(tri, Point{-1.f, -1.f}) == (Array<3>{1, 0, 0}), tt::per_element());
38+
BOOST_TEST(barycentricCoordinates(tri, Point{1.f, -1.f}) == (Array<3>{0, 1, 0}), tt::per_element());
39+
BOOST_TEST(barycentricCoordinates(tri, Point{-1.f, 1.f}) == (Array<3>{0, 0, 1}), tt::per_element());
40+
// mid edges
41+
BOOST_TEST(barycentricCoordinates(tri, Point{0.f, -1.f}) == (Array<3>{0.5f, 0.5f, 0}), tt::per_element());
42+
BOOST_TEST(barycentricCoordinates(tri, Point{-1.f, 0.f}) == (Array<3>{0.5f, 0, 0.5f}), tt::per_element());
43+
BOOST_TEST(barycentricCoordinates(tri, Point{0.f, 0.f}) == (Array<3>{0, 0.5f, 0.5f}), tt::per_element());
44+
// center
45+
BOOST_TEST(barycentricCoordinates(tri, Point{-1.f/3, -1.f/3}) == (Array<3>{1.f/3, 1.f/3, 1.f/3}), tt::tolerance(1e-7f) << tt::per_element());
46+
// off-center
47+
BOOST_TEST(barycentricCoordinates(tri, Point{-0.4f, 0.2f}) == (Array<3>{0.1f, 0.3f, 0.6f}), tt::tolerance(1e-6f) << tt::per_element());
48+
// clang-format on
49+
}
50+
51+
BOOST_AUTO_TEST_CASE(barycentric_tetrahedron)
52+
{
53+
using ArborX::Point;
54+
using ArborX::Details::barycentricCoordinates;
55+
using ArborX::ExperimentalHyperGeometry::Tetrahedron;
56+
57+
Tetrahedron<> tet{{0, 0, 0}, {0, 2, 0}, {2, 0, 0}, {0, 0, 2}};
58+
// clang-format off
59+
// vertices
60+
BOOST_TEST(barycentricCoordinates(tet, Point{0.f, 0.f, 0.f}) == (Array<4>{1, 0, 0, 0}), tt::per_element());
61+
BOOST_TEST(barycentricCoordinates(tet, Point{0.f, 2.f, 0.f}) == (Array<4>{0, 1, 0, 0}), tt::per_element());
62+
BOOST_TEST(barycentricCoordinates(tet, Point{2.f, 0.f, 0.f}) == (Array<4>{0, 0, 1, 0}), tt::per_element());
63+
BOOST_TEST(barycentricCoordinates(tet, Point{0.f, 0.f, 2.f}) == (Array<4>{0, 0, 0, 1}), tt::per_element());
64+
// (some) mid edges
65+
BOOST_TEST(barycentricCoordinates(tet, Point{0.f, 1.f, 0.f}) == (Array<4>{0.5f, 0.5f, 0, 0}), tt::per_element());
66+
BOOST_TEST(barycentricCoordinates(tet, Point{0.f, 1.f, 1.f}) == (Array<4>{0, 0.5f, 0, 0.5f}), tt::per_element());
67+
BOOST_TEST(barycentricCoordinates(tet, Point{1.f, 0.f, 1.f}) == (Array<4>{0, 0, 0.5f, 0.5f}), tt::per_element());
68+
// (some) mid faces
69+
BOOST_TEST(barycentricCoordinates(tet, Point{2.f/3, 2.f/3, 0.f}) == (Array<4>{1.f/3, 1.f/3, 1.f/3, 0}), tt::tolerance(1e-6f) << tt::per_element());
70+
BOOST_TEST(barycentricCoordinates(tet, Point{2.f/3, 2.f/3, 2.f/3}) == (Array<4>{0, 1.f/3, 1.f/3, 1.f/3}), tt::tolerance(1e-7f) << tt::per_element());
71+
// center
72+
BOOST_TEST(barycentricCoordinates(tet, Point{0.5f, 0.5f, 0.5f}) == (Array<4>{0.25f, 0.25f, 0.25f, 0.25f}), tt::tolerance(1e-7f) << tt::per_element());
73+
// clang-format on
74+
}

0 commit comments

Comments
 (0)