From b3e38da5602f908bd7cc101e6de057fbb0f3e2ca Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Thu, 11 Dec 2025 19:08:08 -0600 Subject: [PATCH 01/14] Started Toroid header file with brief description --- src/orange/surf/Toroid.hh | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 src/orange/surf/Toroid.hh diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh new file mode 100644 index 0000000000..6a4a7b03b7 --- /dev/null +++ b/src/orange/surf/Toroid.hh @@ -0,0 +1,34 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file orange/surf/Toroid.hh +//---------------------------------------------------------------------------// +#pragma once + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Z-aligned Elliptical Toroid. + * + * An elliptical toroid is a shape created by revolving an axis-aligned ellipse + * around a central axis. This shape can be used in everything from pipe bends + * to tokamaks in fusion reactors. Possesses a major radius r, and ellipse + * radii a and b, as shown in the below diagram: + * ___ _________ ___ + * / | \ / \ + * / b \ / \ + * | | | | | + * |-a--+ | o-----r--+ | + * \ / \ / + * \ / \ / + * ⁻⁻⁻ ⁻⁻⁻⁻⁻⁻⁻⁻⁻ ⁻⁻⁻ + * + * This torus can be defined with the following quartic equation: + * \f[ + * (x^2 + y^2 + p*y^2 + B_0) - A_0 * (x^2 + y^2) = 0 + * \f] + * where \f[p = a^2/b^2, A_0 = 4*r^2, and B_0 = (r^2-a^2)\f]. + */ +} // namespace celeritas From ed902ef5bab39408c27052a19e818e650e124ad1 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Thu, 11 Dec 2025 22:31:11 -0600 Subject: [PATCH 02/14] Added declarations for toroid functions and 'tor' SurfaceType --- src/orange/OrangeTypes.hh | 1 + src/orange/surf/Toroid.hh | 77 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+) diff --git a/src/orange/OrangeTypes.hh b/src/orange/OrangeTypes.hh index e43efd390f..bdf28e6e7c 100644 --- a/src/orange/OrangeTypes.hh +++ b/src/orange/OrangeTypes.hh @@ -145,6 +145,7 @@ enum class SurfaceType : unsigned char sq, //!< Simple quadric gq, //!< General quadric inv, //!< Involute + tor, //!< Toroid size_ //!< Sentinel value for number of surface types }; diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh index 6a4a7b03b7..d892109c9a 100644 --- a/src/orange/surf/Toroid.hh +++ b/src/orange/surf/Toroid.hh @@ -6,6 +6,10 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/cont/Array.hh" +#include "corecel/cont/Span.hh" +#include "orange/OrangeTypes.hh" + namespace celeritas { //---------------------------------------------------------------------------// @@ -31,4 +35,77 @@ namespace celeritas * \f] * where \f[p = a^2/b^2, A_0 = 4*r^2, and B_0 = (r^2-a^2)\f]. */ +class Toroid +{ + public: + //@{ + //! \name Type aliases + using Intersections = Array; + using StorageSpan = Span; + using Real3 = Array; + //@} + + public: + //// CLASS ATTRIBUTES //// + + // Surface type identifier + static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type() + { + return SurfaceType::tor; + } + + //! Safety distance is calculable w/xy of normal and ellipse safety + //! distance, but this is out of scope at first and might not be trivially + //! calculable + //! https://web.archive.org/web/20170829172516/https://www.spaceroots.org/documents/distance/distance-to-ellipse.pdf + static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; } + + public: + //// CONSTRUCTORS //// + + explicit Toroid(Real3 const& origin, + real_type major_radius, + real_type ellipse_xy_radius, + real_type ellipse_z_radius); + + // Construct from raw data + template + explicit inline CELER_FUNCTION Toroid(Span); + + //// ACCESSORS //// + + //! Center of the toroid (in the donut hole) + CELER_FUNCTION Real3 const& origin() const { return origin_; } + + //! Radius from origin to center of revolved ellipse + CELER_FUNCTION real_type major_radius() const { return r_; } + + //! Radius of revolved ellipse along xy plane + CELER_FUNCTION real_type ellipse_xy_radius() const { return a_; } + + //! Radius of revolved ellipse along z axis + CELER_FUNCTION real_type ellipse_z_radius() const { return b_; } + + //// CALCULATION /// + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const; + + // Calculate all possible straight-line intersections with this surface + inline CELER_FUNCTION Intersections calc_intersections( + Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const; + + // Calculate outward normal at a position + inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const; + + private: + // Location of center of toroid + Real3 origin_; + + // Radii + real_type r_; // Radius from origin to center of revolved ellipse (along + // xy plane) + real_type a_; // Horizontal radius of revolved ellipse (along xy plane) + real_type b_; // Vertical radius of revolved ellipse (along z axis) +}; } // namespace celeritas From 200205c2ee799a0e646479075081bce1fbf35e47 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sat, 13 Dec 2025 01:07:43 -0600 Subject: [PATCH 03/14] Aded toroid cc file --- src/orange/surf/Toroid.cc | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 src/orange/surf/Toroid.cc diff --git a/src/orange/surf/Toroid.cc b/src/orange/surf/Toroid.cc new file mode 100644 index 0000000000..f678c7c9a9 --- /dev/null +++ b/src/orange/surf/Toroid.cc @@ -0,0 +1,34 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file orange/surf/Toroid.cc +//---------------------------------------------------------------------------// +#include "Toroid.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/** + * Construct toroid from origin point and radii + * + * \param origin 3d origin of the toroid. + * \param major_radius Radius from origin to the center of revolved ellipse. + * \param ellipse_xy_radius Radius of ellipse in xy plane, aka 'a'. + * \param ellipse_z_radius Radius of ellipse aligned with z axis, aka 'b'. + */ +Toroid::Toroid(Real3 const& origin, + real_type major_radius, + real_type ellipse_xy_radius, + real_type ellipse_z_radius) + : origin_(origin) + , r_(major_radius) + , a_(ellipse_xy_radius) + , b_(ellipse_z_radius) +{ + CELER_EXPECT(major_radius > 0); + CELER_EXPECT(ellipse_xy_radius > 0); + CELER_EXPECT(ellipse_z_radius > 0); + CELER_EXPECT(major_radius > ellipse_xy_radius); // Degenerate torii +} +} // namespace celeritas From 9eb46499a6a48bdca7eca858426aae35d8f3fb96 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sat, 13 Dec 2025 01:08:37 -0600 Subject: [PATCH 04/14] Added data span getter and span constructor --- src/orange/surf/Toroid.hh | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh index d892109c9a..25c066a0a7 100644 --- a/src/orange/surf/Toroid.hh +++ b/src/orange/surf/Toroid.hh @@ -41,8 +41,8 @@ class Toroid //@{ //! \name Type aliases using Intersections = Array; - using StorageSpan = Span; - using Real3 = Array; + using StorageSpan = Span; + using Real3 = Array; //@} public: @@ -86,6 +86,9 @@ class Toroid //! Radius of revolved ellipse along z axis CELER_FUNCTION real_type ellipse_z_radius() const { return b_; } + //! View of data for type-deleted storage + CELER_FUNCTION StorageSpan data() const { return {&origin_[0], 6}; } + //// CALCULATION /// // Determine the sense of the position relative to this surface @@ -108,4 +111,17 @@ class Toroid real_type a_; // Horizontal radius of revolved ellipse (along xy plane) real_type b_; // Vertical radius of revolved ellipse (along z axis) }; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct from raw data. + */ +template +CELER_FUNCTION Toroid::Toroid(Span data) + : origin_{data[0], data[1], data[2]}, , r_{data[3]}, a_{data[4]}, b_{data[5]} +{ +} + } // namespace celeritas From 7629f3aa415b325f2dedc693c0734ca518d8ad28 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sat, 13 Dec 2025 01:09:48 -0600 Subject: [PATCH 05/14] Started adding test file for Toroid (end of day commit, not finished yet) --- test/orange/surf/Toroid.test.cc | 104 ++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 test/orange/surf/Toroid.test.cc diff --git a/test/orange/surf/Toroid.test.cc b/test/orange/surf/Toroid.test.cc new file mode 100644 index 0000000000..fab3057317 --- /dev/null +++ b/test/orange/surf/Toroid.test.cc @@ -0,0 +1,104 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file orange/surf/Toroid.test.cc +//---------------------------------------------------------------------------// +#include "orange/surf/Toroid.hh" + +#include "corecel/cont/Array.hh" + +#include "SurfaceTestUtils.hh" +#include "celeritas_test.hh" +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// + +using Real3 = Toroid::Real3; +/*! + * Test constructors and span output of toroid class + */ +TEST(ToroidTest, construction) +{ + // Position at 1, 2, 3, major rad 10, xy rad 4, z rad 5 + auto check_props = + [](Toroid const& tor) { + EXPECT_VEC_EQ({1, 2, 3}, tor.origin()); + EXPECT_EQ(10, tor.major_radius()); + EXPECT_EQ(4, tor.ellipse_xy_radius()); + EXPECT_EQ(5, tor.ellipse_z_radius()); + } + + Toroid tor{{1, 2, 3}, 10, 4, 5}; + check_props(tor); + + { + // Reconstruction from data + SCOPED_TRACE("reconstructed"); + Toroid recon{tor.data()}; + check_props(tor); + } +} + +Real3 add(Real3 const a, Real3 const b) +{ + return {a[0] + b[0], a[1] + b[1], a[2] + b[2]}; +} +/*! + * Test sense calculation + */ +TEST(ToroidTest, sense) +{ + Real3 origin{1, 2, 3} Toroid tor{origin, 5, 1, 2}; + Real3 inner_points[] = {{5, 0, 0}, {0, 5, 0}, {5 * 0.707, 5 * 0.707, 1.9}}; + for (Real3 const& point : inner_points) + { + SCOPED_TRACE("Inner point: " + point.to_string()); + EXPECT_EQ(SignedSense::inside, tor.calc_sense(add(point, origin))); + } + + Real3 outer_points[] = {{0, 0, 0}, + {0, 3.9, 0}, + {3.9, 0, 0}, + {-3.9, 0, 0}, + {5, 0, 1.1}, + {6.1, 0, 0}}; + for (Real3 const& point : outer_points) + { + SCOPED_TRACE("Outer point: " + point.to_string()); + EXPECT_EQ(SignedSense::outside, tor.calc_sense(add(point, origin))); + } + + Real3 edge_points[] = {{5.0, 0, 1.0}, {4.0, 0, 0}, {6.0, 0, 0}}; + for (Real3 const& point : edge_points) + { + SCOPED_TRACE("Edge point: " + point.to_string()); + EXPECT_EQ(SignedSense::on, tor.calc_sense(add(point, origin))); + } +} + +/*! + * Test normal vector calculation + */ +TEST(ToroidTest, normal) +{ + Real3 origin{1, 2, 3}; + Toroid tor{origin, 5, 1, 2}; + EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {5.0, 0, 1.0})), + {0, 0, 1.0}); + EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {5.0, 0, -1.0})), + {0, 0, -1.0}); + EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {6.0, 0, 0})), {1.0, 0, 0}); + EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {4.0, 0, 0})), {-1.0, 0, 0}); + EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {0, 6.0, 0})), {0, 1.0, 0}); + EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {0, 4.0, 0})), {0, -1.0, 0}); +} + +/*! + * Test intersection calculation + */ +TEST(ToroidTest, intersect) {} +} // namespace test +} // namespace celeritas From 68f895065f077d510c450d4c78543280a4983432 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sat, 13 Dec 2025 22:38:45 -0600 Subject: [PATCH 06/14] Added ray intersection test cases for toroid --- test/orange/surf/Toroid.test.cc | 71 ++++++++++++++++++++++++++++++++- 1 file changed, 69 insertions(+), 2 deletions(-) diff --git a/test/orange/surf/Toroid.test.cc b/test/orange/surf/Toroid.test.cc index fab3057317..fafb1f3196 100644 --- a/test/orange/surf/Toroid.test.cc +++ b/test/orange/surf/Toroid.test.cc @@ -17,6 +17,34 @@ namespace test //---------------------------------------------------------------------------// using Real3 = Toroid::Real3; +using Intersections = Toroid::Intersections; +//---------------------------------------------------------------------------// +/*! + * Fills a list of fewer than 4 roots with "no real positive root" + */ +Intersections make_roots(std::initializer_list const& inp) +{ + CELER_EXPECT(inp.size() <= Roots{}.size()); + Roots result; + auto iter = std::copy(inp.begin(), inp.end(), result.begin()); + std::fill(iter, result.end(), NumericLimits::infinity()); + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Sorts a given array of four roots and returns the array. + */ +Intersections sorted(Intersections four_roots) +{ + sort(four_roots.begin(), four_roots.end()); + return four_roots; +} + +//---------------------------------------------------------------------------// +// TEST CASES +//---------------------------------------------------------------------------// + /*! * Test constructors and span output of toroid class */ @@ -51,7 +79,8 @@ Real3 add(Real3 const a, Real3 const b) */ TEST(ToroidTest, sense) { - Real3 origin{1, 2, 3} Toroid tor{origin, 5, 1, 2}; + Real3 origin{1, 2, 3}; + Toroid tor{origin, 5, 1, 2}; Real3 inner_points[] = {{5, 0, 0}, {0, 5, 0}, {5 * 0.707, 5 * 0.707, 1.9}}; for (Real3 const& point : inner_points) { @@ -99,6 +128,44 @@ TEST(ToroidTest, normal) /*! * Test intersection calculation */ -TEST(ToroidTest, intersect) {} +TEST(ToroidTest, intersect) +{ + Real3 origin{1, 2, 3}; + Toroid tor{origin, 5, 1, 2}; + + using on = SurfaceState::on; + using off = SurfaceState::off; + + // Ray through center shouldn't hit + Real3 s{add(origin, {0, 0, 2})}; + Real3 u{0, 0, -1}; + EXPECT_VEC_SOFT_EQ(tor.calc_intersections(s, u, off), make_roots({})); + + // Ray inside and out from center should hit exactly once + s = {0, 5, 0}; + u = {0, 1, 0}; + EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), + make_roots({1})); + + // Ray inside towards center should hit 3 times + s = {0, 5, 0}; + u = {0, -1, 0}; + EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), + make_roots({1, 9, 11})); + + // Ray inside towards center again, to have one test that's not nice even + // numbers + s = {0.2, 5.1, 0.1}; + u = {-0.039178047638066676, -0.9990402147707002, 0.019589023819033338}; + Intersections expected = make_roots( + {0.98858498636089176, 10.967709087577969, 9.0295095778886321}); + EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), expected); + + // Ray above torus shouldn't hit torus below it + s = {0.2, 5.1, 2.1}; + u = {1 / 9, 4 / 9, 8 / 9}; + EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), + make_roots({})); +} } // namespace test } // namespace celeritas From 48548c6625d169ef02484618ae6286ee2bfd00d9 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sun, 14 Dec 2025 01:17:23 -0600 Subject: [PATCH 07/14] Unintegrated toroid for now from type lists --- src/orange/OrangeTypes.hh | 1 - src/orange/surf/SurfaceFwd.hh | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/src/orange/OrangeTypes.hh b/src/orange/OrangeTypes.hh index bdf28e6e7c..e43efd390f 100644 --- a/src/orange/OrangeTypes.hh +++ b/src/orange/OrangeTypes.hh @@ -145,7 +145,6 @@ enum class SurfaceType : unsigned char sq, //!< Simple quadric gq, //!< General quadric inv, //!< Involute - tor, //!< Toroid size_ //!< Sentinel value for number of surface types }; diff --git a/src/orange/surf/SurfaceFwd.hh b/src/orange/surf/SurfaceFwd.hh index fcad152c03..b6130e4601 100644 --- a/src/orange/surf/SurfaceFwd.hh +++ b/src/orange/surf/SurfaceFwd.hh @@ -26,6 +26,7 @@ class PlaneAligned; class SimpleQuadric; class Sphere; class SphereCentered; +class Toroid; //---------------------------------------------------------------------------// } // namespace celeritas From efc49d37423ffcef9d7c69c715d4807c21e19c4c Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sun, 14 Dec 2025 01:22:05 -0600 Subject: [PATCH 08/14] Integrated toroid tests, added sense calculation --- src/orange/surf/Toroid.hh | 71 +++++++++++++- test/orange/CMakeLists.txt | 1 + test/orange/surf/Toroid.test.cc | 159 ++++++++++++++++---------------- 3 files changed, 150 insertions(+), 81 deletions(-) diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh index 25c066a0a7..9f63569e76 100644 --- a/src/orange/surf/Toroid.hh +++ b/src/orange/surf/Toroid.hh @@ -6,8 +6,11 @@ //---------------------------------------------------------------------------// #pragma once +#include + #include "corecel/cont/Array.hh" #include "corecel/cont/Span.hh" +#include "corecel/math/Algorithms.hh" #include "orange/OrangeTypes.hh" namespace celeritas @@ -25,6 +28,7 @@ namespace celeritas * / b \ / \ * | | | | | * |-a--+ | o-----r--+ | + * | | | | * \ / \ / * \ / \ / * ⁻⁻⁻ ⁻⁻⁻⁻⁻⁻⁻⁻⁻ ⁻⁻⁻ @@ -40,9 +44,9 @@ class Toroid public: //@{ //! \name Type aliases - using Intersections = Array; + using Intersections = Array; using StorageSpan = Span; - using Real3 = Array; + using Real3 = Array; //@} public: @@ -51,7 +55,7 @@ class Toroid // Surface type identifier static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type() { - return SurfaceType::tor; + return SurfaceType::s; } //! Safety distance is calculable w/xy of normal and ellipse safety @@ -120,8 +124,67 @@ class Toroid */ template CELER_FUNCTION Toroid::Toroid(Span data) - : origin_{data[0], data[1], data[2]}, , r_{data[3]}, a_{data[4]}, b_{data[5]} + : origin_{data[0], data[1], data[2]}, r_{data[3]}, a_{data[4]}, b_{data[5]} +{ +} + +//---------------------------------------------------------------------------// +/** + * Determine the sense of the position relative to this surface. + * + * For a toroid, being inside the toroid (i) counts as inside, outside + * (including in the 'hole' region) (o) as outside, and on the surface exactly + * as on (s). + * ___ _________ ___ + * / \ / \ + * / \ o / \ + * | | | | o + * | | | i s + * \ / \ / + * \ / \ / + * ⁻⁻⁻ ⁻⁻⁻⁻⁻⁻⁻⁻⁻ ⁻⁻⁻ + */ +CELER_FUNCTION SignedSense Toroid::calc_sense(Real3 const& pos) const +{ + auto [x, y, z] = pos; + real_type x0 = x - origin_[0]; + real_type y0 = y - origin_[1]; + real_type z0 = z - origin_[2]; + + real_type val = (ipow<2>(ipow<2>(x0) + ipow<2>(y0) + ipow<2>(z0 * a_ / b_) + + (ipow<2>(r_) - ipow<2>(a_))) + - (4 * ipow<2>(r_)) * (ipow<2>(x0) + ipow<2>(y0))); + if (val < 0) + return SignedSense::inside; + else if (val > 0) + return SignedSense::outside; + else + return SignedSense::on; +} + +//---------------------------------------------------------------------------// +/** + * Calculate all possible straight-line intersections between the given ray and + * this surface. + */ +CELER_FUNCTION auto Toroid::calc_intersections(Real3 const& pos, + Real3 const& dir, + SurfaceState on_surface) const + -> Intersections +{ + return Intersections{no_intersection(), + no_intersection(), + no_intersection(), + no_intersection()}; +} + +//---------------------------------------------------------------------------// +/** + * Calculate outward facing normal at a position on or close to the surface. + */ +CELER_FUNCTION auto Toroid::calc_normal(Real3 const& pos) const -> Real3 { + return Real3{0, 0, 0}; } } // namespace celeritas diff --git a/test/orange/CMakeLists.txt b/test/orange/CMakeLists.txt index 61c61ffa94..457818b236 100644 --- a/test/orange/CMakeLists.txt +++ b/test/orange/CMakeLists.txt @@ -110,6 +110,7 @@ celeritas_add_test(surf/SimpleQuadric.test.cc) celeritas_add_test(surf/Sphere.test.cc) celeritas_add_test(surf/SphereCentered.test.cc) celeritas_add_test(surf/Involute.test.cc) +celeritas_add_test(surf/Toroid.test.cc) # Surface utilities celeritas_add_test(surf/FaceNamer.test.cc) diff --git a/test/orange/surf/Toroid.test.cc b/test/orange/surf/Toroid.test.cc index fafb1f3196..8fd6e586c4 100644 --- a/test/orange/surf/Toroid.test.cc +++ b/test/orange/surf/Toroid.test.cc @@ -4,9 +4,11 @@ //---------------------------------------------------------------------------// //! \file orange/surf/Toroid.test.cc //---------------------------------------------------------------------------// -#include "orange/surf/Toroid.hh" +#include "orange/surf/Toroid.cc" #include "corecel/cont/Array.hh" +#include "corecel/cont/ArrayIO.hh" +#include "orange/surf/Toroid.hh" #include "SurfaceTestUtils.hh" #include "celeritas_test.hh" @@ -22,10 +24,10 @@ using Intersections = Toroid::Intersections; /*! * Fills a list of fewer than 4 roots with "no real positive root" */ -Intersections make_roots(std::initializer_list const& inp) +Intersections make_inters(std::initializer_list const& inp) { - CELER_EXPECT(inp.size() <= Roots{}.size()); - Roots result; + CELER_EXPECT(inp.size() <= Intersections{}.size()); + Intersections result; auto iter = std::copy(inp.begin(), inp.end(), result.begin()); std::fill(iter, result.end(), NumericLimits::infinity()); return result; @@ -51,15 +53,16 @@ Intersections sorted(Intersections four_roots) TEST(ToroidTest, construction) { // Position at 1, 2, 3, major rad 10, xy rad 4, z rad 5 - auto check_props = - [](Toroid const& tor) { - EXPECT_VEC_EQ({1, 2, 3}, tor.origin()); - EXPECT_EQ(10, tor.major_radius()); - EXPECT_EQ(4, tor.ellipse_xy_radius()); - EXPECT_EQ(5, tor.ellipse_z_radius()); - } - - Toroid tor{{1, 2, 3}, 10, 4, 5}; + auto check_props = [](Toroid const& tor) { + Real3 expect{1, 2, 3}; + Real3 actual = tor.origin(); + EXPECT_VEC_EQ(expect, actual); + EXPECT_EQ(10, tor.major_radius()); + EXPECT_EQ(4, tor.ellipse_xy_radius()); + EXPECT_EQ(5, tor.ellipse_z_radius()); + }; + + Toroid tor{Real3{1, 2, 3}, 10, 4, 5}; check_props(tor); { @@ -84,7 +87,7 @@ TEST(ToroidTest, sense) Real3 inner_points[] = {{5, 0, 0}, {0, 5, 0}, {5 * 0.707, 5 * 0.707, 1.9}}; for (Real3 const& point : inner_points) { - SCOPED_TRACE("Inner point: " + point.to_string()); + SCOPED_TRACE("Inner point: " + to_string(point)); EXPECT_EQ(SignedSense::inside, tor.calc_sense(add(point, origin))); } @@ -92,80 +95,82 @@ TEST(ToroidTest, sense) {0, 3.9, 0}, {3.9, 0, 0}, {-3.9, 0, 0}, - {5, 0, 1.1}, + {5, 0, 2.1}, {6.1, 0, 0}}; for (Real3 const& point : outer_points) { - SCOPED_TRACE("Outer point: " + point.to_string()); + SCOPED_TRACE("Outer point: " + to_string(point)); EXPECT_EQ(SignedSense::outside, tor.calc_sense(add(point, origin))); } - Real3 edge_points[] = {{5.0, 0, 1.0}, {4.0, 0, 0}, {6.0, 0, 0}}; + Real3 edge_points[] = {{5.0, 0, 2.0}, {4.0, 0, 0}, {6.0, 0, 0}}; for (Real3 const& point : edge_points) { - SCOPED_TRACE("Edge point: " + point.to_string()); + SCOPED_TRACE("Edge point: " + to_string(point)); EXPECT_EQ(SignedSense::on, tor.calc_sense(add(point, origin))); } } -/*! - * Test normal vector calculation - */ -TEST(ToroidTest, normal) -{ - Real3 origin{1, 2, 3}; - Toroid tor{origin, 5, 1, 2}; - EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {5.0, 0, 1.0})), - {0, 0, 1.0}); - EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {5.0, 0, -1.0})), - {0, 0, -1.0}); - EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {6.0, 0, 0})), {1.0, 0, 0}); - EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {4.0, 0, 0})), {-1.0, 0, 0}); - EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {0, 6.0, 0})), {0, 1.0, 0}); - EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {0, 4.0, 0})), {0, -1.0, 0}); -} - -/*! - * Test intersection calculation - */ -TEST(ToroidTest, intersect) -{ - Real3 origin{1, 2, 3}; - Toroid tor{origin, 5, 1, 2}; - - using on = SurfaceState::on; - using off = SurfaceState::off; - - // Ray through center shouldn't hit - Real3 s{add(origin, {0, 0, 2})}; - Real3 u{0, 0, -1}; - EXPECT_VEC_SOFT_EQ(tor.calc_intersections(s, u, off), make_roots({})); - - // Ray inside and out from center should hit exactly once - s = {0, 5, 0}; - u = {0, 1, 0}; - EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), - make_roots({1})); - - // Ray inside towards center should hit 3 times - s = {0, 5, 0}; - u = {0, -1, 0}; - EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), - make_roots({1, 9, 11})); - - // Ray inside towards center again, to have one test that's not nice even - // numbers - s = {0.2, 5.1, 0.1}; - u = {-0.039178047638066676, -0.9990402147707002, 0.019589023819033338}; - Intersections expected = make_roots( - {0.98858498636089176, 10.967709087577969, 9.0295095778886321}); - EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), expected); - - // Ray above torus shouldn't hit torus below it - s = {0.2, 5.1, 2.1}; - u = {1 / 9, 4 / 9, 8 / 9}; - EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), - make_roots({})); -} +// /*! +// * Test normal vector calculation +// */ +// TEST(ToroidTest, normal) +// { +// Real3 origin{1, 2, 3}; +// Toroid tor{origin, 5, 1, 2}; +// EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {5.0, 0, 1.0})), +// {0, 0, 1.0}); +// EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {5.0, 0, -1.0})), +// {0, 0, -1.0}); +// EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {6.0, 0, 0})), {1.0, 0, +// 0}); EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {4.0, 0, 0})), +// {-1.0, 0, 0}); EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {0, 6.0, +// 0})), {0, 1.0, 0}); EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, +// {0, 4.0, 0})), {0, -1.0, 0}); +// } + +// /*! +// * Test intersection calculation +// */ +// TEST(ToroidTest, intersect) +// { +// Real3 origin{1, 2, 3}; +// Toroid tor{origin, 5, 1, 2}; + +// using on = SurfaceState::on; +// using off = SurfaceState::off; + +// // Ray through center shouldn't hit +// Real3 s{add(origin, {0, 0, 2})}; +// Real3 u{0, 0, -1}; +// EXPECT_VEC_SOFT_EQ(tor.calc_intersections(s, u, off), make_inters({})); + +// // Ray inside and out from center should hit exactly once +// s = {0, 5, 0}; +// u = {0, 1, 0}; +// EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), +// make_inters({1})); + +// // Ray inside towards center should hit 3 times +// s = {0, 5, 0}; +// u = {0, -1, 0}; +// EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), +// make_inters({1, 9, 11})); + +// // Ray inside towards center again, to have one test that's not nice +// even +// // numbers +// s = {0.2, 5.1, 0.1}; +// u = {-0.039178047638066676, -0.9990402147707002, 0.019589023819033338}; +// Intersections expected = make_inters( +// {0.98858498636089176, 10.967709087577969, 9.0295095778886321}); +// EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), expected); + +// // Ray above torus shouldn't hit torus below it +// s = {0.2, 5.1, 2.1}; +// u = {1 / 9, 4 / 9, 8 / 9}; +// EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), +// make_inters({})); +// } } // namespace test } // namespace celeritas From 762e40d8341a27ee850994b29a21d5035d2cef24 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sun, 14 Dec 2025 17:06:47 -0600 Subject: [PATCH 09/14] Added ray intersection algorithm to toroid class --- src/orange/surf/Toroid.hh | 124 +++++++++++++++++++++++++++---- test/orange/surf/Toroid.test.cc | 125 ++++++++++++++++---------------- 2 files changed, 173 insertions(+), 76 deletions(-) diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh index 9f63569e76..6d0e0cf776 100644 --- a/src/orange/surf/Toroid.hh +++ b/src/orange/surf/Toroid.hh @@ -6,11 +6,15 @@ //---------------------------------------------------------------------------// #pragma once +#include #include #include "corecel/cont/Array.hh" +#include "corecel/cont/ArrayIO.hh" #include "corecel/cont/Span.hh" #include "corecel/math/Algorithms.hh" +#include "corecel/math/ArrayUtils.hh" +#include "corecel/math/FerrariSolver.hh" #include "orange/OrangeTypes.hh" namespace celeritas @@ -47,6 +51,8 @@ class Toroid using Intersections = Array; using StorageSpan = Span; using Real3 = Array; + using Real4 = Array; + using Real5 = Array; //@} public: @@ -106,6 +112,7 @@ class Toroid inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const; private: + //// DATA //// // Location of center of toroid Real3 origin_; @@ -114,6 +121,17 @@ class Toroid // xy plane) real_type a_; // Horizontal radius of revolved ellipse (along xy plane) real_type b_; // Vertical radius of revolved ellipse (along z axis) + + //// HELPER FUNCTIONS //// + // Calculate the coefficients of the polynomial for ray intersection + inline CELER_FUNCTION Real5 calc_intersection_polynomial( + Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const; + + // Shorthand to subtract b from a + static CELER_FUNCTION Real3 sub(Real3 const& a, Real3 const& b); + + // Shorthnad to square a number + static CELER_FUNCTION real_type sq(real_type val); }; //---------------------------------------------------------------------------// @@ -146,14 +164,10 @@ CELER_FUNCTION Toroid::Toroid(Span data) */ CELER_FUNCTION SignedSense Toroid::calc_sense(Real3 const& pos) const { - auto [x, y, z] = pos; - real_type x0 = x - origin_[0]; - real_type y0 = y - origin_[1]; - real_type z0 = z - origin_[2]; - - real_type val = (ipow<2>(ipow<2>(x0) + ipow<2>(y0) + ipow<2>(z0 * a_ / b_) - + (ipow<2>(r_) - ipow<2>(a_))) - - (4 * ipow<2>(r_)) * (ipow<2>(x0) + ipow<2>(y0))); + auto [x0, y0, z0] = sub(pos, origin_); + + real_type val = (sq(sq(x0) + sq(y0) + sq(z0 * a_ / b_) + (sq(r_) - sq(a_))) + - (4 * sq(r_)) * (sq(x0) + sq(y0))); if (val < 0) return SignedSense::inside; else if (val > 0) @@ -163,7 +177,7 @@ CELER_FUNCTION SignedSense Toroid::calc_sense(Real3 const& pos) const } //---------------------------------------------------------------------------// -/** +/*! * Calculate all possible straight-line intersections between the given ray and * this surface. */ @@ -172,19 +186,99 @@ CELER_FUNCTION auto Toroid::calc_intersections(Real3 const& pos, SurfaceState on_surface) const -> Intersections { - return Intersections{no_intersection(), - no_intersection(), - no_intersection(), - no_intersection()}; + Real5 abcde = calc_intersection_polynomial(pos, dir, on_surface); + std::cout << "Polynomial: " << to_string(abcde) << "\n"; + FerrariSolver solve{}; // Default tolerance + Intersections roots; + + if (on_surface == SurfaceState::on) + { + auto [a, b, c, d, e] = abcde; + roots = solve(Real4{a, b, c, d}); + } + else + { + roots = solve(abcde); + } + + return roots; } //---------------------------------------------------------------------------// -/** +/*! + * Calculate the coefficients of the polynomial corresponding to the given + * ray's intersections with the toroid. + * + * Written referencing Graphics Gems II. + */ +CELER_FUNCTION auto +Toroid::calc_intersection_polynomial(Real3 const& pos, + Real3 const& dir, + SurfaceState on_surface) const -> Real5 +{ + auto [x0, y0, z0] = sub(pos, origin_); + auto [ax, ay, az] = make_unit_vector(dir); + std::cout << "pos: " << to_string(pos) + << "dir: " << to_string(make_unit_vector(dir)) << "\n"; + + // Intermediate terms + real_type p = sq(a_) / sq(b_); + real_type A0 = 4 * sq(r_); + real_type B0 = sq(r_) - sq(a_); + std::cout << "p: " << p << ", A0: " << A0 << ", B0: " << B0 << "\n"; + + real_type f = 1 - sq(az); + real_type g = f + p * sq(az); + real_type l = 2 * (x0 * ax + y0 * ay); + real_type t = sq(x0) + sq(y0); + real_type q = A0 / sq(g); + real_type m = (l + 2 * p * z0 * az) / g; + real_type u = (t + p * sq(z0) + B0) / g; + + std::cout << "f: " << f << ", g: " << g << ", l: " << l << ", t: " << t + << ", q: " << q << ", m: " << m << ", u: " << u << "\n"; + + // Polynomial coefficients, i.e. cn*x^n + real_type c4 = 1; + real_type c3 = 2 * m; + real_type c2 = sq(m) + 2 * u - q * f; + real_type c1 = 2 * m * u - q * l; + real_type c0 = on_surface == SurfaceState::on ? 0 : sq(u) - q * t; + // Potential refinement of c0 if close to 0? + + return Real5{c4, c3, c2, c1, c0}; +} + +//---------------------------------------------------------------------------// +/*! * Calculate outward facing normal at a position on or close to the surface. */ CELER_FUNCTION auto Toroid::calc_normal(Real3 const& pos) const -> Real3 { - return Real3{0, 0, 0}; + auto [x0, y0, z0] = sub(pos, origin_); + + real_type d = std::sqrt(sq(x0) + sq(y0)); + real_type f = 2 * (d - r_) / (d * sq(a_)); + Real3 n{x0 * f, y0 * f, (2 * z0) / (sq(b_))}; + return make_unit_vector(n); +} + +//---------------------------------------------------------------------------// +/*! + * Shorthand for power macro for readability, as squares are used frequently + */ +CELER_FUNCTION real_type Toroid::sq(real_type val) +{ + return ipow<2>(val); +} + +//---------------------------------------------------------------------------// +/*! + * Shorthand for subtracting vector b from vector a + */ +CELER_FUNCTION Real3 Toroid::sub(Real3 const& a, Real3 const& b) +{ + return Real3{a[0] - b[0], a[1] - b[1], a[2] - b[2]}; } } // namespace celeritas diff --git a/test/orange/surf/Toroid.test.cc b/test/orange/surf/Toroid.test.cc index 8fd6e586c4..64f30a20c2 100644 --- a/test/orange/surf/Toroid.test.cc +++ b/test/orange/surf/Toroid.test.cc @@ -111,66 +111,69 @@ TEST(ToroidTest, sense) } } -// /*! -// * Test normal vector calculation -// */ -// TEST(ToroidTest, normal) -// { -// Real3 origin{1, 2, 3}; -// Toroid tor{origin, 5, 1, 2}; -// EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {5.0, 0, 1.0})), -// {0, 0, 1.0}); -// EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {5.0, 0, -1.0})), -// {0, 0, -1.0}); -// EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {6.0, 0, 0})), {1.0, 0, -// 0}); EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {4.0, 0, 0})), -// {-1.0, 0, 0}); EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, {0, 6.0, -// 0})), {0, 1.0, 0}); EXPECT_VEC_SOFT_EQ(tor.calc_normal(add(origin, -// {0, 4.0, 0})), {0, -1.0, 0}); -// } - -// /*! -// * Test intersection calculation -// */ -// TEST(ToroidTest, intersect) -// { -// Real3 origin{1, 2, 3}; -// Toroid tor{origin, 5, 1, 2}; - -// using on = SurfaceState::on; -// using off = SurfaceState::off; - -// // Ray through center shouldn't hit -// Real3 s{add(origin, {0, 0, 2})}; -// Real3 u{0, 0, -1}; -// EXPECT_VEC_SOFT_EQ(tor.calc_intersections(s, u, off), make_inters({})); - -// // Ray inside and out from center should hit exactly once -// s = {0, 5, 0}; -// u = {0, 1, 0}; -// EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), -// make_inters({1})); - -// // Ray inside towards center should hit 3 times -// s = {0, 5, 0}; -// u = {0, -1, 0}; -// EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), -// make_inters({1, 9, 11})); - -// // Ray inside towards center again, to have one test that's not nice -// even -// // numbers -// s = {0.2, 5.1, 0.1}; -// u = {-0.039178047638066676, -0.9990402147707002, 0.019589023819033338}; -// Intersections expected = make_inters( -// {0.98858498636089176, 10.967709087577969, 9.0295095778886321}); -// EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), expected); - -// // Ray above torus shouldn't hit torus below it -// s = {0.2, 5.1, 2.1}; -// u = {1 / 9, 4 / 9, 8 / 9}; -// EXPECT_VEC_SOFT_EQ(sorted(tor.calc_intersections(s, u, off)), -// make_inters({})); -// } +/*! + * Test normal vector calculation + */ +TEST(ToroidTest, normal) +{ + Real3 origin{1, 2, 3}; + Toroid tor{origin, 5, 1, 2}; + EXPECT_VEC_SOFT_EQ((Real3{0, 0, 1}), + tor.calc_normal(add(origin, {5, 0, 2}))); + EXPECT_VEC_SOFT_EQ((Real3{0, 0, -1}), + tor.calc_normal(add(origin, {5, 0, -2}))); + EXPECT_VEC_SOFT_EQ((Real3{1, 0, 0}), + tor.calc_normal(add(origin, {6, 0, 0}))); + EXPECT_VEC_SOFT_EQ((Real3{-1, 0, 0}), + tor.calc_normal(add(origin, {4, 0, 0}))); + EXPECT_VEC_SOFT_EQ((Real3{0, 1, 0}), + tor.calc_normal(add(origin, {0, 6, 0}))); +} + +/*! + * Test intersection calculation + */ +TEST(ToroidTest, intersect) +{ + Real3 origin{1, 2, 3}; + Toroid tor{origin, 5, 1, 2}; + + SurfaceState on = SurfaceState::on; + SurfaceState off = SurfaceState::off; + + // Ray through center shouldn't hit + Real3 s{add(origin, {0, 0, 2})}; + Real3 u{0, 0, -1}; + EXPECT_VEC_SOFT_EQ(make_inters({}), + sorted(tor.calc_intersections(s, u, off))); + + // Ray inside and out from center should hit exactly once + s = add(origin, {0, 5, 0}); + u = {0, 1, 0}; + EXPECT_VEC_SOFT_EQ(make_inters({1}), + sorted(tor.calc_intersections(s, u, off))); + + // Ray inside towards center should hit 3 times + s = add(origin, {0, 5, 0}); + u = {0, -1, 0}; + EXPECT_VEC_SOFT_EQ(make_inters({1, 9, 11}), + sorted(tor.calc_intersections(s, u, off))); + + // Ray inside towards center again, to have one test that's not nice even + // numbers + s = add(origin, {0.2, 5.1, 0.1}); + u = {-0.039178047638066676, -0.9990402147707002, 0.019589023819033338}; + Intersections expected = sorted(make_inters( + {1.1022820700552722, 11.093380864669665, 9.1154137268987121})); + EXPECT_VEC_NEAR(expected, + sorted(tor.calc_intersections(s, u, off)), + 1e-4); // Not a precision test, that comes later + + // Ray above torus shouldn't hit torus below it + s = add(origin, {0.2, 5.1, 2.1}); + u = {1.0 / 9, 4.0 / 9, 8.0 / 9}; + EXPECT_VEC_SOFT_EQ(make_inters({}), + sorted(tor.calc_intersections(s, u, off))); +} } // namespace test } // namespace celeritas From 4f86208e0824de7c794cc985629e4fe3136448a5 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sun, 14 Dec 2025 17:13:56 -0600 Subject: [PATCH 10/14] Removed print statements from debugging ray intersection --- src/orange/surf/Toroid.hh | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh index 6d0e0cf776..27cf2bc138 100644 --- a/src/orange/surf/Toroid.hh +++ b/src/orange/surf/Toroid.hh @@ -7,7 +7,6 @@ #pragma once #include -#include #include "corecel/cont/Array.hh" #include "corecel/cont/ArrayIO.hh" @@ -187,7 +186,6 @@ CELER_FUNCTION auto Toroid::calc_intersections(Real3 const& pos, -> Intersections { Real5 abcde = calc_intersection_polynomial(pos, dir, on_surface); - std::cout << "Polynomial: " << to_string(abcde) << "\n"; FerrariSolver solve{}; // Default tolerance Intersections roots; @@ -218,14 +216,11 @@ Toroid::calc_intersection_polynomial(Real3 const& pos, { auto [x0, y0, z0] = sub(pos, origin_); auto [ax, ay, az] = make_unit_vector(dir); - std::cout << "pos: " << to_string(pos) - << "dir: " << to_string(make_unit_vector(dir)) << "\n"; // Intermediate terms real_type p = sq(a_) / sq(b_); real_type A0 = 4 * sq(r_); real_type B0 = sq(r_) - sq(a_); - std::cout << "p: " << p << ", A0: " << A0 << ", B0: " << B0 << "\n"; real_type f = 1 - sq(az); real_type g = f + p * sq(az); @@ -235,9 +230,6 @@ Toroid::calc_intersection_polynomial(Real3 const& pos, real_type m = (l + 2 * p * z0 * az) / g; real_type u = (t + p * sq(z0) + B0) / g; - std::cout << "f: " << f << ", g: " << g << ", l: " << l << ", t: " << t - << ", q: " << q << ", m: " << m << ", u: " << u << "\n"; - // Polynomial coefficients, i.e. cn*x^n real_type c4 = 1; real_type c3 = 2 * m; From 4e272097a2a7c547432210a028ec7c2d23f570f4 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sun, 14 Dec 2025 17:29:03 -0600 Subject: [PATCH 11/14] Fixed typo in storage span declaraiton (said 6 instead of correct 4) --- src/orange/surf/Toroid.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh index 27cf2bc138..d78c475356 100644 --- a/src/orange/surf/Toroid.hh +++ b/src/orange/surf/Toroid.hh @@ -96,7 +96,7 @@ class Toroid CELER_FUNCTION real_type ellipse_z_radius() const { return b_; } //! View of data for type-deleted storage - CELER_FUNCTION StorageSpan data() const { return {&origin_[0], 6}; } + CELER_FUNCTION StorageSpan data() const { return {&origin_[0], 4}; } //// CALCULATION /// From 310b59ef71266a83660fdca391785829a3b8dd49 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sun, 14 Dec 2025 18:20:03 -0600 Subject: [PATCH 12/14] Added Toroid IO functions --- src/orange/surf/SurfaceIO.cc | 9 +++++++++ src/orange/surf/SurfaceIO.hh | 2 ++ 2 files changed, 11 insertions(+) diff --git a/src/orange/surf/SurfaceIO.cc b/src/orange/surf/SurfaceIO.cc index 7a0cc78352..01597873cc 100644 --- a/src/orange/surf/SurfaceIO.cc +++ b/src/orange/surf/SurfaceIO.cc @@ -22,6 +22,7 @@ #include "SimpleQuadric.hh" // IWYU pragma: associated #include "Sphere.hh" // IWYU pragma: associated #include "SphereCentered.hh" // IWYU pragma: associated +#include "Toroid.hh" //IWYU pragma: associated namespace celeritas { @@ -124,6 +125,14 @@ std::ostream& operator<<(std::ostream& os, SphereCentered const& s) os << "Sphere: r=" << std::sqrt(s.radius_sq()); return os; } +//---------------------------------------------------------------------------// +std::ostream& operator<<(std::ostream& os, Toroid const& tor) +{ + os << "Toroid: r=" << tor.major_radius() + << ", a=" << tor.ellipse_xy_radius() << ", b=" << tor.ellipse_z_radius() + << " at " << tor.origin(); + return os; +} //---------------------------------------------------------------------------// #undef ORANGE_INSTANTIATE_SHAPE_STREAM diff --git a/src/orange/surf/SurfaceIO.hh b/src/orange/surf/SurfaceIO.hh index 95d77b9cb5..9c9367ecaf 100644 --- a/src/orange/surf/SurfaceIO.hh +++ b/src/orange/surf/SurfaceIO.hh @@ -40,6 +40,8 @@ std::ostream& operator<<(std::ostream&, SimpleQuadric const&); std::ostream& operator<<(std::ostream&, Sphere const&); std::ostream& operator<<(std::ostream&, SphereCentered const&); + +std::ostream& operator<<(std::ostream&, Toroid const&); //!@} //---------------------------------------------------------------------------// } // namespace celeritas From 7bdfd75a13b79b37cf4d0006203a7600e1808fac Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sun, 14 Dec 2025 20:14:36 -0600 Subject: [PATCH 13/14] Added Toroid.cc to CMakeLists, fixed static private functions sq and sub not being declared inline properly --- src/orange/CMakeLists.txt | 1 + src/orange/surf/Toroid.hh | 4 ++-- test/orange/surf/Toroid.test.cc | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/orange/CMakeLists.txt b/src/orange/CMakeLists.txt index 7bd46d75b8..43d479ab64 100644 --- a/src/orange/CMakeLists.txt +++ b/src/orange/CMakeLists.txt @@ -81,6 +81,7 @@ list(APPEND SOURCES surf/SurfaceClipper.cc surf/SurfaceIO.cc surf/SurfaceSimplifier.cc + surf/Toroid.cc surf/VariantSurface.cc surf/detail/SurfaceTransformer.cc surf/detail/SurfaceTranslator.cc diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh index d78c475356..6d44f3aa17 100644 --- a/src/orange/surf/Toroid.hh +++ b/src/orange/surf/Toroid.hh @@ -127,10 +127,10 @@ class Toroid Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const; // Shorthand to subtract b from a - static CELER_FUNCTION Real3 sub(Real3 const& a, Real3 const& b); + static inline CELER_FUNCTION Real3 sub(Real3 const& a, Real3 const& b); // Shorthnad to square a number - static CELER_FUNCTION real_type sq(real_type val); + static inline CELER_FUNCTION real_type sq(real_type val); }; //---------------------------------------------------------------------------// diff --git a/test/orange/surf/Toroid.test.cc b/test/orange/surf/Toroid.test.cc index 64f30a20c2..5fc5a30f4c 100644 --- a/test/orange/surf/Toroid.test.cc +++ b/test/orange/surf/Toroid.test.cc @@ -4,7 +4,7 @@ //---------------------------------------------------------------------------// //! \file orange/surf/Toroid.test.cc //---------------------------------------------------------------------------// -#include "orange/surf/Toroid.cc" +#include "orange/surf/Toroid.hh" #include "corecel/cont/Array.hh" #include "corecel/cont/ArrayIO.hh" From 109dd8408da267d6bcf47259cdf6c9715c28b573 Mon Sep 17 00:00:00 2001 From: Owen Strong Date: Sun, 14 Dec 2025 20:57:07 -0600 Subject: [PATCH 14/14] Made surface type more clearly not implemented yet (since Toroid is not integrated) --- src/orange/surf/Toroid.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh index 6d44f3aa17..b012f4bd0b 100644 --- a/src/orange/surf/Toroid.hh +++ b/src/orange/surf/Toroid.hh @@ -58,9 +58,9 @@ class Toroid //// CLASS ATTRIBUTES //// // Surface type identifier - static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type() + static SurfaceType surface_type() { - return SurfaceType::s; + CELER_NOT_IMPLEMENTED("runtime toroid"); } //! Safety distance is calculable w/xy of normal and ellipse safety