surf/detail/QuadricCylConverter.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/QuadricCylConverter.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <optional>
0011
0012 #include "corecel/math/SoftEqual.hh"
0013
0014 #include "../CylAligned.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 cylinder.
0024  *
0025  * The simple quadric must have already undergone normalization (one
0026  * second-order term approximately zero, the others positive).
0027  */
0028 class QuadricCylConverter
0029 {
0030   public:
0031     // Construct with tolerance
0032     inline QuadricCylConverter(real_type tol);
0033
0034     // Try converting to a cylinder with this orientation
0035     template<Axis T>
0036     std::optional<CylAligned<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 QuadricCylConverter::QuadricCylConverter(real_type tol) : soft_equal_{tol} {}
0050
0051 //---------------------------------------------------------------------------//
0052 /*!
0053  * Try converting to a cylinder with this orientation.
0054  *
0055  * Cone:
0056  * \verbatim
0057     (x - x_0)^2 + (y - y_0)^2 - r^2 = 0
0058  * \endverbatim
0059  * Expanded:
0060  * \verbatim
0061            x^2 +      y^2
0062     - 2x_0 x   - 2y_0 y
0063     + (x_0^2 + y_0^2 - r^2) = 0
0064  * \endverbatim
0065  * SQ:
0066  * \verbatim
0067    a x^2 + (a + \epsilon) y^2 + e x + f y + h = 0
0068  * \endverbatim
0069  * Normalized (\c epsilon = 0)
0070  * \verbatim
0071    x^2 + y^2 + e/b x + f/b y + g/b z  + h/b = 0
0072  * \endverbatim
0073  * Match terms:
0074  * \verbatim
0075     -2 x_0 = e/b --> x_0 = e / (-2 * b)
0076     -2 y_0 = f/b --> y_0 = f / (-2 * b)
0077    (x_0^2 + y_0^2 - r^2) = h/b
0078       -> r^2 = x_0^2 + y_0^2 - h/b
0079  * \endverbatim
0080  */
0081 template<Axis T>
0082 std::optional<CylAligned<T>>
0083 QuadricCylConverter::operator()(AxisTag<T>, SimpleQuadric const& sq) const
0084 {
0085     // Other coordinate system
0086     constexpr auto U = CylAligned<T>::u_axis();
0087     constexpr auto V = CylAligned<T>::v_axis();
0088
0089     auto second = sq.second();
0090     if (!soft_equal_(0, second[to_int(T)]))
0091     {
0092         // Not the zero component we're looking for
0093         return {};
0094     }
0095     if (!soft_equal_(0, sq.first()[to_int(T)]))
0096     {
0097         // Not an axis-aligned cylinder
0098         return {};
0099     }
0100     if (!soft_equal_(second[to_int(U)], second[to_int(V)]))
0101     {
0102         // Not a *circular* cylinder
0103         return {};
0104     }
0105     CELER_ASSERT(second[to_int(U)] > 0 && second[to_int(V)] > 0);
0106
0107     // Normalize so U, V second-order coefficients are 1
0108     auto const inv_norm = 2 / (second[to_int(U)] + second[to_int(V)]);
0109
0110     // Calculate origin from first-order coefficients
0111     Real3 origin{0, 0, 0};
0112     origin[to_int(U)] = real_type{-0.5} * inv_norm * sq.first()[to_int(U)];
0113     origin[to_int(V)] = real_type{-0.5} * inv_norm * sq.first()[to_int(V)];
0114
0115     real_type radius_sq = ipow<2>(origin[to_int(U)])
0116                           + ipow<2>(origin[to_int(V)]) - sq.zeroth() * inv_norm;
0117
0118     if (radius_sq <= 0)
0119     {
0120         // No real solution
0121         return {};
0122     }
0123
0124     // Clear potential signed zeros before returning
0125     origin += real_type{0};
0126     return CylAligned<T>::from_radius_sq(origin, radius_sq);
0127 }
0128
0129 //---------------------------------------------------------------------------//
0130 }  // namespace detail
0131 }  // namespace celeritas