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