surf/detail/QuadricConeConverter.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/detail/QuadricConeConverter.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <optional>
0011
0012 #include "corecel/math/SoftEqual.hh"
0013
0014 #include "../ConeAligned.hh"
0015 #include "../SimpleQuadric.hh"
0016
0017 namespace celeritas
0018 {
0019 namespace detail
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Try to convert a simple quadric to a cone.
0024  *
0025  * The simple quadric must have already undergone normalization (one
0026  * second-order term less than zero, two positive).
0027  */
0028 class QuadricConeConverter
0029 {
0030   public:
0031     // Construct with tolerance
0032     inline QuadricConeConverter(real_type tol);
0033
0034     // Try converting to a cone with this orientation
0035     template<Axis T>
0036     std::optional<ConeAligned<T>>
0037     operator()(AxisTag<T>, SimpleQuadric const& sq) const;
0038
0039   private:
0040     SoftEqual<> soft_equal_;
0041 };
0042
0043 //---------------------------------------------------------------------------//
0044 // INLINE DEFINITIONS
0045 //---------------------------------------------------------------------------//
0046 /*!
0047  * Construct with tolerance.
0048  */
0049 QuadricConeConverter::QuadricConeConverter(real_type tol) : soft_equal_{tol} {}
0050
0051 //---------------------------------------------------------------------------//
0052 /*!
0053  * Try converting to a cone with this orientation.
0054  *
0055  * Cone:
0056  * \verbatim
0057     -t^2 (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = 0
0058  * \endverbatim
0059  * Expanded:
0060  * \verbatim
0061            -t^2 x^2    + y^2    + z^2
0062     + 2 t^2 x_0 x - 2y_0 y - 2z_0 z
0063     + (- t^2 x_0^2 + y_0^2 + z_0^2) = 0
0064  * \endverbatim
0065
0066  * SQ:
0067  * \verbatim
0068    a x^2 + b y^2 + (b + \epsilon) z^2 + e x + f y + g z  + h = 0
0069  * \endverbatim
0070  * Normalized (\c epsilon = 0)
0071  * \verbatim
0072    a/b x^2 + y^2 + z^2 + e/b x + f/b y + g/b z  + h/b = 0
0073  * \endverbatim
0074  * Match terms:
0075  * \verbatim
0076         -t^2  = a/b
0077     2 t^2 x_0 = e/b --> x_0 = e / (2 * a)
0078        -2 y_0 = f/b --> y_0 = f / (-2 * b)
0079        -2 z_0 = g/b --> z_0 = g / (-2 * b)
0080    - t^2 x_0^2 + y_0^2 + z_0^2 = h/b    [if not, it's a hyperboloid!]
0081  * \endverbatim
0082  */
0083 template<Axis T>
0084 std::optional<ConeAligned<T>>
0085 QuadricConeConverter::operator()(AxisTag<T>, SimpleQuadric const& sq) const
0086 {
0087     // Other coordinate system
0088     constexpr auto U = ConeAligned<T>::u_axis();
0089     constexpr auto V = ConeAligned<T>::v_axis();
0090
0091     auto second = sq.second();
0092     if (!(second[to_int(T)] < 0))
0093     {
0094         // Not the negative component
0095         return {};
0096     }
0097     if (!soft_equal_(second[to_int(U)], second[to_int(V)]))
0098     {
0099         // Not a circular cone
0100         return {};
0101     }
0102     CELER_ASSERT(second[to_int(U)] > 0 && second[to_int(V)] > 0);
0103
0104     // Normalize so U, V second-order coefficients are 1
0105     real_type const norm = (second[to_int(U)] + second[to_int(V)]) / 2;
0106     real_type const tsq = -second[to_int(T)] / norm;
0107
0108     // Calculate origin from first-order coefficients
0109     Real3 origin = make_array(sq.first());
0110     origin[to_int(T)] /= -2 * second[to_int(T)];
0111     origin[to_int(U)] /= -2 * norm;
0112     origin[to_int(V)] /= -2 * norm;
0113
0114     real_type const expected_h_b = -tsq * ipow<2>(origin[to_int(T)])
0115                                    + ipow<2>(origin[to_int(U)])
0116                                    + ipow<2>(origin[to_int(V)]);
0117     if (!soft_equal_(expected_h_b, sq.zeroth() / norm))
0118     {
0119         // Leftover constant: it's a hyperboloid!
0120         return {};
0121     }
0122
0123     // Clear potential signed zeros before returning
0124     origin += real_type{0};
0125     return ConeAligned<T>::from_tangent_sq(origin, tsq);
0126 }
0127
0128 //---------------------------------------------------------------------------//
0129 }  // namespace detail
0130 }  // namespace celeritas