orange/surf/CylAligned.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2023-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file orange/surf/CylAligned.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Types.hh"
0011 #include "corecel/cont/Array.hh"
0012 #include "corecel/cont/Span.hh"
0013 #include "corecel/math/ArrayUtils.hh"
0014 #include "orange/OrangeTypes.hh"
0015
0016 #include "detail/QuadraticSolver.hh"
0017
0018 namespace celeritas
0019 {
0020 //---------------------------------------------------------------------------//
0021 template<Axis T>
0022 class CylCentered;
0023
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Axis-aligned cylinder.
0027  *
0028  * The cylinder is centered about the template parameter Axis.
0029  *
0030  * For a cylinder along the x axis:
0031  * \f[
0032     (y - y_0)^2 + (z - z_0)^2 - R^2 = 0
0033    \f]
0034  *
0035  * The axis-aligned cylinder is used to calculate a tolerance for the
0036  * second-order coefficient on the quadratic equation, where a small value is
0037  * approximately zero.
0038  *
0039  * Basing the intersection tolerance on a particle traveling nearly parallel to
0040  * the cylinder, the angle \f$\theta \f$
0041  * between the particle's direction and the center axis of the cylinder needed
0042  * to give an intersection distance of \f$ 1/\delta \f$ when a distance of 1
0043  * from the cylinder (or equivalently, an intersection of 1 when a distance
0044  * \f$ \delta \f$ away) is
0045  * \f[
0046    \sin \theta = \frac{1}{1/\delta} \;.
0047    \f]
0048  * Letting \f$ \mu = \cos \theta = \Omega \cdot t \f$ we have \f[
0049    \mu^2 + \delta^2 = 1 \;.
0050    \f]
0051  * The quadratic intersection for an axis-aligned cylinder where the particle
0052  * is at \f$ R + \delta \f$ is
0053  * \f[
0054    0 = ax^2 + bx + c
0055      = (1 - \mu^2)x^2 + 2 \sqrt{1 - \mu^2} (R + 1) x + (R + 1)^2 - R^2
0056   \f]
0057  * so the quadric coefficient when a particle is effectively parallel to a
0058  * cylinder is:
0059   \f[
0060   a = (1 - \mu^2) = \delta^2 \;.
0061   \f]
0062  *
0063  * Because \em a is calculated by subtraction, this puts a hard limit on the
0064  * accuracy of the intersection distance: the default value of a=1e-10
0065  * corresponds to a tolerance of 1e-5, not the default tolerance 1e-8.
0066  */
0067 template<Axis T>
0068 class CylAligned
0069 {
0070   public:
0071     //@{
0072     //! \name Type aliases
0073     using Intersections = Array<real_type, 2>;
0074     using StorageSpan = Span<real_type const, 3>;
0075     //@}
0076
0077   private:
0078     static constexpr Axis U{T == Axis::x ? Axis::y : Axis::x};
0079     static constexpr Axis V{T == Axis::z ? Axis::y : Axis::z};
0080
0081   public:
0082     //// CLASS ATTRIBUTES ////
0083
0084     // Surface type identifier
0085     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type();
0086
0087     //! Safety is intersection along surface normal
0088     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; }
0089
0090     //!@{
0091     //! Axes
0092     static CELER_CONSTEXPR_FUNCTION Axis t_axis() { return T; }
0093     static CELER_CONSTEXPR_FUNCTION Axis u_axis() { return U; }
0094     static CELER_CONSTEXPR_FUNCTION Axis v_axis() { return V; }
0095     //!@}
0096
0097   public:
0098     //// CONSTRUCTORS ////
0099
0100     // Construct with square of radius for simplification
0101     static CylAligned from_radius_sq(Real3 const& origin, real_type rsq);
0102
0103     // Construct with radius
0104     inline CELER_FUNCTION CylAligned(Real3 const& origin, real_type radius);
0105
0106     // Construct from raw data
0107     template<class R>
0108     explicit inline CELER_FUNCTION CylAligned(Span<R, StorageSpan::extent>);
0109
0110     // Promote from a centered axis-aligned cylinder
0111     explicit CylAligned(CylCentered<T> const& other) noexcept;
0112
0113     //// ACCESSORS ////
0114
0115     //! Get the origin vector along the 'u' axis
0116     CELER_FUNCTION real_type origin_u() const { return origin_u_; }
0117
0118     //! Get the origin vector along the 'v' axis
0119     CELER_FUNCTION real_type origin_v() const { return origin_v_; }
0120
0121     //! Get the square of the radius
0122     CELER_FUNCTION real_type radius_sq() const { return radius_sq_; }
0123
0124     //! Get a view to the data for type-deleted storage
0125     CELER_FUNCTION StorageSpan data() const { return {&origin_u_, 3}; }
0126
0127     // Helper function to get the origin as a 3-vector
0128     Real3 calc_origin() const;
0129
0130     //// CALCULATION ////
0131
0132     // Determine the sense of the position relative to this surface
0133     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0134
0135     // Calculate all possible straight-line intersections with this surface
0136     inline CELER_FUNCTION Intersections calc_intersections(
0137         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0138
0139     // Calculate outward normal at a position
0140     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0141
0142   private:
0143     // Off-axis location
0144     real_type origin_u_;
0145     real_type origin_v_;
0146
0147     // Square of the radius
0148     real_type radius_sq_;
0149
0150     //! Private default constructor for manual construction
0151     CylAligned() = default;
0152 };
0153
0154 //---------------------------------------------------------------------------//
0155 // TYPE ALIASES
0156 //---------------------------------------------------------------------------//
0157
0158 using CylX = CylAligned<Axis::x>;
0159 using CylY = CylAligned<Axis::y>;
0160 using CylZ = CylAligned<Axis::z>;
0161
0162 //---------------------------------------------------------------------------//
0163 // INLINE DEFINITIONS
0164 //---------------------------------------------------------------------------//
0165 /*!
0166  * Surface type identifier.
0167  */
0168 template<Axis T>
0169 CELER_CONSTEXPR_FUNCTION SurfaceType CylAligned<T>::surface_type()
0170 {
0171     return T == Axis::x   ? SurfaceType::cx
0172            : T == Axis::y ? SurfaceType::cy
0173            : T == Axis::z ? SurfaceType::cz
0174                           : SurfaceType::size_;
0175 }
0176
0177 //---------------------------------------------------------------------------//
0178 /*!
0179  * Construct from origin and radius.
0180  */
0181 template<Axis T>
0182 CELER_FUNCTION CylAligned<T>::CylAligned(Real3 const& origin, real_type radius)
0183     : origin_u_{origin[to_int(U)]}
0184     , origin_v_{origin[to_int(V)]}
0185     , radius_sq_{ipow<2>(radius)}
0186 {
0187     CELER_EXPECT(radius > 0);
0188 }
0189
0190 //---------------------------------------------------------------------------//
0191 /*!
0192  * Construct from raw data.
0193  */
0194 template<Axis T>
0195 template<class R>
0196 CELER_FUNCTION CylAligned<T>::CylAligned(Span<R, StorageSpan::extent> data)
0197     : origin_u_{data[0]}, origin_v_{data[1]}, radius_sq_{data[2]}
0198 {
0199 }
0200
0201 //---------------------------------------------------------------------------//
0202 /*!
0203  * Determine the sense of the position relative to this surface.
0204  */
0205 template<Axis T>
0206 CELER_FUNCTION SignedSense CylAligned<T>::calc_sense(Real3 const& pos) const
0207 {
0208     real_type const u = pos[to_int(U)] - origin_u_;
0209     real_type const v = pos[to_int(V)] - origin_v_;
0210
0211     return real_to_sense(ipow<2>(u) + ipow<2>(v) - radius_sq_);
0212 }
0213
0214 //---------------------------------------------------------------------------//
0215 /*!
0216  * Calculate all possible straight-line intersections with this surface.
0217  */
0218 template<Axis T>
0219 CELER_FUNCTION auto
0220 CylAligned<T>::calc_intersections(Real3 const& pos,
0221                                   Real3 const& dir,
0222                                   SurfaceState on_surface) const -> Intersections
0223 {
0224     // 1 - \omega \dot e
0225     real_type const a = 1 - ipow<2>(dir[to_int(T)]);
0226
0227     if (a < ipow<2>(Tolerance<>::sqrt_quadratic()))
0228     {
0229         // No intersection if we're traveling along the cylinder axis
0230         return {no_intersection(), no_intersection()};
0231     }
0232
0233     real_type const u = pos[to_int(U)] - origin_u_;
0234     real_type const v = pos[to_int(V)] - origin_v_;
0235
0236     // b/2 = \omega \dot (x - x_0)
0237     detail::QuadraticSolver solve_quadric(
0238         a, dir[to_int(U)] * u + dir[to_int(V)] * v);
0239     if (on_surface == SurfaceState::on)
0240     {
0241         // Solve degenerate case (c=0)
0242         return solve_quadric();
0243     }
0244
0245     // c = (x - x_0) \dot (x - x_0) - R * R
0246     return solve_quadric(ipow<2>(u) + ipow<2>(v) - radius_sq_);
0247 }
0248
0249 //---------------------------------------------------------------------------//
0250 /*!
0251  * Calculate outward normal at a position.
0252  */
0253 template<Axis T>
0254 CELER_FUNCTION Real3 CylAligned<T>::calc_normal(Real3 const& pos) const
0255 {
0256     Real3 norm{0, 0, 0};
0257
0258     norm[to_int(U)] = pos[to_int(U)] - origin_u_;
0259     norm[to_int(V)] = pos[to_int(V)] - origin_v_;
0260
0261     return make_unit_vector(norm);
0262 }
0263
0264 //---------------------------------------------------------------------------//
0265 }  // namespace celeritas