surf/detail/QuadricSphereConverter.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/QuadricSphereConverter.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <optional>
0011
0012 #include "corecel/math/SoftEqual.hh"
0013
0014 #include "../Plane.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 sphere.
0024  *
0025  * The simple quadric must have already undergone normalization (zero or one
0026  * negative signs in the second order terms).
0027  */
0028 class QuadricSphereConverter
0029 {
0030   public:
0031     // Construct with tolerance
0032     inline QuadricSphereConverter(real_type tol);
0033
0034     // Try converting to a sphere
0035     std::optional<Sphere> operator()(SimpleQuadric const& sq) const;
0036
0037   private:
0038     SoftEqual<> soft_equal_;
0039 };
0040
0041 //---------------------------------------------------------------------------//
0042 // INLINE DEFINITIONS
0043 //---------------------------------------------------------------------------//
0044 /*!
0045  * Construct with tolerance.
0046  */
0047 QuadricSphereConverter::QuadricSphereConverter(real_type tol)
0048     : soft_equal_{tol}
0049 {
0050 }
0051
0052 //---------------------------------------------------------------------------//
0053 /*!
0054  * Try converting to a sphere with this orientation.
0055  *
0056  * \verbatim
0057    x^2 + y^2 + z^2
0058        - 2x_0 - 2y_0 - 2z_0
0059        + x_0^2 + y_0^2 + z_0^2 - r^2 = 0
0060  * \endverbatim
0061  * SQ:
0062  * \verbatim
0063    a x^2 + a y^2 + a z^2 + e x + f y + g z  + h = 0
0064  * \endverbatim
0065  * so
0066  * \verbatim
0067     -2 x_0 = e/b --> x_0 = e / (-2 * a)
0068     -2 y_0 = f/b --> y_0 = f / (-2 * a)
0069     -2 z_0 = g/b --> z_0 = g / (-2 * a)
0070    x_0^2 + y_0^2 + z_0^2 - r^2 = h/b --> r^2 = (origin . origin - h/b)
0071  * \endverbatim
0072  *
0073  * All the coefficients should be positive or nearly so.
0074  */
0075 std::optional<Sphere>
0076 QuadricSphereConverter::operator()(SimpleQuadric const& sq) const
0077 {
0078     CELER_EXPECT(std::all_of(
0079         sq.second().begin(), sq.second().end(), [this](real_type v) {
0080             return v >= 0 || soft_equal_(v, 0);
0081         }));
0082     constexpr auto X = to_int(Axis::x);
0083     constexpr auto Y = to_int(Axis::y);
0084     constexpr auto Z = to_int(Axis::z);
0085
0086     auto second = sq.second();
0087     if (!soft_equal_(second[X], second[Y])
0088         || !soft_equal_(second[X], second[Z]))
0089     {
0090         // Coefficients aren't equal: it's some sort of ellipsoid
0091         return {};
0092     }
0093
0094     real_type const inv_norm = 3 / (second[0] + second[1] + second[2]);
0095     CELER_ASSERT(inv_norm > 0);
0096
0097     Real3 origin = make_array(sq.first());
0098     origin *= real_type{-0.5} * inv_norm;
0099
0100     real_type radius_sq = dot_product(origin, origin) - sq.zeroth() * inv_norm;
0101     if (radius_sq <= 0)
0102     {
0103         // No real solution
0104         return {};
0105     }
0106
0107     // Clear potential signed zeros before returning
0108     origin += real_type{0};
0109     return Sphere::from_radius_sq(origin, radius_sq);
0110 }
0111
0112 //---------------------------------------------------------------------------//
0113 }  // namespace detail
0114 }  // namespace celeritas