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/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 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 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 diff --git a/src/orange/surf/Toroid.hh b/src/orange/surf/Toroid.hh new file mode 100644 index 0000000000..b012f4bd0b --- /dev/null +++ b/src/orange/surf/Toroid.hh @@ -0,0 +1,276 @@ +//------------------------------- -*- 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 + +#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 +{ +//---------------------------------------------------------------------------// +/*! + * 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]. + */ +class Toroid +{ + public: + //@{ + //! \name Type aliases + using Intersections = Array; + using StorageSpan = Span; + using Real3 = Array; + using Real4 = Array; + using Real5 = Array; + //@} + + public: + //// CLASS ATTRIBUTES //// + + // Surface type identifier + static SurfaceType surface_type() + { + CELER_NOT_IMPLEMENTED("runtime toroid"); + } + + //! 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_; } + + //! View of data for type-deleted storage + CELER_FUNCTION StorageSpan data() const { return {&origin_[0], 4}; } + + //// 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: + //// DATA //// + // 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) + + //// 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 inline CELER_FUNCTION Real3 sub(Real3 const& a, Real3 const& b); + + // Shorthnad to square a number + static inline CELER_FUNCTION real_type sq(real_type val); +}; + +//---------------------------------------------------------------------------// +// 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]} +{ +} + +//---------------------------------------------------------------------------// +/** + * 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 [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) + 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 +{ + Real5 abcde = calc_intersection_polynomial(pos, dir, on_surface); + 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); + + // Intermediate terms + real_type p = sq(a_) / sq(b_); + real_type A0 = 4 * sq(r_); + real_type B0 = sq(r_) - sq(a_); + + 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; + + // 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 +{ + 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/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 new file mode 100644 index 0000000000..5fc5a30f4c --- /dev/null +++ b/test/orange/surf/Toroid.test.cc @@ -0,0 +1,179 @@ +//------------------------------- -*- 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 "corecel/cont/ArrayIO.hh" +#include "orange/surf/Toroid.hh" + +#include "SurfaceTestUtils.hh" +#include "celeritas_test.hh" +namespace celeritas +{ +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_inters(std::initializer_list const& inp) +{ + 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; +} + +//---------------------------------------------------------------------------// +/*! + * 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 + */ +TEST(ToroidTest, construction) +{ + // Position at 1, 2, 3, major rad 10, xy rad 4, z rad 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); + + { + // 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: " + to_string(point)); + 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, 2.1}, + {6.1, 0, 0}}; + for (Real3 const& point : outer_points) + { + SCOPED_TRACE("Outer point: " + to_string(point)); + EXPECT_EQ(SignedSense::outside, tor.calc_sense(add(point, origin))); + } + + 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: " + 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((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