orange/surf/GeneralQuadric.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/GeneralQuadric.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 class SimpleQuadric;
0024
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * General quadric surface.
0028  *
0029  * General quadrics that cannot be simplified to other ORANGE surfaces include
0030  * hyperboloids and paraboloids; and non-axis-aligned cylinders, ellipsoids,
0031  * and cones.
0032  *
0033  * \f[
0034     ax^2 + by^2 + cz^2 + dxy + eyz + fzx + gx + hy + iz + j = 0
0035    \f]
0036  *
0037  * Note that some formulations of a general quadric include a factor of 2 for
0038  * the g/h/i terms.
0039  */
0040 class GeneralQuadric
0041 {
0042   public:
0043     //@{
0044     //! Type aliases
0045     using Intersections = Array<real_type, 2>;
0046     using StorageSpan = Span<real_type const, 10>;
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::gq;
0056     }
0057
0058     //! Safety is *not* the nearest intersection along the surface "normal"
0059     static CELER_CONSTEXPR_FUNCTION bool simple_safety() { return false; }
0060
0061   public:
0062     //// CONSTRUCTORS ////
0063
0064     // Construct from coefficients
0065     explicit GeneralQuadric(Real3 const& abc,
0066                             Real3 const& def,
0067                             Real3 const& ghi,
0068                             real_type j);
0069
0070     // Construct from raw data
0071     template<class R>
0072     explicit inline CELER_FUNCTION GeneralQuadric(Span<R, StorageSpan::extent>);
0073
0074     // Promote from a simple quadric
0075     explicit GeneralQuadric(SimpleQuadric const& other) noexcept(
0076         !CELERITAS_DEBUG);
0077
0078     //// ACCESSORS ////
0079
0080     //! Second-order terms
0081     CELER_FUNCTION SpanConstReal3 second() const { return {&a_, 3}; }
0082
0083     //! Cross terms (xy, yz, zx)
0084     CELER_FUNCTION SpanConstReal3 cross() const { return {&d_, 3}; }
0085
0086     //! First-order terms
0087     CELER_FUNCTION SpanConstReal3 first() const { return {&g_, 3}; }
0088
0089     //! Zeroth-order term
0090     CELER_FUNCTION real_type zeroth() const { return j_; }
0091
0092     //! Get a view to the data for type-deleted storage
0093     CELER_FUNCTION StorageSpan data() const { return {&a_, 10}; }
0094
0095     //// CALCULATION ////
0096
0097     // Determine the sense of the position relative to this surface
0098     inline CELER_FUNCTION SignedSense calc_sense(Real3 const& pos) const;
0099
0100     // Calculate all possible straight-line intersections with this surface
0101     inline CELER_FUNCTION Intersections calc_intersections(
0102         Real3 const& pos, Real3 const& dir, SurfaceState on_surface) const;
0103
0104     // Calculate outward normal at a position
0105     inline CELER_FUNCTION Real3 calc_normal(Real3 const& pos) const;
0106
0107   private:
0108     // Second-order terms (a, b, c)
0109     real_type a_, b_, c_;
0110     // Second-order cross terms (d, e, f)
0111     real_type d_, e_, f_;
0112     // First-order terms (g, h, i)
0113     real_type g_, h_, i_;
0114     // Constant term
0115     real_type j_;
0116 };
0117
0118 //---------------------------------------------------------------------------//
0119 // INLINE DEFINITIONS
0120 //---------------------------------------------------------------------------//
0121 /*!
0122  * Construct from raw data.
0123  */
0124 template<class R>
0125 CELER_FUNCTION GeneralQuadric::GeneralQuadric(Span<R, StorageSpan::extent> data)
0126     : a_(data[0])
0127     , b_(data[1])
0128     , c_(data[2])
0129     , d_(data[3])
0130     , e_(data[4])
0131     , f_(data[5])
0132     , g_(data[6])
0133     , h_(data[7])
0134     , i_(data[8])
0135     , j_(data[9])
0136 {
0137 }
0138
0139 //---------------------------------------------------------------------------//
0140 /*!
0141  * Determine the sense of the position relative to this surface.
0142  */
0143 CELER_FUNCTION SignedSense GeneralQuadric::calc_sense(Real3 const& pos) const
0144 {
0145     real_type const x = pos[0];
0146     real_type const y = pos[1];
0147     real_type const z = pos[2];
0148
0149     real_type result = (a_ * x + d_ * y + f_ * z + g_) * x
0150                        + (b_ * y + e_ * z + h_) * y + (c_ * z + i_) * z + j_;
0151
0152     return real_to_sense(result);
0153 }
0154
0155 //---------------------------------------------------------------------------//
0156 /*!
0157  * Calculate all possible straight-line intersections with this surface.
0158  */
0159 CELER_FUNCTION auto
0160 GeneralQuadric::calc_intersections(Real3 const& pos,
0161                                    Real3 const& dir,
0162                                    SurfaceState on_surface) const -> Intersections
0163 {
0164     real_type const x = pos[0];
0165     real_type const y = pos[1];
0166     real_type const z = pos[2];
0167     real_type const u = dir[0];
0168     real_type const v = dir[1];
0169     real_type const w = dir[2];
0170
0171     // Quadratic values
0172     real_type a = (a_ * u + d_ * v) * u + (b_ * v + e_ * w) * v
0173                   + (c_ * w + f_ * u) * w;
0174     real_type b = (2 * a_ * x + d_ * y + f_ * z + g_) * u
0175                   + (2 * b_ * y + d_ * x + e_ * z + h_) * v
0176                   + (2 * c_ * z + e_ * y + f_ * x + i_) * w;
0177     real_type c = ((a_ * x + d_ * y + g_) * x + (b_ * y + e_ * z + h_) * y
0178                    + (c_ * z + f_ * x + i_) * z + j_);
0179
0180     return detail::QuadraticSolver::solve_general(a, b / 2, c, on_surface);
0181 }
0182
0183 //---------------------------------------------------------------------------//
0184 /*!
0185  * Calculate outward normal at a position.
0186  */
0187 CELER_FUNCTION Real3 GeneralQuadric::calc_normal(Real3 const& pos) const
0188 {
0189     real_type const x = pos[0];
0190     real_type const y = pos[1];
0191     real_type const z = pos[2];
0192
0193     Real3 norm;
0194     norm[0] = 2 * a_ * x + d_ * y + f_ * z + g_;
0195     norm[1] = 2 * b_ * y + d_ * x + e_ * z + h_;
0196     norm[2] = 2 * c_ * z + e_ * y + f_ * x + i_;
0197
0198     return make_unit_vector(norm);
0199 }
0200
0201 //---------------------------------------------------------------------------//
0202 }  // namespace celeritas