orange/surf/SimpleQuadric.hh

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