orange/transform/SignedPermutation.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2023 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/transform/SignedPermutation.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/cont/EnumArray.hh"
0011 #include "corecel/cont/Range.hh"
0012 #include "corecel/cont/Span.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/math/Turn.hh"
0015 #include "geocel/Types.hh"
0016 #include "orange/OrangeTypes.hh"
0017
0018 namespace celeritas
0019 {
0020 //---------------------------------------------------------------------------//
0021 /*!
0022  * Apply a rotation that remaps and possibly flips signs.
0023  *
0024  * A signed permutation matrix is a special matrix that has only one entry with
0025  * the value of \f$\pm1\f$ in each row and each column. This is a specialized
0026  * rotation matrix that, when applied to a vector, simply exchanges the
0027  * locations and/or flips the signs of the vector entries.
0028  *
0029  * This class stores a special version of the daughter-to-parent rotation
0030  * matrix:
0031  * \f[
0032  \mathbf{R} = \begin{bmatrix}
0033   \mathbf{e}_x \\ \hline
0034   \mathbf{e}_y \\ \hline
0035   \mathbf{e}_z
0036   \end{bmatrix}
0037   \f]
0038  * where \f$ \mathbf{e}_u \f$ has exactly one entry with a value \f$ \pm 1
0039  \f$ and the other entries are zero.
0040  *
0041  * The underlying storage are a compressed series of bits in little-endian form
0042  * that indicate the positions of the nonzero entry followed by the sign:
0043  * \verbatim
0044    [flip z'][z' axis][flip y'][y' axis][flip x'][x' axis]
0045          8   7    6        5   4     3        2   1   0  bit position
0046  * \endverbatim
0047  * This allows the "rotate up" to simply copy one value at a time into a new
0048  * position, and optionally flip the sign of the result.
0049  *
0050  * Construction of this class takes a length 3 array of \c SignedAxis values.
0051  * The sign is a '+' or '-' character and the axis is the position of the
0052  * nonzero component in that row.
0053  *
0054  * Note that the unsigned integer type (UIntT) is at least 16 bits, which is
0055  * sufficient accurately round-trip through a floating point value.
0056  */
0057 class SignedPermutation
0058 {
0059   public:
0060     //!@{
0061     //! \name Type aliases
0062     using SignedAxis = std::pair<char, Axis>;
0063     using SignedAxes = EnumArray<Axis, SignedAxis>;
0064     using StorageSpan = Span<real_type const, 1>;
0065     using DataArray = Array<real_type, 1>;
0066     using UIntT = short unsigned int;
0067     //!@}
0068
0069   public:
0070     // Construct with an identity permutation
0071     SignedPermutation();
0072
0073     // Construct with the permutation vector
0074     explicit SignedPermutation(SignedAxes permutation);
0075
0076     // Construct inline from storage
0077     explicit inline CELER_FUNCTION SignedPermutation(StorageSpan s);
0078
0079     //// ACCESSORS ////
0080
0081     //! Opaque type for hashing and comparison
0082     UIntT value() const { return compressed_; }
0083
0084     // Reconstruct the permutation
0085     SignedAxes permutation() const;
0086
0087     // Get a view to the data for type-deleted storage
0088     DataArray data() const;
0089
0090     //// CALCULATION ////
0091
0092     // Transform from daughter to parent
0093     [[nodiscard]] inline CELER_FUNCTION Real3
0094     transform_up(Real3 const& pos) const;
0095
0096     // Transform from parent to daughter
0097     [[nodiscard]] inline CELER_FUNCTION Real3
0098     transform_down(Real3 const& parent_pos) const;
0099
0100     // Rotate from daughter to parent
0101     [[nodiscard]] inline CELER_FUNCTION Real3 rotate_up(Real3 const& dir) const;
0102
0103     // Rotate from parent to daughter
0104     [[nodiscard]] inline CELER_FUNCTION Real3
0105     rotate_down(Real3 const& parent_dir) const;
0106
0107   private:
0108     //// DATA ////
0109
0110     UIntT compressed_;
0111
0112     //// FUNCTIONS ////
0113
0114     // Maximum compressed integer value used for bounds checking
0115     static CELER_CONSTEXPR_FUNCTION UIntT max_value() { return (1 << 9) - 1; }
0116 };
0117
0118 //---------------------------------------------------------------------------//
0119 // FREE FUNCTIONS
0120 //---------------------------------------------------------------------------//
0121
0122 // Make a permutation by rotating about the given axis
0123 SignedPermutation make_permutation(Axis ax, QuarterTurn qtr);
0124
0125 //!@{
0126 //! Host-only comparators
0127 inline bool operator==(SignedPermutation const& a, SignedPermutation const& b)
0128 {
0129     return a.value() == b.value();
0130 }
0131
0132 inline bool operator!=(SignedPermutation const& a, SignedPermutation const& b)
0133 {
0134     return !(a == b);
0135 }
0136 //!@}
0137
0138 //---------------------------------------------------------------------------//
0139 // INLINE DEFINITIONS
0140 //---------------------------------------------------------------------------//
0141 /*!
0142  * Construct inline from storage.
0143  */
0144 CELER_FUNCTION SignedPermutation::SignedPermutation(StorageSpan s)
0145     : compressed_{static_cast<UIntT>(s[0])}
0146 {
0147     CELER_EXPECT(s[0] >= 0 && s[0] <= static_cast<real_type>(max_value()));
0148 }
0149
0150 //---------------------------------------------------------------------------//
0151 /*!
0152  * Transform from daughter to parent.
0153  */
0154 CELER_FORCEINLINE_FUNCTION Real3
0155 SignedPermutation::transform_up(Real3 const& pos) const
0156 {
0157     return this->rotate_up(pos);
0158 }
0159
0160 //---------------------------------------------------------------------------//
0161 /*!
0162  * Transform from parent to daughter.
0163  */
0164 CELER_FORCEINLINE_FUNCTION Real3
0165 SignedPermutation::transform_down(Real3 const& parent_pos) const
0166 {
0167     return this->rotate_down(parent_pos);
0168 }
0169
0170 //---------------------------------------------------------------------------//
0171 /*!
0172  * Rotate from daughter to parent.
0173  */
0174 CELER_FORCEINLINE_FUNCTION Real3 SignedPermutation::rotate_up(Real3 const& d) const
0175 {
0176     Real3 result;
0177
0178     UIntT temp = compressed_;
0179     for (auto ax : range(Axis::size_))
0180     {
0181         // Copy to new axis
0182         unsigned int new_ax = temp & 0b11;
0183         result[to_int(ax)] = d[new_ax];
0184         temp >>= 2;
0185         if (temp & 0b1)
0186         {
0187             // Flip the bit, avoiding signed zero
0188             result[to_int(ax)] = negate(result[to_int(ax)]);
0189         }
0190         temp >>= 1;
0191     }
0192
0193     return result;
0194 }
0195
0196 //---------------------------------------------------------------------------//
0197 /*!
0198  * Rotate from parent to daughter.
0199  */
0200 CELER_FORCEINLINE_FUNCTION Real3
0201 SignedPermutation::rotate_down(Real3 const& d) const
0202 {
0203     Real3 result;
0204
0205     UIntT temp = compressed_;
0206     for (auto ax : range(Axis::size_))
0207     {
0208         // Copy to new axis
0209         unsigned int new_ax = temp & 0b11;
0210         result[new_ax] = d[to_int(ax)];
0211         temp >>= 2;
0212         if (temp & 0b1)
0213         {
0214             // Flip the bit, avoiding signed zero
0215             result[new_ax] = negate(result[new_ax]);
0216         }
0217         temp >>= 1;
0218     }
0219
0220     return result;
0221 }
0222
0223 //---------------------------------------------------------------------------//
0224 }  // namespace celeritas