orange/surf/Toroid.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file orange/surf/Toroid.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <cmath>
0010
0011 #include "corecel/cont/Array.hh"
0012 #include "corecel/cont/Span.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/math/ArrayOperators.hh"
0015 #include "corecel/math/ArrayUtils.hh"
0016 #include "corecel/math/FerrariSolver.hh"
0017 #include "orange/OrangeTypes.hh"
0018 #include "orange/SenseUtils.hh"
0019
0020 namespace celeritas
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Z-aligned Elliptical Toroid.
0025  *
0026  * An elliptical toroid is a shape created by revolving an axis-aligned ellipse
0027  * around a central axis. This shape can be used in everything from pipe bends
0028  * to tokamaks in fusion reactors. It possesses a major radius \em r, and
0029  * ellipse radii \em a and \em b, as shown in the below diagram:
0030  *
0031  * \verbatim
0032        ___   _________   ___
0033      /  |  \           /     \
0034     /   b   \         /       \
0035    |    |    |       |         |
0036    |-a--+    |   o-----r--+    |
0037    |         |       |         |
0038     \       /         \       /
0039      \...../...........\...../
0040
0041    \endverbatim
0042  *
0043  * This torus can be defined with the following quartic equation:
0044  * \f[
0045  *   (x^2 + y^2 + p*y^2 + B_0) - A_0 * (x^2 + y^2) = 0
0046  * \f]
0047  * where \f$p = a^2/b^2 \f$,  \f$A_0 = 4 r^2 \f$, and \f$ B_0 = (r^2-a^2)\f$.
0048  *
0049  * The intersection is calculated using an algorithm from
0050  * \cite{arvo-graphicsgems-1995}.
0051  */
0052 class Toroid
0053 {
0054   public:
0055     //@{
0056     //! \name Type aliases
0057     using Intersections = Array<real_type, 4>;
0058     using StorageSpan = Span<real_type const, 6>;
0059     using Real3 = Array<real_type, 3>;
0060     using Real4 = Array<real_type, 4>;
0061     using Real5 = Array<real_type, 5>;
0062     //@}
0063
0064   public:
0065     //// CLASS ATTRIBUTES ////
0066
0067     // Surface type identifier
0068     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type()
0069     {
0070         return SurfaceType::tor;
0071     }
0072
0073     //! Safety distance is calculable w/xy of normal and ellipse safety
0074     //! distance, but this is out of scope at first and might not be trivially
0075     //! calculable
0076     //! https://web.archive.org/web/20170829172516/https://www.spaceroots.org/documents/distance/distance-to-ellipse.pdf
0077     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; }
0078
0079   public:
0080     //// CONSTRUCTORS ////
0081
0082     explicit Toroid(Real3 const& origin,
0083                     real_type major_radius,
0084                     real_type ellipse_xy_radius,
0085                     real_type ellipse_z_radius);
0086
0087     // Construct from raw data
0088     template<class R>
0089     explicit inline CELER_FUNCTION Toroid(Span<R, StorageSpan::extent>);
0090
0091     //// ACCESSORS ////
0092
0093     //! Center of the toroid (in the donut hole)
0094     CELER_FUNCTION Real3 const& origin() const { return origin_; }
0095
0096     //! Radius from origin to center of revolved ellipse
0097     CELER_FUNCTION real_type major_radius() const { return r_; }
0098
0099     //! Radius of revolved ellipse along xy plane
0100     CELER_FUNCTION real_type ellipse_xy_radius() const { return a_; }
0101
0102     //! Radius of revolved ellipse along z axis
0103     CELER_FUNCTION real_type ellipse_z_radius() const { return b_; }
0104
0105     //! View of data for type-deleted storage
0106     CELER_FUNCTION StorageSpan data() const
0107     {
0108         static_assert(sizeof(*this) == sizeof(real_type) * 6);
0109         return {&origin_[0], 6};
0110     }
0111
0112     //// CALCULATION ///
0113
0114     // Determine the sense of the position relative to this surface
0115     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0116
0117     // Calculate all possible straight-line intersections with this surface
0118     inline CELER_FUNCTION Intersections calc_intersections(
0119         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0120
0121     // Calculate outward normal at a position
0122     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0123
0124   private:
0125     //// DATA ////
0126
0127     //! Location of center of toroid
0128     Real3 origin_;
0129
0130     // Radii
0131     //! Radius from origin to center of revolved ellipse (along xy plane)
0132     real_type r_;
0133     //! Horizontal radius of revolved ellipse (along xy plane)
0134     real_type a_;
0135     //! Vertical radius of revolved ellipse (along z axis)
0136     real_type b_;
0137
0138     //// HELPER FUNCTIONS ////
0139
0140     // Calculate the coefficients of the polynomial for ray intersection
0141     inline CELER_FUNCTION Real5 calc_intersection_polynomial(
0142         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0143
0144     // Shorthand to square a number
0145     static CELER_CONSTEXPR_FUNCTION real_type sq(real_type val);
0146 };
0147
0148 //---------------------------------------------------------------------------//
0149 // INLINE DEFINITIONS
0150 //---------------------------------------------------------------------------//
0151 /*!
0152  * Construct from raw data.
0153  */
0154 template<class R>
0155 CELER_FUNCTION Toroid::Toroid(Span<R, StorageSpan::extent> data)
0156     : origin_{data[0], data[1], data[2]}, r_{data[3]}, a_{data[4]}, b_{data[5]}
0157 {
0158 }
0159
0160 //---------------------------------------------------------------------------//
0161 /*!
0162  * Determine the sense of the position relative to this surface.
0163  *
0164  * For a toroid, being inside the toroid (i) counts as inside, outside
0165  * (including in the 'hole' region) (o) as outside, and on the surface exactly
0166  * as on (s).
0167  * \verbatim
0168        ___   _________   ___
0169      /     \           /     \
0170     /       \     o   /       \
0171    |         |       |         | o
0172    |         |       |    i    s
0173     \       /         \       /
0174      \_____/___________\_____/
0175    \endverbatim
0176  */
0177 CELER_FUNCTION SignedSense Toroid::calc_sense(Real3 const& pos) const
0178 {
0179     auto [x0, y0, z0] = pos - origin_;
0180
0181     real_type val = (sq(sq(x0) + sq(y0) + sq(z0 * a_ / b_) + (sq(r_) - sq(a_)))
0182                      - (4 * sq(r_)) * (sq(x0) + sq(y0)));
0183
0184     return real_to_sense(val);
0185 }
0186
0187 //---------------------------------------------------------------------------//
0188 /*!
0189  * Calculate all possible straight-line intersections between the given ray and
0190  * this surface.
0191  */
0192 CELER_FUNCTION auto Toroid::calc_intersections(
0193     Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const
0194     -> Intersections
0195 {
0196     Real5 abcde = calc_intersection_polynomial(pos, dir, on_surface);
0197     FerrariSolver solve{};  // Default tolerance
0198     Intersections roots;
0199
0200     if (on_surface == SurfaceState::on)
0201     {
0202         auto [a, b, c, d, e] = abcde;
0203         roots = solve(Real4{a, b, c, d});
0204     }
0205     else
0206     {
0207         roots = solve(abcde);
0208     }
0209
0210     return roots;
0211 }
0212
0213 //---------------------------------------------------------------------------//
0214 /*!
0215  * Calculate the coefficients of the polynomial corresponding to the given
0216  * ray's intersections with the toroid.
0217  *
0218  * Written referencing Graphics Gems II\cite{arvo-graphicsgems-1995}.
0219  */
0220 CELER_FUNCTION auto Toroid::calc_intersection_polynomial(
0221     Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const -> Real5
0222 {
0223     auto [x0, y0, z0] = pos - origin_;
0224     auto [ax, ay, az] = make_unit_vector(dir);
0225
0226     // Intermediate terms
0227     real_type p = sq(a_) / sq(b_);
0228
0229     real_type f = 1 - sq(az);
0230     real_type g = f + p * sq(az);
0231     real_type h = 2 * (x0 * ax + y0 * ay);
0232     real_type t = sq(x0) + sq(y0);
0233     real_type q = 4 * sq(r_) / sq(g);
0234     real_type m = (h + 2 * p * z0 * az) / g;
0235     real_type u = (t + p * sq(z0) + sq(r_) - sq(a_)) / g;
0236
0237     // Polynomial coefficients, i.e. cn*x^n
0238     real_type c4 = 1;
0239     real_type c3 = 2 * m;
0240     real_type c2 = sq(m) + 2 * u - q * f;
0241     real_type c1 = 2 * m * u - q * h;
0242     real_type c0 = on_surface == SurfaceState::on ? 0 : sq(u) - q * t;
0243     // Potential refinement of c0 if close to 0?
0244
0245     return Real5{c4, c3, c2, c1, c0};
0246 }
0247
0248 //---------------------------------------------------------------------------//
0249 /*!
0250  * Calculate outward facing normal at a position on or close to the surface.
0251  */
0252 CELER_FUNCTION auto Toroid::calc_normal(Real3 const& pos) const -> Real3
0253 {
0254     auto [x0, y0, z0] = pos - origin_;
0255
0256     real_type d = hypot(x0, y0);
0257     real_type f = 2 * (d - r_) / (d * sq(a_));
0258     Real3 n{x0 * f, y0 * f, (2 * z0) / (sq(b_))};
0259     return make_unit_vector(n);
0260 }
0261
0262 //---------------------------------------------------------------------------//
0263 /*!
0264  * Shorthand for power macro for readability, as squares are used frequently
0265  */
0266 CELER_CONSTEXPR_FUNCTION real_type Toroid::sq(real_type val)
0267 {
0268     return ipow<2>(val);
0269 }
0270
0271 }  // namespace celeritas