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.  The quadric is ill-defined if all non-constants are zero.
0041  *
0042  * Even though this equation can be scaled arbitrarily, the quadratic solver
0043  * is sensitive to the scale of the coefficients to avoid numerical precision
0044  * loss. If present, the first-order coefficients should be on the length scale
0045  * of the problem (i.e., translations avoid catastrophic precision loss).
0046  * In most non-pathological cases, the scale of the second-order components
0047  * should be unity, as is the case when promoting spheres, cones, and
0048  * cylinders.
0049  */
0050 class SimpleQuadric
0051 {
0052   public:
0053     //@{
0054     //! \name Type aliases
0055     using Intersections = Array<real_type, 2>;
0056     using StorageSpan = Span<real_type const, 7>;
0057     using SpanConstReal3 = Span<real_type const, 3>;
0058     //@}
0059
0060     //// CLASS ATTRIBUTES ////
0061
0062     //! Surface type identifier
0063     static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type()
0064     {
0065         return SurfaceType::sq;
0066     }
0067
0068     //! Safety is *not* the intersection along surface normal
0069     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; }
0070
0071   public:
0072     //// CONSTRUCTORS ////
0073
0074     // Construct from coefficients
0075     SimpleQuadric(Real3 const& abc, Real3 const& def, real_type g);
0076
0077     // Construct from raw data
0078     template<class R>
0079     explicit inline CELER_FUNCTION SimpleQuadric(Span<R, StorageSpan::extent>);
0080
0081     // Promote from a plane
0082     explicit SimpleQuadric(Plane const& other) noexcept(!CELERITAS_DEBUG);
0083
0084     // Promote from an axis-aligned cylinder
0085     template<Axis T>
0086     explicit SimpleQuadric(CylAligned<T> const& other) noexcept;
0087
0088     // Promote from a sphere
0089     explicit SimpleQuadric(Sphere const& other) noexcept(!CELERITAS_DEBUG);
0090
0091     // Promote from a cone
0092     template<Axis T>
0093     explicit SimpleQuadric(ConeAligned<T> const& other) noexcept;
0094
0095     //// ACCESSORS ////
0096
0097     //! Second-order terms
0098     CELER_FUNCTION SpanConstReal3 second() const { return {&a_, 3}; }
0099
0100     //! First-order terms
0101     CELER_FUNCTION SpanConstReal3 first() const { return {&d_, 3}; }
0102
0103     //! Zeroth-order term
0104     CELER_FUNCTION real_type zeroth() const { return g_; }
0105
0106     //! Get a view to the data for type-deleted storage
0107     CELER_FUNCTION StorageSpan data() const { return {&a_, 7}; }
0108
0109     //// CALCULATION ////
0110
0111     // Determine the sense of the position relative to this surface
0112     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0113
0114     // Calculate all possible straight-line intersections with this surface
0115     inline CELER_FUNCTION Intersections calc_intersections(
0116         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0117
0118     // Calculate outward normal at a position
0119     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0120
0121   private:
0122     // Second-order terms (a, b, c)
0123     real_type a_, b_, c_;
0124     // First-order terms (d, e, f)
0125     real_type d_, e_, f_;
0126     // Constant term
0127     real_type g_;
0128 };
0129
0130 //---------------------------------------------------------------------------//
0131 // INLINE DEFINITIONS
0132 //---------------------------------------------------------------------------//
0133 /*!
0134  * Construct from raw data.
0135  */
0136 template<class R>
0137 CELER_FUNCTION SimpleQuadric::SimpleQuadric(Span<R, StorageSpan::extent> data)
0138     : a_{data[0]}
0139     , b_{data[1]}
0140     , c_{data[2]}
0141     , d_{data[3]}
0142     , e_{data[4]}
0143     , f_{data[5]}
0144     , g_{data[6]}
0145 {
0146 }
0147
0148 //---------------------------------------------------------------------------//
0149 /*!
0150  * Determine the sense of the position relative to this surface.
0151  */
0152 CELER_FUNCTION SignedSense SimpleQuadric::calc_sense(Real3 const& pos) const
0153 {
0154     real_type const x = pos[to_int(Axis::x)];
0155     real_type const y = pos[to_int(Axis::y)];
0156     real_type const z = pos[to_int(Axis::z)];
0157
0158     return real_to_sense((a_ * ipow<2>(x) + b_ * ipow<2>(y) + c_ * ipow<2>(z))
0159                          + (d_ * x + e_ * y + f_ * z) + (g_));
0160 }
0161
0162 //---------------------------------------------------------------------------//
0163 /*!
0164  * Calculate all possible straight-line intersections with this surface.
0165  */
0166 CELER_FUNCTION auto
0167 SimpleQuadric::calc_intersections(Real3 const& pos,
0168                                   Real3 const& dir,
0169                                   SurfaceState on_surface) const
0170     -> Intersections
0171 {
0172     real_type const x = pos[to_int(Axis::x)];
0173     real_type const y = pos[to_int(Axis::y)];
0174     real_type const z = pos[to_int(Axis::z)];
0175     real_type const u = dir[to_int(Axis::x)];
0176     real_type const v = dir[to_int(Axis::y)];
0177     real_type const w = dir[to_int(Axis::z)];
0178
0179     // Quadratic values
0180     real_type a = (a_ * u) * u + (b_ * v) * v + (c_ * w) * w;
0181     real_type b = (2 * a_ * x + d_) * u + (2 * b_ * y + e_) * v
0182                   + (2 * c_ * z + f_) * w;
0183     real_type c = (a_ * x + d_) * x + (b_ * y + e_) * y + (c_ * z + f_) * z
0184                   + g_;
0185
0186     return detail::QuadraticSolver::solve_general(a, b / 2, c, on_surface);
0187 }
0188
0189 //---------------------------------------------------------------------------//
0190 /*!
0191  * Calculate outward normal at a position.
0192  */
0193 CELER_FUNCTION Real3 SimpleQuadric::calc_normal(Real3 const& pos) const
0194 {
0195     real_type const x = pos[to_int(Axis::x)];
0196     real_type const y = pos[to_int(Axis::y)];
0197     real_type const z = pos[to_int(Axis::z)];
0198
0199     Real3 norm{2 * a_ * x + d_, 2 * b_ * y + e_, 2 * c_ * z + f_};
0200     return make_unit_vector(norm);
0201 }
0202
0203 //---------------------------------------------------------------------------//
0204 }  // namespace celeritas