surf/detail/QuadricPlaneConverter.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/QuadricPlaneConverter.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <algorithm>
0011
0012 #include "corecel/Assert.hh"
0013 #include "corecel/math/SoftEqual.hh"
0014
0015 #include "../Plane.hh"
0016 #include "../SimpleQuadric.hh"
0017
0018 namespace celeritas
0019 {
0020 namespace detail
0021 {
0022 //-------------------------------------------------------------------------//
0023 /*!
0024  * Convert a simple quadric to a plane.
0025  */
0026 class QuadricPlaneConverter
0027 {
0028   public:
0029     // Construct with tolerance
0030     inline QuadricPlaneConverter(real_type tol);
0031
0032     // Convert to a plane
0033     Plane operator()(SimpleQuadric const& sq) const;
0034
0035   private:
0036     SoftZero<> soft_zero_;
0037 };
0038
0039 //---------------------------------------------------------------------------//
0040 // INLINE DEFINITIONS
0041 //---------------------------------------------------------------------------//
0042 /*!
0043  * Construct with tolerance.
0044  */
0045 QuadricPlaneConverter::QuadricPlaneConverter(real_type tol) : soft_zero_{tol}
0046 {
0047 }
0048
0049 //---------------------------------------------------------------------------//
0050 /*!
0051  * Convert to a plane.
0052  *
0053  * Simple quadric is calculated as
0054  * dx + ey + fz + g = 0, but plane is
0055  * ax + bx + cz - d = 0
0056  * so we need to reverse the sign of the scalar component. We also need to
0057  * normalize with the same factor.
0058  */
0059 Plane QuadricPlaneConverter::operator()(SimpleQuadric const& sq) const
0060 {
0061     CELER_EXPECT(
0062         std::all_of(sq.second().begin(), sq.second().end(), soft_zero_));
0063     CELER_EXPECT(
0064         !std::all_of(sq.first().begin(), sq.first().end(), soft_zero_));
0065     // Second-order coefficients are zero: return a plane
0066     auto n = make_array(sq.first());
0067
0068     real_type norm_factor = 1 / celeritas::norm(n);
0069     n *= norm_factor;
0070
0071     real_type d = -sq.zeroth() * norm_factor;
0072
0073     return Plane{n, d};
0074 }
0075
0076 //---------------------------------------------------------------------------//
0077 }  // namespace detail
0078 }  // namespace celeritas