orange/transform/Transformation.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/Transformation.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <algorithm>
0010
0011 #include "corecel/cont/Span.hh"
0012 #include "corecel/math/ArrayOperators.hh"
0013 #include "geocel/Types.hh"
0014 #include "orange/MatrixUtils.hh"
0015 #include "orange/OrangeTypes.hh"
0016
0017 namespace celeritas
0018 {
0019 class Translation;
0020 class SignedPermutation;
0021
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Apply transformations with rotation and/or reflection.
0025  *
0026  * \note The nomenclature in this class assumes the translation vector and
0027  * rotation matrix given represent "daughter-to-parent"! This is because we
0028  * think of rotations being with respect to the daughter's origin rather than
0029  * the parent's.
0030  *
0031  * This class enables transforms between daughter and parent coordinate
0032  * system. The transfer from a daughter into a parent system ("up" in a
0033  * hierarchy of universes) is
0034  * \f[
0035    \mathbf{r}_p = \mathbf{R}\mathbf{r}_d + \mathbf{t}\:,
0036  * \f]
0037  * Where the subscripts \e p,d refer to the parent and daughter coordinate
0038  * systems, respectively.  The vector \b t is a translation vector.  To go
0039  * from the parent into the daughter system ("down" in a universe hierarchy) we
0040  * apply the inverse:
0041  * \f[
0042    \mathbf{r}_d = \mathbf{R}^T(\mathbf{r}_p - \mathbf{t})\:.
0043  * \f]
0044  * where the transpose of \b R is equal to its inverse because the matrix is
0045  * unitary.
0046  *
0047  * The rotation matrix is indexed with C ordering, [i][j]. If a rotation
0048  * matrix, it should be a orthonormal with a determinant is 1 if not reflecting
0049  * (proper) or -1 if reflecting (improper). A transformation that applies a
0050  * scaling has non-unit eigenvalues.
0051  *
0052  * It is the caller's job to ensure a user-provided low-precision rotation
0053  * matrix is orthonormal: see \c celeritas::orthonormalize . (Add \c
0054  * CELER_VALIDATE to the calling code if constructing a transformation matrix
0055  * from user input or a suspect source.)
0056  *
0057  * \todo Scaling is not yet implemented correctly.
0058  */
0059 class Transformation
0060 {
0061   public:
0062     //@{
0063     //! \name Type aliases
0064     using StorageSpan = Span<real_type const, 12>;
0065     using Mat3 = SquareMatrixReal3;
0066     //@}
0067
0068     //! Calculated properties about the transformation
0069     struct Properties
0070     {
0071         bool reflects{false};  //!< Improper: applies a reflection
0072         bool scales{false};  //!< Applies a scale factor
0073     };
0074
0075     //! Transformation type identifier
0076     static CELER_CONSTEXPR_FUNCTION TransformType transform_type()
0077     {
0078         return TransformType::transformation;
0079     }
0080
0081   public:
0082     // Construct by inverting a parent-to-daughter transformation
0083     static Transformation from_inverse(Mat3 const& rot, Real3 const& trans);
0084
0085     //// CONSTRUCTORS ////
0086
0087     // Construct and check the input
0088     Transformation(Mat3 const& rot, Real3 const& trans);
0089
0090     // Construct from an identity transformation + translation
0091     Transformation();
0092
0093     // Promote from a translation
0094     explicit Transformation(Translation const&);
0095
0096     // Promote from a signed permutation
0097     explicit Transformation(SignedPermutation const&);
0098
0099     // Construct inline from storage
0100     explicit inline CELER_FUNCTION Transformation(StorageSpan);
0101
0102     //// ACCESSORS ////
0103
0104     //! Rotation matrix
0105     CELER_FORCEINLINE_FUNCTION Mat3 const& rotation() const { return rot_; }
0106
0107     //! Translation vector
0108     CELER_FORCEINLINE_FUNCTION Real3 const& translation() const
0109     {
0110         return tra_;
0111     }
0112
0113     //! Get a view to the data for type-deleted storage
0114     CELER_FUNCTION StorageSpan data() const { return {&rot_[0][0], 12}; }
0115
0116     //// CALCULATION ////
0117
0118     // Transform from daughter to parent
0119     [[nodiscard]] inline CELER_FUNCTION Real3
0120     transform_up(Real3 const& pos) const;
0121
0122     // Transform from parent to daughter
0123     [[nodiscard]] inline CELER_FUNCTION Real3
0124     transform_down(Real3 const& parent_pos) const;
0125
0126     // Rotate from daughter to parent
0127     [[nodiscard]] inline CELER_FUNCTION Real3 rotate_up(Real3 const& dir) const;
0128
0129     // Rotate from parent to daughter
0130     [[nodiscard]] inline CELER_FUNCTION Real3
0131     rotate_down(Real3 const& parent_dir) const;
0132
0133     // Calculate the inverse during preprocessing
0134     Transformation calc_inverse() const;
0135
0136     // Calculate properties about the matrix
0137     Properties calc_properties() const;
0138
0139   private:
0140     Mat3 rot_;
0141     Real3 tra_;
0142 };
0143
0144 //---------------------------------------------------------------------------//
0145 //!@{
0146 //! Host-only comparators
0147 inline bool operator==(Transformation const& a, Transformation const& b)
0148 {
0149     auto a_data = a.data();
0150     return std::equal(a_data.begin(), a_data.end(), b.data().begin());
0151 }
0152
0153 inline bool operator!=(Transformation const& a, Transformation const& b)
0154 {
0155     return !(a == b);
0156 }
0157 //!@}
0158
0159 //---------------------------------------------------------------------------//
0160 // INLINE DEFINITIONS
0161 //---------------------------------------------------------------------------//
0162 /*!
0163  * Construct inline from storage.
0164  */
0165 CELER_FUNCTION Transformation::Transformation(StorageSpan s)
0166     : rot_{Real3{s[0], s[1], s[2]},
0167            Real3{s[3], s[4], s[5]},
0168            Real3{s[6], s[7], s[8]}}
0169     , tra_{s[9], s[10], s[11]}
0170 {
0171 }
0172
0173 //---------------------------------------------------------------------------//
0174 /*!
0175  * Transform from daughter to parent.
0176  *
0177  * Apply the rotation matrix, add the translation.
0178  */
0179 CELER_FUNCTION Real3 Transformation::transform_up(Real3 const& pos) const
0180 {
0181     return gemv(real_type{1}, rot_, pos, real_type{1}, tra_);
0182 }
0183
0184 //---------------------------------------------------------------------------//
0185 /*!
0186  * Transform from parent to daughter.
0187  *
0188  * Subtract the translation, then apply the inverse of the rotation matrix (its
0189  * transpose).
0190  */
0191 CELER_FUNCTION Real3 Transformation::transform_down(Real3 const& pos) const
0192 {
0193     return gemv(matrix::transpose, rot_, pos - tra_);
0194 }
0195
0196 //---------------------------------------------------------------------------//
0197 /*!
0198  * Rotate from daughter to parent.
0199  */
0200 CELER_FUNCTION Real3 Transformation::rotate_up(Real3 const& d) const
0201 {
0202     return gemv(rot_, d);
0203 }
0204
0205 //---------------------------------------------------------------------------//
0206 /*!
0207  * Rotate from parent to daughter.
0208  */
0209 CELER_FUNCTION Real3 Transformation::rotate_down(Real3 const& d) const
0210 {
0211     return gemv(matrix::transpose, rot_, d);
0212 }
0213
0214 //---------------------------------------------------------------------------//
0215 }  // namespace celeritas