orange/surf/CylCentered.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-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/CylCentered.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Macros.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 /*!
0022  * Axis-aligned cylinder centered about the origin.
0023  *
0024  * The cylinder is centered along an Axis template parameter.
0025  *
0026  * For a cylinder along the x axis:
0027  * \f[
0028     y^2 + z^2 - R^2 = 0
0029    \f]
0030  *
0031  * This is an optimization of the Cyl. The motivations are:
0032  * - Many geometries have units with concentric cylinders centered about the
0033  *   origin, so having this as a special case reduces the memory usage of those
0034  *   units (improved memory localization).
0035  * - Cylindrical mesh geometries have lots of these cylinders, so efficient
0036  *   tracking through its cells should make this optimization worthwhile.
0037  */
0038 template<Axis T>
0039 class CylCentered
0040 {
0041   public:
0042     //@{
0043     //! \name Type aliases
0044     using Intersections = Array<real_type, 2>;
0045     using StorageSpan = Span<real_type const, 1>;
0046     //@}
0047
0048   private:
0049     static constexpr Axis U{T == Axis::x ? Axis::y : Axis::x};
0050     static constexpr Axis V{T == Axis::z ? Axis::y : Axis::z};
0051
0052   public:
0053     //// CLASS ATTRIBUTES ////
0054
0055     // Surface type identifier
0056     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type();
0057
0058     //! Safety is intersection along surface normal
0059     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return true; }
0060
0061     //!@{
0062     //! Axes
0063     static CELER_CONSTEXPR_FUNCTION Axis t_axis() { return T; }
0064     static CELER_CONSTEXPR_FUNCTION Axis u_axis() { return U; }
0065     static CELER_CONSTEXPR_FUNCTION Axis v_axis() { return V; }
0066     //!@}
0067
0068   public:
0069     //// CONSTRUCTORS ////
0070
0071     // Construct with square of radius for simplification
0072     static inline CylCentered from_radius_sq(real_type rsq);
0073
0074     // Construct with radius
0075     explicit inline CELER_FUNCTION CylCentered(real_type radius);
0076
0077     // Construct from raw data
0078     template<class R>
0079     explicit inline CELER_FUNCTION CylCentered(Span<R, StorageSpan::extent>);
0080
0081     //// ACCESSORS ////
0082
0083     //! Get the square of the radius
0084     CELER_FUNCTION real_type radius_sq() const { return radius_sq_; }
0085
0086     //! Get a view to the data for type-deleted storage
0087     CELER_FUNCTION StorageSpan data() const { return {&radius_sq_, 1}; }
0088
0089     //// CALCULATION ////
0090
0091     // Determine the sense of the position relative to this surface
0092     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0093
0094     // Calculate all possible straight-line intersections with this surface
0095     inline CELER_FUNCTION Intersections calc_intersections(
0096         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0097
0098     // Calculate outward normal at a position
0099     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0100
0101   private:
0102     //! Square of cylinder radius
0103     real_type radius_sq_;
0104
0105     //! Private default constructor for manual construction
0106     CylCentered() = default;
0107 };
0108
0109 //---------------------------------------------------------------------------//
0110 // TYPE ALIASES
0111 //---------------------------------------------------------------------------//
0112
0113 using CCylX = CylCentered<Axis::x>;
0114 using CCylY = CylCentered<Axis::y>;
0115 using CCylZ = CylCentered<Axis::z>;
0116
0117 //---------------------------------------------------------------------------//
0118 // INLINE DEFINITIONS
0119 //---------------------------------------------------------------------------//
0120 /*!
0121  * Surface type identifier.
0122  */
0123 template<Axis T>
0124 CELER_CONSTEXPR_FUNCTION SurfaceType CylCentered<T>::surface_type()
0125 {
0126     return T == Axis::x   ? SurfaceType::cxc
0127            : T == Axis::y ? SurfaceType::cyc
0128            : T == Axis::z ? SurfaceType::czc
0129                           : SurfaceType::size_;
0130 }
0131
0132 //---------------------------------------------------------------------------//
0133 /*!
0134  * Construct from the square of the radius.
0135  *
0136  * This is used for surface simplification.
0137  */
0138 template<Axis T>
0139 CylCentered<T> CylCentered<T>::from_radius_sq(real_type rsq)
0140 {
0141     CELER_EXPECT(rsq > 0);
0142
0143     CylCentered<T> result;
0144     result.radius_sq_ = rsq;
0145     return result;
0146 }
0147
0148 //---------------------------------------------------------------------------//
0149 /*!
0150  * Construct with radius.
0151  */
0152 template<Axis T>
0153 CELER_FUNCTION CylCentered<T>::CylCentered(real_type radius)
0154     : radius_sq_(ipow<2>(radius))
0155 {
0156     CELER_EXPECT(radius > 0);
0157 }
0158
0159 //---------------------------------------------------------------------------//
0160 /*!
0161  * Construct from raw data.
0162  */
0163 template<Axis T>
0164 template<class R>
0165 CELER_FUNCTION CylCentered<T>::CylCentered(Span<R, StorageSpan::extent> data)
0166     : radius_sq_(data[0])
0167 {
0168 }
0169
0170 //---------------------------------------------------------------------------//
0171 /*!
0172  * Determine the sense of the position relative to this surface.
0173  */
0174 template<Axis T>
0175 CELER_FUNCTION SignedSense CylCentered<T>::calc_sense(Real3 const& pos) const
0176 {
0177     real_type const u = pos[to_int(U)];
0178     real_type const v = pos[to_int(V)];
0179
0180     return real_to_sense(ipow<2>(u) + ipow<2>(v) - radius_sq_);
0181 }
0182
0183 //---------------------------------------------------------------------------//
0184 /*!
0185  * Calculate all possible straight-line intersections with this surface.
0186  */
0187 template<Axis T>
0188 CELER_FUNCTION auto
0189 CylCentered<T>::calc_intersections(Real3 const& pos,
0190                                    Real3 const& dir,
0191                                    SurfaceState on_surface) const -> Intersections
0192 {
0193     // 1 - (\omega \dot t)^2 where t is axis of cylinder
0194     real_type const a = 1 - ipow<2>(dir[to_int(T)]);
0195
0196     if (a < ipow<2>(Tolerance<>::sqrt_quadratic()))
0197     {
0198         // No intersection if we're traveling along the cylinder axis
0199         return {no_intersection(), no_intersection()};
0200     }
0201
0202     real_type const u = pos[to_int(U)];
0203     real_type const v = pos[to_int(V)];
0204
0205     // b/2 = \omega \dot (x - x_0)
0206     detail::QuadraticSolver solve_quadric(
0207         a, dir[to_int(U)] * u + dir[to_int(V)] * v);
0208     if (on_surface == SurfaceState::on)
0209     {
0210         // Solve degenerate case (c=0)
0211         return solve_quadric();
0212     }
0213
0214     // c = (x - x_0) \dot (x - x_0) - R * R
0215     return solve_quadric(ipow<2>(u) + ipow<2>(v) - radius_sq_);
0216 }
0217
0218 //---------------------------------------------------------------------------//
0219 /*!
0220  * Calculate outward normal at a position.
0221  */
0222 template<Axis T>
0223 CELER_FUNCTION Real3 CylCentered<T>::calc_normal(Real3 const& pos) const
0224 {
0225     Real3 norm{0, 0, 0};
0226
0227     norm[to_int(U)] = pos[to_int(U)];
0228     norm[to_int(V)] = pos[to_int(V)];
0229
0230     return make_unit_vector(norm);
0231 }
0232
0233 //---------------------------------------------------------------------------//
0234 }  // namespace celeritas