corecel/math/ArraySoftUnit.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 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 corecel/math/ArraySoftUnit.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/Types.hh"
0013 #include "corecel/cont/Array.hh"
0014
0015 #include "detail/SoftEqualTraits.hh"
0016
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Test for being approximately a unit vector.
0022  *
0023  * Consider a unit vector \em v with a small perturbation along a unit vector
0024  * \em e : \f[
0025    \vec v + \epsilon \vec e
0026   \f]
0027  * The magnitude squared is
0028  * \f[
0029   m^2 = (v + \epsilon e) \cdot (v + \epsilon e)
0030    = v \cdot v + 2 \epsilon v \cdot e +  \epsilon^2 e \cdot e
0031    = 1 + 2 \epsilon v \cdot e + \epsilon^2
0032  \f]
0033  *
0034  * Since \f[ |v \cdot e|  <= |v||e| = 1 \f] by the triangle inequality,
0035  * then the magnitude squared of a perturbed unit vector is bounded
0036  * \f[
0037   m^2 = 1 \pm 2 \epsilon + \epsilon^2
0038   \f]
0039  *
0040  * Instead of calculating the square of the tolerance we loosely bound with
0041  * another epsilon.
0042  */
0043 template<class T = ::celeritas::real_type>
0044 class ArraySoftUnit
0045 {
0046   public:
0047     //!@{
0048     //! \name Type aliases
0049     using value_type = T;
0050     //!@}
0051
0052   public:
0053     // Construct with explicit tolerance
0054     CELER_FUNCTION inline ArraySoftUnit(value_type tol);
0055
0056     // Construct with default tolerance
0057     CELER_CONSTEXPR_FUNCTION ArraySoftUnit();
0058
0059     // Calculate whether the array is nearly a unit vector
0060     template<::celeritas::size_type N>
0061     CELER_CONSTEXPR_FUNCTION bool operator()(Array<T, N> const& arr) const;
0062
0063   private:
0064     value_type tol_;
0065 };
0066
0067 //---------------------------------------------------------------------------//
0068 // TEMPLATE DEDUCTION GUIDES
0069 //---------------------------------------------------------------------------//
0070 template<class T>
0071 CELER_FUNCTION ArraySoftUnit(T) -> ArraySoftUnit<T>;
0072
0073 //---------------------------------------------------------------------------//
0074 // FREE FUNCTIONS
0075 //---------------------------------------------------------------------------//
0076 // Test for being approximately a unit vector
0077 template<class T, size_type N>
0078 CELER_CONSTEXPR_FUNCTION bool is_soft_unit_vector(Array<T, N> const& v);
0079
0080 //---------------------------------------------------------------------------//
0081 // INLINE DEFINITIONS
0082 //---------------------------------------------------------------------------//
0083 /*!
0084  * Construct with explicit tolereance.
0085  */
0086 template<class T>
0087 CELER_FUNCTION ArraySoftUnit<T>::ArraySoftUnit(T tol) : tol_{3 * tol}
0088 {
0089     CELER_EXPECT(tol_ > 0);
0090 }
0091
0092 //---------------------------------------------------------------------------//
0093 /*!
0094  * Construct with default tolereance.
0095  */
0096 template<class T>
0097 CELER_CONSTEXPR_FUNCTION ArraySoftUnit<T>::ArraySoftUnit()
0098     : tol_{3 * detail::SoftEqualTraits<T>::rel_prec()}
0099 {
0100 }
0101
0102 //---------------------------------------------------------------------------//
0103 /*!
0104  * Calculate whether the array is nearly a unit vector.
0105  *
0106  * The calculation below is equivalent to
0107  * \code
0108  * return SoftEqual{tol, tol}(1, dot_product(arr, arr));
0109  * \endcode.
0110  */
0111 template<class T>
0112 template<::celeritas::size_type N>
0113 CELER_CONSTEXPR_FUNCTION bool
0114 ArraySoftUnit<T>::operator()(Array<T, N> const& arr) const
0115 {
0116     T length_sq{};
0117     for (size_type i = 0; i != N; ++i)
0118     {
0119         length_sq = std::fma(arr[i], arr[i], length_sq);
0120     }
0121     return std::fabs(length_sq - 1) < tol_ * std::fmax(real_type(1), length_sq);
0122 }
0123
0124 //---------------------------------------------------------------------------//
0125 //! Test with default tolerance for being a unit vector
0126 template<class T, size_type N>
0127 CELER_CONSTEXPR_FUNCTION bool is_soft_unit_vector(Array<T, N> const& v)
0128 {
0129     return ArraySoftUnit<T>{}(v);
0130 }
0131
0132 //---------------------------------------------------------------------------//
0133 }  // namespace celeritas