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