orange/surf/SimpleQuadric.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/SimpleQuadric.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 ConeAligned;
0023 template<Axis T>
0024 class CylAligned;
0025 class Plane;
0026 class Sphere;
0027
0028 //---------------------------------------------------------------------------//
0029 /*!
0030  * General quadric expression but with no off-axis terms.
0031  *
0032  * Stored:
0033  * \f[
0034    ax^2 + by^2 + cz^2 + dx + ey + fz + g = 0
0035   \f]
0036  *
0037  * This can represent axis-aligned hyperboloids, ellipsoids, elliptical
0038  * cylinders, etc.
0039  */
0040 class SimpleQuadric
0041 {
0042   public:
0043     //@{
0044     //! \name Type aliases
0045     using Intersections = Array<real_type, 2>;
0046     using StorageSpan = Span<real_type const, 7>;
0047     using SpanConstReal3 = Span<real_type const, 3>;
0048     //@}
0049
0050     //// CLASS ATTRIBUTES ////
0051
0052     //! Surface type identifier
0053     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type()
0054     {
0055         return SurfaceType::sq;
0056     }
0057
0058     //! Safety is *not* the intersection along surface normal
0059     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; }
0060
0061   public:
0062     //// CONSTRUCTORS ////
0063
0064     // Construct from coefficients
0065     SimpleQuadric(Real3 const& abc, Real3 const& def, real_type g);
0066
0067     // Construct from raw data
0068     template<class R>
0069     explicit inline CELER_FUNCTION SimpleQuadric(Span<R, StorageSpan::extent>);
0070
0071     // Promote from a plane
0072     explicit SimpleQuadric(Plane const& other) noexcept;
0073
0074     // Promote from an axis-aligned cylinder
0075     template<Axis T>
0076     explicit SimpleQuadric(CylAligned<T> const& other) noexcept;
0077
0078     // Promote from a sphere
0079     explicit SimpleQuadric(Sphere const& other) noexcept;
0080
0081     // Promote from a cone
0082     template<Axis T>
0083     explicit SimpleQuadric(ConeAligned<T> const& other) noexcept;
0084
0085     //// ACCESSORS ////
0086
0087     //! Second-order terms
0088     CELER_FUNCTION SpanConstReal3 second() const { return {&a_, 3}; }
0089
0090     //! First-order terms
0091     CELER_FUNCTION SpanConstReal3 first() const { return {&d_, 3}; }
0092
0093     //! Zeroth-order term
0094     CELER_FUNCTION real_type zeroth() const { return g_; }
0095
0096     //! Get a view to the data for type-deleted storage
0097     CELER_FUNCTION StorageSpan data() const { return {&a_, 7}; }
0098
0099     //// CALCULATION ////
0100
0101     // Determine the sense of the position relative to this surface
0102     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0103
0104     // Calculate all possible straight-line intersections with this surface
0105     inline CELER_FUNCTION Intersections calc_intersections(
0106         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0107
0108     // Calculate outward normal at a position
0109     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0110
0111   private:
0112     // Second-order terms (a, b, c)
0113     real_type a_, b_, c_;
0114     // First-order terms (d, e, f)
0115     real_type d_, e_, f_;
0116     // Constant term
0117     real_type g_;
0118 };
0119
0120 //---------------------------------------------------------------------------//
0121 // INLINE DEFINITIONS
0122 //---------------------------------------------------------------------------//
0123 /*!
0124  * Construct from raw data.
0125  */
0126 template<class R>
0127 CELER_FUNCTION SimpleQuadric::SimpleQuadric(Span<R, StorageSpan::extent> data)
0128     : a_{data[0]}
0129     , b_{data[1]}
0130     , c_{data[2]}
0131     , d_{data[3]}
0132     , e_{data[4]}
0133     , f_{data[5]}
0134     , g_{data[6]}
0135 {
0136 }
0137
0138 //---------------------------------------------------------------------------//
0139 /*!
0140  * Determine the sense of the position relative to this surface.
0141  */
0142 CELER_FUNCTION SignedSense SimpleQuadric::calc_sense(Real3 const& pos) const
0143 {
0144     real_type const x = pos[to_int(Axis::x)];
0145     real_type const y = pos[to_int(Axis::y)];
0146     real_type const z = pos[to_int(Axis::z)];
0147
0148     return real_to_sense((a_ * ipow<2>(x) + b_ * ipow<2>(y) + c_ * ipow<2>(z))
0149                          + (d_ * x + e_ * y + f_ * z) + (g_));
0150 }
0151
0152 //---------------------------------------------------------------------------//
0153 /*!
0154  * Calculate all possible straight-line intersections with this surface.
0155  */
0156 CELER_FUNCTION auto
0157 SimpleQuadric::calc_intersections(Real3 const& pos,
0158                                   Real3 const& dir,
0159                                   SurfaceState on_surface) const -> Intersections
0160 {
0161     real_type const x = pos[to_int(Axis::x)];
0162     real_type const y = pos[to_int(Axis::y)];
0163     real_type const z = pos[to_int(Axis::z)];
0164     real_type const u = dir[to_int(Axis::x)];
0165     real_type const v = dir[to_int(Axis::y)];
0166     real_type const w = dir[to_int(Axis::z)];
0167
0168     // Quadratic values
0169     real_type a = (a_ * u) * u + (b_ * v) * v + (c_ * w) * w;
0170     real_type b = (2 * a_ * x + d_) * u + (2 * b_ * y + e_) * v
0171                   + (2 * c_ * z + f_) * w;
0172     real_type c = (a_ * x + d_) * x + (b_ * y + e_) * y + (c_ * z + f_) * z
0173                   + g_;
0174
0175     return detail::QuadraticSolver::solve_general(a, b / 2, c, on_surface);
0176 }
0177
0178 //---------------------------------------------------------------------------//
0179 /*!
0180  * Calculate outward normal at a position.
0181  */
0182 CELER_FUNCTION Real3 SimpleQuadric::calc_normal(Real3 const& pos) const
0183 {
0184     real_type const x = pos[to_int(Axis::x)];
0185     real_type const y = pos[to_int(Axis::y)];
0186     real_type const z = pos[to_int(Axis::z)];
0187
0188     Real3 norm{2 * a_ * x + d_, 2 * b_ * y + e_, 2 * c_ * z + f_};
0189     return make_unit_vector(norm);
0190 }
0191
0192 //---------------------------------------------------------------------------//
0193 }  // namespace celeritas