orange/surf/Sphere.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/Sphere.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/Algorithms.hh"
0014 #include "corecel/math/ArrayUtils.hh"
0015
0016 #include "detail/QuadraticSolver.hh"
0017
0018 namespace celeritas
0019 {
0020 //---------------------------------------------------------------------------//
0021 class SphereCentered;
0022
0023 //---------------------------------------------------------------------------//
0024 /*!
0025  * Sphere centered at an arbitrary point.
0026  *
0027  * \f[
0028    (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - r^2 = 0
0029  * \f]
0030  */
0031 class Sphere
0032 {
0033   public:
0034     //@{
0035     //! Type aliases
0036     using Intersections = Array<real_type, 2>;
0037     using StorageSpan = Span<real_type const, 4>;
0038     //@}
0039
0040     //// CLASS ATTRIBUTES ////
0041
0042     //! Surface type identifier
0043     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type()
0044     {
0045         return SurfaceType::s;
0046     }
0047
0048     //! Safety is intersection along surface normal
0049     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return true; }
0050
0051   public:
0052     //// CONSTRUCTORS ////
0053
0054     // Construct with square of radius for simplification
0055     static Sphere from_radius_sq(Real3 const& origin, real_type rsq);
0056
0057     // Construct with origin and radius
0058     inline CELER_FUNCTION Sphere(Real3 const& origin, real_type radius);
0059
0060     // Construct from raw data
0061     template<class R>
0062     explicit inline CELER_FUNCTION Sphere(Span<R, StorageSpan::extent>);
0063
0064     // Promote from a centered sphere
0065     explicit Sphere(SphereCentered const& other) noexcept;
0066
0067     //// ACCESSORS ////
0068
0069     //! Center position
0070     CELER_FUNCTION Real3 const& origin() const { return origin_; }
0071
0072     //! Square of the radius
0073     CELER_FUNCTION real_type radius_sq() const { return radius_sq_; }
0074
0075     //! Get a view to the data for type-deleted storage
0076     CELER_FUNCTION StorageSpan data() const { return {origin_.data(), 4}; }
0077
0078     //// CALCULATION ////
0079
0080     // Determine the sense of the position relative to this surface
0081     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0082
0083     // Calculate all possible straight-line intersections with this surface
0084     inline CELER_FUNCTION Intersections calc_intersections(
0085         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0086
0087     // Calculate outward normal at a position
0088     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0089
0090   private:
0091     // Spatial position
0092     Real3 origin_;
0093     // Square of the radius
0094     real_type radius_sq_;
0095
0096     //! Private default constructor for manual construction
0097     Sphere() = default;
0098 };
0099
0100 //---------------------------------------------------------------------------//
0101 // INLINE DEFINITIONS
0102 //---------------------------------------------------------------------------//
0103 /*!
0104  * Construct with sphere origin and radius.
0105  */
0106 CELER_FUNCTION Sphere::Sphere(Real3 const& origin, real_type radius)
0107     : origin_(origin), radius_sq_(ipow<2>(radius))
0108 {
0109     CELER_EXPECT(radius > 0);
0110 }
0111
0112 //---------------------------------------------------------------------------//
0113 /*!
0114  * Construct from raw data.
0115  */
0116 template<class R>
0117 CELER_FUNCTION Sphere::Sphere(Span<R, StorageSpan::extent> data)
0118     : origin_{data[0], data[1], data[2]}, radius_sq_{data[3]}
0119 {
0120 }
0121
0122 //---------------------------------------------------------------------------//
0123 /*!
0124  * Determine the sense of the position relative to this surface.
0125  */
0126 CELER_FUNCTION SignedSense Sphere::calc_sense(Real3 const& pos) const
0127 {
0128     Real3 tpos{pos[0] - origin_[0], pos[1] - origin_[1], pos[2] - origin_[2]};
0129
0130     return real_to_sense(dot_product(tpos, tpos) - radius_sq_);
0131 }
0132
0133 //---------------------------------------------------------------------------//
0134 /*!
0135  * Calculate all possible straight-line intersections with this surface.
0136  */
0137 CELER_FUNCTION auto
0138 Sphere::calc_intersections(Real3 const& pos,
0139                            Real3 const& dir,
0140                            SurfaceState on_surface) const -> Intersections
0141 {
0142     Real3 tpos{pos[0] - origin_[0], pos[1] - origin_[1], pos[2] - origin_[2]};
0143
0144     detail::QuadraticSolver solve_quadric(real_type(1), dot_product(tpos, dir));
0145     if (on_surface == SurfaceState::off)
0146     {
0147         return solve_quadric(dot_product(tpos, tpos) - radius_sq_);
0148     }
0149     else
0150     {
0151         return solve_quadric();
0152     }
0153 }
0154
0155 //---------------------------------------------------------------------------//
0156 /*!
0157  * Calculate outward normal at a position.
0158  */
0159 CELER_FUNCTION Real3 Sphere::calc_normal(Real3 const& pos) const
0160 {
0161     return make_unit_vector(
0162         Real3{pos[0] - origin_[0], pos[1] - origin_[1], pos[2] - origin_[2]});
0163 }
0164
0165 //---------------------------------------------------------------------------//
0166 }  // namespace celeritas