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