orange/surf/ConeAligned.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/ConeAligned.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 /*!
0022  * Axis-aligned cone (infinite and double-sheeted).
0023  *
0024  * For a cone parallel to the x axis:
0025  * \f[
0026     (y - y_0)^2 + (z - z_0)^2 - t^2 (x - x_0)^2 = 0
0027    \f]
0028
0029  * where \em t is the tangent of the opening angle (\f$r/h\f$ for a finite cone
0030  * with radius \em r and height \em h).
0031  */
0032 template<Axis T>
0033 class ConeAligned
0034 {
0035   public:
0036     //@{
0037     //! \name Type aliases
0038     using Intersections = Array<real_type, 2>;
0039     using StorageSpan = Span<real_type const, 4>;
0040     //@}
0041
0042   private:
0043     static constexpr Axis U{T == Axis::x ? Axis::y : Axis::x};
0044     static constexpr Axis V{T == Axis::z ? Axis::y : Axis::z};
0045
0046   public:
0047     //// CLASS ATTRIBUTES ////
0048
0049     // Surface type identifier
0050     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type();
0051
0052     //! Safety is intersection along surface normal
0053     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; }
0054
0055     //!@{
0056     //! Axes
0057     static CELER_CONSTEXPR_FUNCTION Axis t_axis() { return T; }
0058     static CELER_CONSTEXPR_FUNCTION Axis u_axis() { return U; }
0059     static CELER_CONSTEXPR_FUNCTION Axis v_axis() { return V; }
0060     //!@}
0061
0062   public:
0063     //// CONSTRUCTORS ////
0064
0065     // Construct with square of tangent for simplification
0066     static ConeAligned from_tangent_sq(Real3 const& origin, real_type tsq);
0067
0068     // Construct from origin and tangent of the angle of its opening
0069     inline CELER_FUNCTION ConeAligned(Real3 const& origin, real_type tangent);
0070
0071     // Construct from raw data
0072     template<class R>
0073     explicit inline CELER_FUNCTION ConeAligned(Span<R, StorageSpan::extent>);
0074
0075     //// ACCESSORS ////
0076
0077     //! Get the origin position along the normal axis
0078     CELER_FUNCTION Real3 const& origin() const { return origin_; }
0079
0080     //! Get the square of the tangent of the opening angle
0081     CELER_FUNCTION real_type tangent_sq() const { return tsq_; }
0082
0083     //! Get a view to the data for type-deleted storage
0084     CELER_FUNCTION StorageSpan data() const { return {&origin_[0], 4}; }
0085
0086     //// CALCULATION ////
0087
0088     // Determine the sense of the position relative to this surface
0089     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0090
0091     // Calculate all possible straight-line intersections with this surface
0092     inline CELER_FUNCTION Intersections calc_intersections(
0093         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0094
0095     // Calculate outward normal at a position
0096     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0097
0098   private:
0099     // Location of the vanishing point
0100     Real3 origin_;
0101
0102     // Quadric value
0103     real_type tsq_;
0104
0105     //! Private default constructor for manual construction
0106     ConeAligned() = default;
0107 };
0108
0109 //---------------------------------------------------------------------------//
0110 // TYPE ALIASES
0111 //---------------------------------------------------------------------------//
0112
0113 using ConeX = ConeAligned<Axis::x>;
0114 using ConeY = ConeAligned<Axis::y>;
0115 using ConeZ = ConeAligned<Axis::z>;
0116
0117 //---------------------------------------------------------------------------//
0118 // INLINE DEFINITIONS
0119 //---------------------------------------------------------------------------//
0120 /*!
0121  * Surface type identifier.
0122  */
0123 template<Axis T>
0124 CELER_CONSTEXPR_FUNCTION SurfaceType ConeAligned<T>::surface_type()
0125 {
0126     return T == Axis::x   ? SurfaceType::kx
0127            : T == Axis::y ? SurfaceType::ky
0128            : T == Axis::z ? SurfaceType::kz
0129                           : SurfaceType::size_;
0130 }
0131
0132 //---------------------------------------------------------------------------//
0133 /*!
0134  * Construct from origin and tangent of the angle of its opening.
0135  *
0136  * Given a finite cone, the tangent is the ratio of its base radius to its
0137  * height.
0138  *
0139  * In the cone, below, the tangent of the inner angle
0140  * \f$ \tan(\theta) = r/h \f$ is the second argument.
0141  * \verbatim
0142      r
0143    +-------*
0144    |   _--^
0145  h |_--
0146    O
0147    \endverbatim
0148  */
0149 template<Axis T>
0150 CELER_FUNCTION
0151 ConeAligned<T>::ConeAligned(Real3 const& origin, real_type tangent)
0152     : origin_{origin}, tsq_{ipow<2>(tangent)}
0153 {
0154     CELER_EXPECT(tangent > 0);
0155 }
0156
0157 //---------------------------------------------------------------------------//
0158 /*!
0159  * Construct from raw data.
0160  */
0161 template<Axis T>
0162 template<class R>
0163 CELER_FUNCTION ConeAligned<T>::ConeAligned(Span<R, StorageSpan::extent> data)
0164     : origin_{data[0], data[1], data[2]}, tsq_{data[3]}
0165 {
0166 }
0167
0168 //---------------------------------------------------------------------------//
0169 /*!
0170  * Determine the sense of the position relative to this surface.
0171  */
0172 template<Axis T>
0173 CELER_FUNCTION SignedSense ConeAligned<T>::calc_sense(Real3 const& pos) const
0174 {
0175     real_type const x = pos[to_int(T)] - origin_[to_int(T)];
0176     real_type const y = pos[to_int(U)] - origin_[to_int(U)];
0177     real_type const z = pos[to_int(V)] - origin_[to_int(V)];
0178
0179     return real_to_sense((-tsq_ * ipow<2>(x)) + ipow<2>(y) + ipow<2>(z));
0180 }
0181
0182 //---------------------------------------------------------------------------//
0183 /*!
0184  * Calculate all possible straight-line intersections with this surface.
0185  *
0186  * \f[
0187     (y - yc)^2 + (z - zc)^2 - t^2 * (x - xc)^2 = 0
0188    \f]
0189  */
0190 template<Axis T>
0191 CELER_FUNCTION auto
0192 ConeAligned<T>::calc_intersections(Real3 const& pos,
0193                                    Real3 const& dir,
0194                                    SurfaceState on_surface) const -> Intersections
0195 {
0196     // Expand translated positions into 'xyz' coordinate system
0197     real_type const x = pos[to_int(T)] - origin_[to_int(T)];
0198     real_type const y = pos[to_int(U)] - origin_[to_int(U)];
0199     real_type const z = pos[to_int(V)] - origin_[to_int(V)];
0200
0201     real_type const u = dir[to_int(T)];
0202     real_type const v = dir[to_int(U)];
0203     real_type const w = dir[to_int(V)];
0204
0205     // Scaled direction
0206     real_type a = (-tsq_ * ipow<2>(u)) + ipow<2>(v) + ipow<2>(w);
0207     real_type half_b = (-tsq_ * x * u) + (y * v) + (z * w);
0208     real_type c = (-tsq_ * ipow<2>(x)) + ipow<2>(y) + ipow<2>(z);
0209
0210     return detail::QuadraticSolver::solve_general(a, half_b, c, on_surface);
0211 }
0212
0213 //---------------------------------------------------------------------------//
0214 /*!
0215  * Calculate outward normal at a position.
0216  */
0217 template<Axis T>
0218 CELER_FUNCTION Real3 ConeAligned<T>::calc_normal(Real3 const& pos) const
0219 {
0220     Real3 norm;
0221     for (auto i = to_int(Axis::x); i < to_int(Axis::size_); ++i)
0222     {
0223         norm[i] = pos[i] - origin_[i];
0224     }
0225     norm[to_int(T)] *= -tsq_;
0226
0227     return make_unit_vector(norm);
0228 }
0229
0230 //---------------------------------------------------------------------------//
0231 }  // namespace celeritas