diff --git a/src/corecel/cont/Span.hh b/src/corecel/cont/Span.hh index b28ac04074..850f8a4d74 100644 --- a/src/corecel/cont/Span.hh +++ b/src/corecel/cont/Span.hh @@ -177,7 +177,7 @@ template constexpr std::size_t Span::extent; //---------------------------------------------------------------------------// -// HELPER FUNCTIONS +// FREE FUNCTIONS //---------------------------------------------------------------------------// //! Get a mutable fixed-size view to an array template @@ -218,5 +218,19 @@ CELER_FUNCTION Span make_span(T const& cont) return {cont.data(), cont.size()}; } +//---------------------------------------------------------------------------// +//! Construct an array from a fixed-size span +template +CELER_CONSTEXPR_FUNCTION Array, N> +make_array(Span const& s) +{ + Array, N> result; + for (std::size_t i = 0; i < N; ++i) + { + result[i] = s[i]; + } + return result; +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/orange/CMakeLists.txt b/src/orange/CMakeLists.txt index 5cd05bb7bf..89a6879ca8 100644 --- a/src/orange/CMakeLists.txt +++ b/src/orange/CMakeLists.txt @@ -22,6 +22,7 @@ list(APPEND SOURCES detail/RectArrayInserter.cc detail/UnitInserter.cc surf/SurfaceIO.cc + transform/Transformation.cc ) if(CELERITAS_USE_JSON) diff --git a/src/orange/MatrixUtils.cc b/src/orange/MatrixUtils.cc index 9e9025c6dd..78c6550294 100644 --- a/src/orange/MatrixUtils.cc +++ b/src/orange/MatrixUtils.cc @@ -7,7 +7,10 @@ //---------------------------------------------------------------------------// #include "MatrixUtils.hh" +#include + #include "corecel/math/ArrayUtils.hh" +#include "corecel/math/SoftEqual.hh" using Mat3 = celeritas::SquareMatrixReal3; @@ -66,11 +69,12 @@ gemm(SquareMatrix const& a, SquareMatrix const& b) /*! * Normalize and orthogonalize a small, dense matrix. * - * \return normalization (zero if already orthonormal) - * * This is used for constructing rotation matrices from user-given matrices * that may only have a few digits of precision (e.g. were read from an XML * file). It uses the modified Gram-Schmidt orthogonalization algorithm. + * + * If debug assertions are enabled, the normality of the resulting matrix will + * be checked. A singular matrix will fail. */ template void orthonormalize(SquareMatrix* mat) @@ -95,19 +99,8 @@ void orthonormalize(SquareMatrix* mat) } } - // Check result for orthonormality: if the matrix is singular this will - // fail - if (CELERITAS_DEBUG) - { - for (size_type i = 0; i != N; ++i) - { - for (size_type ip = 0; ip != i; ++ip) - { - CELER_ENSURE(soft_equal(dot_product((*mat)[i], (*mat)[ip]), - static_cast(i == ip ? 1 : 0))); - } - } - } + // Check result for orthonormality + CELER_ENSURE(soft_equal(std::fabs(determinant(*mat)), T{1})); } //---------------------------------------------------------------------------// diff --git a/src/orange/MatrixUtils.hh b/src/orange/MatrixUtils.hh index 26674e9d65..2ff2986453 100644 --- a/src/orange/MatrixUtils.hh +++ b/src/orange/MatrixUtils.hh @@ -15,6 +15,18 @@ namespace celeritas { +namespace matrix +{ +//---------------------------------------------------------------------------// +//!@{ +//! Policy tags +struct TransposePolicy +{ +}; +inline constexpr TransposePolicy transpose{}; +//!@} +} // namespace matrix + //---------------------------------------------------------------------------// // Apply a matrix to an array template @@ -25,7 +37,18 @@ inline CELER_FUNCTION Array gemv(T alpha, Array const& y); //---------------------------------------------------------------------------// -//! Apply a matrix to an array without scaling or addition +// Apply the transpose of a matrix to an array +template +inline CELER_FUNCTION Array gemv(matrix::TransposePolicy, + T alpha, + SquareMatrix const& a, + Array const& x, + T beta, + Array const& y); + +//---------------------------------------------------------------------------// +//!@{ +//! Apply a matrix or its transpose to an array, without scaling or addition template inline CELER_FUNCTION Array gemv(SquareMatrix const& a, Array const& x) @@ -33,6 +56,13 @@ gemv(SquareMatrix const& a, Array const& x) return gemv(T{1}, a, x, T{0}, x); } +template +inline CELER_FUNCTION Array +gemv(matrix::TransposePolicy, SquareMatrix const& a, Array const& x) +{ + return gemv(matrix::transpose, T{1}, a, x, T{0}, x); +} +//!@} //---------------------------------------------------------------------------// // Host-only declarations // (double and float (and some int) for N=3 are instantiated in MatrixUtils.cc) @@ -67,7 +97,7 @@ SquareMatrixReal3 make_rotation(Axis ax, Turn rev, SquareMatrixReal3 const&); * z \gets \alpha A x + \beta y * \f] * - * This should be equivalent to BLAS' GEMV without the option to transpose. All + * This should be equivalent to BLAS' GEMV without transposition. All * matrix orderings are C-style: mat[i][j] is for row i, column j . * * \warning This implementation is limited and slow. @@ -91,5 +121,41 @@ CELER_FUNCTION Array gemv(T alpha, return result; } +//---------------------------------------------------------------------------// +/*! + * Naive transposed generalized matrix-vector multiply. + * + * \f[ + * z \gets \alpha A^T x + \beta y + * \f] + * + * This should be equivalent to BLAS' GEMV with the 't' option. All + * matrix orderings are C-style: mat[i][j] is for row i, column j . + * + * \warning This implementation is limited and slow. + */ +template +CELER_FUNCTION Array gemv(matrix::TransposePolicy, + T alpha, + SquareMatrix const& a, + Array const& x, + T beta, + Array const& y) +{ + Array result; + for (size_type i = 0; i != N; ++i) + { + result[i] = beta * y[i]; + } + for (size_type j = 0; j != N; ++j) + { + for (size_type i = 0; i != N; ++i) + { + result[i] += alpha * (a[j][i] * x[j]); + } + } + return result; +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/orange/OrangeData.hh b/src/orange/OrangeData.hh index 0c255a44b4..18823967e2 100644 --- a/src/orange/OrangeData.hh +++ b/src/orange/OrangeData.hh @@ -283,7 +283,7 @@ struct OrangeParamsData Items bih_leaf_nodes; Items daughters; - Items translations; + Items translations; UniverseIndexerData universe_indexer_data; diff --git a/src/orange/OrangeTrackView.hh b/src/orange/OrangeTrackView.hh index 7a887858bc..874d69d805 100644 --- a/src/orange/OrangeTrackView.hh +++ b/src/orange/OrangeTrackView.hh @@ -14,9 +14,9 @@ #include "OrangeData.hh" #include "OrangeTypes.hh" -#include "Translator.hh" #include "detail/LevelStateAccessor.hh" #include "detail/UniverseIndexer.hh" +#include "transform/Translation.hh" #include "univ/SimpleUnitTracker.hh" #include "univ/UniverseTypeTraits.hh" #include "univ/detail/Types.hh" @@ -263,9 +263,8 @@ OrangeTrackView::operator=(Initializer_t const& init) if (daughter_id) { auto const& daughter = params_.daughters[daughter_id]; - auto const& trans = params_.translations[daughter.translation_id]; - TranslatorDown td(trans); - local.pos = td(local.pos); + Translation trans{params_.translations[daughter.translation_id]}; + local.pos = trans.transform_down(local.pos); uid = daughter.universe_id; ++level; @@ -544,10 +543,8 @@ CELER_FUNCTION void OrangeTrackView::move_internal(Real3 const& pos) auto daughter_id = tracker.daughter(lsa.vol()); CELER_ASSERT(daughter_id); auto const& daughter = params_.daughters[daughter_id]; - auto const& trans = params_.translations[daughter.translation_id]; - - TranslatorDown td(trans); - local_pos = td(pos); + Translation trans{params_.translations[daughter.translation_id]}; + local_pos = trans.transform_down(pos); } } @@ -625,8 +622,7 @@ CELER_FUNCTION void OrangeTrackView::cross_boundary() auto daughter = params_.daughters[daughter_id]; // Get the translator at the parent level, in order to translate into // daughter - TranslatorDown translator( - params_.translations[daughter.translation_id]); + Translation trans{params_.translations[daughter.translation_id]}; // Make the current level the daughter level ++level; @@ -634,7 +630,7 @@ CELER_FUNCTION void OrangeTrackView::cross_boundary() auto tracker = this->make_tracker(universe_id); // Create local state on the daughter level - local.pos = translator(local.pos); + local.pos = trans.transform_down(local.pos); local.volume = {}; local.surface = {}; local.temp_sense = this->make_temp_sense(); diff --git a/src/orange/OrangeTypes.hh b/src/orange/OrangeTypes.hh index 44b8aac948..1776a265e0 100644 --- a/src/orange/OrangeTypes.hh +++ b/src/orange/OrangeTypes.hh @@ -64,11 +64,8 @@ using SimpleUnitId = OpaqueId; //! Opaque index for rectilinear array data using RectArrayId = OpaqueId; -//! Translation of a single embedded universe -using Translation = Real3; - //! Identifier for a translation of a single embedded universe -using TranslationId = OpaqueId; +using TranslationId = OpaqueId; //! Identifier for a relocatable set of volumes using UniverseId = OpaqueId; @@ -124,6 +121,17 @@ enum class SurfaceType : unsigned char size_ //!< Sentinel value for number of surface types }; +//---------------------------------------------------------------------------// +/*! + * Enumeration for mapping transform implementations to integers. + */ +enum class TransformType : unsigned char +{ + translation, //!< Translation only + transformation, //!< Translation plus rotation + size_ +}; + //---------------------------------------------------------------------------// /*! * Enumeration for type-deleted universe storage. diff --git a/src/orange/Translator.hh b/src/orange/Translator.hh deleted file mode 100644 index e7c730035d..0000000000 --- a/src/orange/Translator.hh +++ /dev/null @@ -1,106 +0,0 @@ -//----------------------------------*-C++-*----------------------------------// -// Copyright 2021-2023 UT-Battelle, LLC, and other Celeritas developers. -// See the top-level COPYRIGHT file for details. -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file orange/Translator.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "corecel/cont/Range.hh" -#include "orange/OrangeTypes.hh" - -namespace celeritas -{ -//---------------------------------------------------------------------------// -/*! - * Translate points from a parent's reference frame into the daughter. - * - * The input "translation" is the transformation applied to a daughter universe - * to place it in the parent. The daughter is a "lower" level compared to the - * parent. - */ -class TranslatorDown -{ - public: - // Construct with the parent-to-daughter translation - explicit inline CELER_FUNCTION TranslatorDown(Real3 const& translation); - - // Translate a single point - CELER_FORCEINLINE_FUNCTION Real3 operator()(Real3 const& parent) const; - - private: - Real3 const& translation_; -}; - -//---------------------------------------------------------------------------// -/*! - * Translate points from a daughter's reference frame "up" into the parent. - * - * The input "translation" is the transformation applied to a daughter universe - * to place it in the parent. The daughter is a "lower" level compared to the - * parent. - */ -class TranslatorUp -{ - public: - // Construct with the parent-to-daughter translation - inline CELER_FUNCTION TranslatorUp(Real3 const& translation); - - // Translate a single point - CELER_FORCEINLINE_FUNCTION Real3 operator()(Real3 const& parent) const; - - private: - Real3 const& translation_; -}; - -//---------------------------------------------------------------------------// -// INLINE DEFINITIONS -//---------------------------------------------------------------------------// -/*! - * Construct with translation. - */ -CELER_FUNCTION TranslatorDown::TranslatorDown(Real3 const& translation) - : translation_(translation) -{ -} - -//---------------------------------------------------------------------------// -/*! - * Translate a single point. - */ -CELER_FUNCTION Real3 TranslatorDown::operator()(Real3 const& parent) const -{ - Real3 daughter; - for (int i : range(3)) - { - daughter[i] = parent[i] - translation_[i]; - } - return daughter; -} - -//---------------------------------------------------------------------------// -/*! - * Construct with translation. - */ -CELER_FUNCTION TranslatorUp::TranslatorUp(Real3 const& translation) - : translation_(translation) -{ -} - -//---------------------------------------------------------------------------// -/*! - * Translate a single point. - */ -CELER_FUNCTION Real3 TranslatorUp::operator()(Real3 const& parent) const -{ - Real3 daughter; - for (int i : range(3)) - { - daughter[i] = parent[i] + translation_[i]; - } - return daughter; -} - -//---------------------------------------------------------------------------// -} // namespace celeritas diff --git a/src/orange/construct/OrangeInput.hh b/src/orange/construct/OrangeInput.hh index 7a65b0420f..f012fecf77 100644 --- a/src/orange/construct/OrangeInput.hh +++ b/src/orange/construct/OrangeInput.hh @@ -16,7 +16,6 @@ #include "orange/BoundingBox.hh" #include "orange/OrangeData.hh" #include "orange/OrangeTypes.hh" -#include "orange/Translator.hh" namespace celeritas { @@ -81,7 +80,7 @@ struct VolumeInput struct DaughterInput { UniverseId universe_id; - Translation translation; + Real3 translation; }; //---------------------------------------------------------------------------// diff --git a/src/orange/construct/OrangeInputIO.json.cc b/src/orange/construct/OrangeInputIO.json.cc index 88b5ab842b..c84b74b13d 100644 --- a/src/orange/construct/OrangeInputIO.json.cc +++ b/src/orange/construct/OrangeInputIO.json.cc @@ -10,6 +10,7 @@ #include #include #include +#include #include #include @@ -26,21 +27,6 @@ #include "orange/OrangeTypes.hh" #include "orange/construct/OrangeInput.hh" -namespace -{ - -//---------------------------------------------------------------------------// -/*! - * Create a Translation object from a Span into a vector of translation data. - */ -celeritas::Translation -make_translation(celeritas::Span const& trans) -{ - CELER_EXPECT(trans.size() == 3); - return celeritas::Translation{trans[0], trans[1], trans[2]}; -} -} // namespace - namespace celeritas { namespace @@ -111,6 +97,20 @@ std::vector parse_logic(char const* c) return result; } +//---------------------------------------------------------------------------// +/*! + * Get the i'th slice of a span of data. + */ +template +decltype(auto) slice(Span data, size_type i) +{ + CELER_ASSERT(N * (i + 1) <= data.size()); + Array, N> result; + std::copy_n(data.data() + i * N, N, result.begin()); + return result; +} + +//---------------------------------------------------------------------------// } // namespace //---------------------------------------------------------------------------// @@ -209,10 +209,8 @@ void from_json(nlohmann::json const& j, UnitInput& value) UnitInput::MapVolumeDaughter daughter_map; for (auto i : range(parent_cells.size())) { - daughter_map[LocalVolumeId{parent_cells[i]}] - = {UniverseId{daughters[i]}, - make_translation( - Span(translations.data() + 3 * i, 3))}; + daughter_map[LocalVolumeId{parent_cells[i]}] = { + UniverseId{daughters[i]}, slice<3>(make_span(translations), i)}; } value.daughter_map = std::move(daughter_map); @@ -248,8 +246,7 @@ void from_json(nlohmann::json const& j, RectArrayInput& value) for (auto i : range(daughters.size())) { value.daughters.push_back({UniverseId{daughters[i]}, - make_translation(Span( - translations.data() + 3 * i, 3))}); + slice<3>(make_span(translations), i)}); } } } diff --git a/src/orange/transform/Transformation.cc b/src/orange/transform/Transformation.cc new file mode 100644 index 0000000000..77db14859f --- /dev/null +++ b/src/orange/transform/Transformation.cc @@ -0,0 +1,50 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file orange/transform/Transformation.cc +//---------------------------------------------------------------------------// +#include "Transformation.hh" + +#include + +#include "corecel/math/SoftEqual.hh" +#include "orange/MatrixUtils.hh" + +#include "Translation.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct as an identity transform. + */ +Transformation::Transformation() : Transformation{Translation{}} {} + +//---------------------------------------------------------------------------// +/*! + * Construct and check the input. + * + * The input rotation matrix should be an orthonormal matrix without + * reflecting, thus having a determinant of unity. It is the caller's job to + * ensure a user-provided low-precision matrix is orthonormalized. + */ +Transformation::Transformation(Mat3 const& rot, Real3 const& trans) + : rot_(rot), tra_(trans) +{ + CELER_EXPECT(soft_equal(determinant(rot_), real_type(1))); +} + +//---------------------------------------------------------------------------// +/*! + * Promote from a translation. + */ +Transformation::Transformation(Translation const& tr) + : rot_{Real3{1, 0, 0}, Real3{0, 1, 0}, Real3{0, 0, 1}} + , tra_{tr.translation()} +{ +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/transform/Transformation.hh b/src/orange/transform/Transformation.hh new file mode 100644 index 0000000000..f4f69389c5 --- /dev/null +++ b/src/orange/transform/Transformation.hh @@ -0,0 +1,163 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file orange/Transformation.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/cont/Span.hh" +#include "corecel/math/ArrayOperators.hh" +#include "orange/MatrixUtils.hh" +#include "orange/OrangeTypes.hh" +#include "orange/Types.hh" + +namespace celeritas +{ +class Translation; +//---------------------------------------------------------------------------// +/*! + * Apply transformations with rotation. + * + * \note The nomenclature in this class assumes the translation vector and + * rotation matrix given represent "daughter-to-parent"! This is because we + * think of rotations being with respect to the daughter's origin rather than + * the parent's. + * + * This class enables transforms between daughter and parent coordinate + * system. The transfer from a daughter into a parent system ("up" in a + * hierarchy of universes) is + * \f[ + \mathbf{r}_p = \mathbf{R}\mathbf{r}_d + \mathbf{t}\:, + * \f] + * Where the subscripts \e p,d refer to the parent and daughter coordinate + * systems, respectively. The vector \b t is a translation vector. To go + * from the parent into the daughter system ("down" in a universe hierarchy) we + * apply the inverse: + * \f[ + \mathbf{r}_d = \mathbf{R}^T(\mathbf{r}_p - \mathbf{t})\:. + * \f] + * where the transpose of \b R is equal to its inverse because the matrix is + * unitary. + * + * The rotation matrix is indexed with C ordering, [i][j]. + */ +class Transformation +{ + public: + //@{ + //! \name Type aliases + using StorageSpan = Span; + using Mat3 = SquareMatrixReal3; + //@} + + //! Transformation type identifier + static CELER_CONSTEXPR_FUNCTION TransformType transform_type() + { + return TransformType::transformation; + } + + public: + //// CONSTRUCTORS //// + + // Construct and check the input + Transformation(Mat3 const& rot, Real3 const& trans); + + // Construct from an identity transformation + translation + Transformation(); + + // Promote from a translation + explicit Transformation(Translation const&); + + // Construct inline from storage + explicit inline CELER_FUNCTION Transformation(StorageSpan); + + //// ACCESSORS //// + + //! Rotation matrix + CELER_FORCEINLINE_FUNCTION Mat3 const& rotation() const { return rot_; } + + //! Translation vector + CELER_FORCEINLINE_FUNCTION Real3 const& translation() const + { + return tra_; + } + + //! Get a view to the data for type-deleted storage + CELER_FUNCTION StorageSpan data() const { return {&rot_[0][0], 12}; } + + //// CALCULATION //// + + // Transform from daughter to parent + inline CELER_FUNCTION Real3 transform_up(Real3 const& pos) const; + + // Transform from parent to daughter + inline CELER_FUNCTION Real3 transform_down(Real3 const& parent_pos) const; + + // Rotate from daughter to parent + inline CELER_FUNCTION Real3 rotate_up(Real3 const& dir) const; + + // Rotate from parent to daughter + inline CELER_FUNCTION Real3 rotate_down(Real3 const& parent_dir) const; + + private: + Mat3 rot_; + Real3 tra_; +}; + +//---------------------------------------------------------------------------// +/*! + * Construct inline from storage. + */ +CELER_FUNCTION Transformation::Transformation(StorageSpan s) + : rot_{Real3{s[0], s[1], s[2]}, + Real3{s[3], s[4], s[5]}, + Real3{s[6], s[7], s[8]}} + , tra_{s[9], s[10], s[11]} +{ +} + +//---------------------------------------------------------------------------// +/*! + * Transform from daughter to parent. + * + * Apply the rotation matrix, add the translation. + */ +CELER_FUNCTION Real3 Transformation::transform_up(Real3 const& pos) const +{ + return gemv(real_type{1}, rot_, pos, real_type{1}, tra_); +} + +//---------------------------------------------------------------------------// +/*! + * Transform from parent to daughter. + * + * Subtract the translation, then apply the inverse of the rotation matrix (its + * transpose). + */ +CELER_FUNCTION Real3 Transformation::transform_down(Real3 const& pos) const +{ + return gemv(matrix::transpose, rot_, pos - tra_); +} + +//---------------------------------------------------------------------------// +/*! + * Rotate from daughter to parent. + */ +CELER_FUNCTION Real3 Transformation::rotate_up(Real3 const& d) const +{ + return gemv(rot_, d); +} + +//---------------------------------------------------------------------------// +/*! + * Rotate from parent to daughter. + */ +CELER_FUNCTION Real3 Transformation::rotate_down(Real3 const& d) const +{ + return gemv(matrix::transpose, rot_, d); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/transform/Translation.hh b/src/orange/transform/Translation.hh new file mode 100644 index 0000000000..7632811e4d --- /dev/null +++ b/src/orange/transform/Translation.hh @@ -0,0 +1,125 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021-2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file orange/Translator.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/cont/Span.hh" +#include "corecel/math/Algorithms.hh" +#include "corecel/math/ArrayOperators.hh" +#include "orange/OrangeTypes.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Apply a translation (no rotation). + * + * The input translation is specified by translating the daughter's coordinate + * into the parent. + */ +class Translation +{ + public: + //@{ + //! \name Type aliases + using StorageSpan = Span; + //@} + + //! Transform type identifier + static CELER_CONSTEXPR_FUNCTION TransformType transform_type() + { + return TransformType::translation; + } + + public: + //! Construct with an identity translation + explicit CELER_FUNCTION Translation() : tra_{0, 0, 0} {} + + //! Construct with the translation vector + explicit CELER_FUNCTION Translation(Real3 const& trans) : tra_(trans) {} + + // Construct inline from storage + explicit inline CELER_FUNCTION Translation(StorageSpan s); + + //// ACCESSORS //// + + //! Translation vector + CELER_FORCEINLINE_FUNCTION Real3 const& translation() const + { + return tra_; + } + + //! Get a view to the data for type-deleted storage + CELER_FUNCTION StorageSpan data() const { return {&tra_[0], 3}; } + + //// CALCULATION //// + + // Transform from daughter to parent + inline CELER_FUNCTION Real3 transform_up(Real3 const& pos) const; + + // Transform from parent to daughter + inline CELER_FUNCTION Real3 transform_down(Real3 const& parent_pos) const; + + // Rotate from daughter to parent (identity) + inline CELER_FUNCTION Real3 const& rotate_up(Real3 const& d) const; + + //! Rotate from parent to daughter (identity) + inline CELER_FUNCTION Real3 const& rotate_down(Real3 const& d) const; + + private: + Real3 tra_; +}; + +//---------------------------------------------------------------------------// +/*! + * Construct inline from storage. + */ +CELER_FUNCTION Translation::Translation(StorageSpan s) : tra_{s[0], s[1], s[2]} +{ +} + +//---------------------------------------------------------------------------// +/*! + * Transform from daughter to parent. + */ +CELER_FORCEINLINE_FUNCTION Real3 Translation::transform_up(Real3 const& pos) const +{ + return pos + tra_; +} + +//---------------------------------------------------------------------------// +/*! + * Transform from parent to daughter. + */ +CELER_FORCEINLINE_FUNCTION Real3 +Translation::transform_down(Real3 const& parent_pos) const +{ + return parent_pos - tra_; +} + +//---------------------------------------------------------------------------// +/*! + * Rotate from daughter to parent (identity). + */ +CELER_FORCEINLINE_FUNCTION Real3 const& +Translation::rotate_up(Real3 const& d) const +{ + return d; +} + +//---------------------------------------------------------------------------// +/*! + * Rotate from parent to daughter (identity). + */ +CELER_FORCEINLINE_FUNCTION Real3 const& +Translation::rotate_down(Real3 const& d) const +{ + return d; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3c0b7e608f..9a91bf6a60 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -203,13 +203,18 @@ celeritas_add_test(orange/BoundingBox.test.cc) celeritas_add_test(orange/BoundingBoxUtils.test.cc) celeritas_add_test(orange/MatrixUtils.test.cc) celeritas_add_test(orange/Orange.test.cc) -celeritas_add_test(orange/Translator.test.cc) # Base detail celeritas_add_test(orange/detail/BIHBuilder.test.cc) celeritas_add_test(orange/detail/BIHUtils.test.cc) celeritas_add_test(orange/detail/UniverseIndexer.test.cc) +#-------------------------------------# +# Transforms +set(CELERITASTEST_PREFIX orange/transform) +celeritas_add_test(orange/transform/Transformation.test.cc) +celeritas_add_test(orange/transform/Translation.test.cc) + #-------------------------------------# # Surfaces set(CELERITASTEST_PREFIX orange/surf) diff --git a/test/orange/MatrixUtils.test.cc b/test/orange/MatrixUtils.test.cc index 3347e2dd31..3022e3b610 100644 --- a/test/orange/MatrixUtils.test.cc +++ b/test/orange/MatrixUtils.test.cc @@ -50,6 +50,21 @@ TEST_F(MatrixUtilsTest, gemv) EXPECT_EQ((2 * -1) + 10 * (6 * 1 + 7 * 2 + 8 * 3), result[0]); EXPECT_EQ((2 * -2) + 10 * (6 * 0 + 7 * 1 + 8 * 2), result[1]); EXPECT_EQ((2 * -3) + 10 * (6 * 4 + 7 * 0 + 8 * 1), result[2]); + + result = gemv(matrix::transpose, 10, mat, {6, 7, 8}, 2, {-1, -2, -3}); + EXPECT_EQ((2 * -1) + 10 * (6 * 1 + 7 * 0 + 8 * 4), result[0]); + EXPECT_EQ((2 * -2) + 10 * (6 * 2 + 7 * 1 + 8 * 0), result[1]); + EXPECT_EQ((2 * -3) + 10 * (6 * 3 + 7 * 2 + 8 * 1), result[2]); + + result = gemv(mat, {6, 7, 8}); + EXPECT_EQ((6 * 1 + 7 * 2 + 8 * 3), result[0]); + EXPECT_EQ((6 * 0 + 7 * 1 + 8 * 2), result[1]); + EXPECT_EQ((6 * 4 + 7 * 0 + 8 * 1), result[2]); + + result = gemv(matrix::transpose, mat, {6, 7, 8}); + EXPECT_EQ((6 * 1 + 7 * 0 + 8 * 4), result[0]); + EXPECT_EQ((6 * 2 + 7 * 1 + 8 * 0), result[1]); + EXPECT_EQ((6 * 3 + 7 * 2 + 8 * 1), result[2]); } //---------------------------------------------------------------------------// diff --git a/test/orange/Translator.test.cc b/test/orange/Translator.test.cc deleted file mode 100644 index 7a5aae5762..0000000000 --- a/test/orange/Translator.test.cc +++ /dev/null @@ -1,39 +0,0 @@ -//----------------------------------*-C++-*----------------------------------// -// Copyright 2021-2023 UT-Battelle, LLC, and other Celeritas developers. -// See the top-level COPYRIGHT file for details. -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file orange/Translator.test.cc -//---------------------------------------------------------------------------// -#include "orange/Translator.hh" - -#include "celeritas_test.hh" - -namespace celeritas -{ -namespace test -{ -//---------------------------------------------------------------------------// -class TranslatorTest : public Test -{ - protected: - Real3 translation_{1, 2, 3}; -}; - -TEST_F(TranslatorTest, down) -{ - TranslatorDown translate(translation_); - - EXPECT_VEC_SOFT_EQ((Real3{.1, .2, .3}), translate(Real3{1.1, 2.2, 3.3})); -} - -TEST_F(TranslatorTest, up) -{ - TranslatorUp translate(translation_); - - EXPECT_VEC_SOFT_EQ((Real3{1.1, 2.2, 3.3}), translate(Real3{.1, .2, .3})); -} - -//---------------------------------------------------------------------------// -} // namespace test -} // namespace celeritas diff --git a/test/orange/transform/Transformation.test.cc b/test/orange/transform/Transformation.test.cc new file mode 100644 index 0000000000..bdddd95a5e --- /dev/null +++ b/test/orange/transform/Transformation.test.cc @@ -0,0 +1,99 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file orange/transform/Transformation.test.cc +//---------------------------------------------------------------------------// +#include "orange/transform/Transformation.hh" + +#include "orange/MatrixUtils.hh" +#include "orange/transform/Translation.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// + +class TransformationTest : public ::celeritas::test::Test +{ + protected: + template + std::vector flattened(SquareMatrix const& inp) const + { + std::vector result; + + for (auto& row : inp) + { + result.insert(result.end(), row.begin(), row.end()); + } + return result; + } +}; + +TEST_F(TransformationTest, construction) +{ + static double const identity[] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + { + SCOPED_TRACE("default construct"); + Transformation tr; + EXPECT_VEC_SOFT_EQ(identity, flattened(tr.rotation())); + EXPECT_VEC_SOFT_EQ((Real3{0, 0, 0}), tr.translation()); + } + { + SCOPED_TRACE("promotion"); + Transformation tr{Translation{Real3{1, 2, 3}}}; + EXPECT_VEC_SOFT_EQ(identity, flattened(tr.rotation())); + EXPECT_VEC_SOFT_EQ((Real3{1, 2, 3}), tr.translation()); + } + { + SCOPED_TRACE("serialize"); + Transformation tr{make_rotation(Axis::z, Turn{0.125}), {1, 2, 3}}; + Transformation tr2{tr.data()}; + EXPECT_EQ(tr.translation(), tr2.translation()); + EXPECT_EQ(tr.rotation(), tr2.rotation()); + } +} + +TEST_F(TransformationTest, transform) +{ + { + Transformation tr{make_rotation(Axis::z, Turn{0.25}), Real3{0, 0, 1}}; + // Daughter to parent: rotate quarter turn around Z, then add 1 to Z + EXPECT_VEC_EQ((Real3{-3, 2, 1}), tr.transform_up({2, 3, 0})); + // Parent to daughter: subtract, then rotate back + EXPECT_VEC_EQ((Real3{2, 3, 0}), tr.transform_down({-3, 2, 1})); + } + { + Transformation tr{ + make_rotation( + Axis::x, + native_value_to(std::acos(-0.5)), + make_rotation(Axis::y, native_value_to(std::acos(0.2)))), + Real3{1.1, -0.5, 3.2}}; + + auto const daughter = Real3{-3.4, 2.1, 0.4}; + auto const parent = tr.transform_up(daughter); + EXPECT_VEC_NEAR( + (Real3{0.81191836, -4.5042777, 3.31300032}), parent, 1e-6); + EXPECT_VEC_SOFT_EQ(daughter, tr.transform_down(parent)); + } +} + +TEST_F(TransformationTest, rotate) +{ + { + Transformation tr{make_rotation(Axis::z, Turn{0.25}), Real3{0, 0, 1}}; + // Daughter to parent: rotate quarter turn around Z + EXPECT_VEC_EQ((Real3{0, 1, 0}), tr.rotate_up({1, 0, 0})); + // Parent to daughter: rotate back + EXPECT_VEC_EQ((Real3{1, 0, 0}), tr.rotate_down({0, 1, 0})); + } +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/orange/transform/Translation.test.cc b/test/orange/transform/Translation.test.cc new file mode 100644 index 0000000000..6e9d2cbc0f --- /dev/null +++ b/test/orange/transform/Translation.test.cc @@ -0,0 +1,52 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021-2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file orange/transform/Translation.test.cc +//---------------------------------------------------------------------------// +#include "orange/transform/Translation.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +class TranslatorTest : public Test +{ +}; + +TEST_F(TranslatorTest, translation) +{ + Translation tr(Real3{1, 2, 3}); + + EXPECT_VEC_SOFT_EQ((Real3{.1, .2, .3}), + tr.transform_down(Real3{1.1, 2.2, 3.3})); + EXPECT_VEC_SOFT_EQ((Real3{1.1, 2.2, 3.3}), + tr.transform_up(Real3{.1, .2, .3})); +} + +TEST_F(TranslatorTest, rotation) +{ + Translation tr(Real3{1, 2, 3}); + + Real3 dir{0, 0, 1}; + EXPECT_VEC_SOFT_EQ(dir, tr.rotate_down(dir)); + EXPECT_VEC_SOFT_EQ(dir, tr.rotate_up(dir)); +} + +TEST_F(TranslatorTest, serialization) +{ + Translation tr{Real3{3, 2, 1}}; + + EXPECT_VEC_EQ((Real3{3, 2, 1}), tr.data()); + + Translation tr2(tr.data()); + EXPECT_VEC_EQ((Real3{3, 2, 1}), tr2.translation()); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas