Skip to content

Commit

Permalink
Add specialized "permutation" transform (#1077)
Browse files Browse the repository at this point in the history
* WIP: signed permutation

* Simplify rotation matrix construction

* Define integer quarter turn helper class

* Add signed permutation class

* Rename pos/neg to remove ambiguity

* Use explicit characters and add determinant error checking

* Add documentation
  • Loading branch information
sethrj committed Jan 20, 2024
1 parent b1e52f2 commit 33eea0f
Show file tree
Hide file tree
Showing 11 changed files with 566 additions and 29 deletions.
2 changes: 2 additions & 0 deletions src/corecel/cont/EnumArray.hh
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@ namespace celeritas
*
* The enum *must* be a zero-indexed contiguous enumeration with a \c size_
* enumeration as its last value.
*
* \todo The template parameters are reversed!!!
*/
template<class E, class T>
struct EnumArray
Expand Down
33 changes: 33 additions & 0 deletions src/corecel/math/Turn.hh
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
//---------------------------------------------------------------------------//
#pragma once

#include <cmath>

#include "corecel/Types.hh"

#include "Algorithms.hh"
Expand All @@ -23,10 +25,23 @@ struct TwoPi
static char const* label() { return "tr"; }
};

//---------------------------------------------------------------------------//
//! Unit for pi/2 radians
struct HalfPi
{
static real_type value() { return static_cast<real_type>(m_pi / 2); }
//! Text label for output
static char const* label() { return "qtr"; }
};

//---------------------------------------------------------------------------//
//! Quantity denoting a full turn
using Turn = Quantity<TwoPi, real_type>;

//---------------------------------------------------------------------------//
//! Quantity for an integer number of turns for axis swapping
using QuarterTurn = Quantity<HalfPi, int>;

//---------------------------------------------------------------------------//
//!@{
//! Special overrides for math functions for more precise arithmetic
Expand All @@ -44,6 +59,24 @@ CELER_FORCEINLINE_FUNCTION void sincos(Turn r, real_type* sinv, real_type* cosv)
{
return sincospi(r.value() * 2, sinv, cosv);
}

CELER_FORCEINLINE_FUNCTION int cos(QuarterTurn r)
{
constexpr int cosval[] = {1, 0, -1, 0};
return cosval[std::abs(r.value()) % 4];
}

CELER_FORCEINLINE_FUNCTION int sin(QuarterTurn r)
{
// Define in terms of the symmetric "cos"
return cos(QuarterTurn{r.value() - 1});
}

CELER_FORCEINLINE_FUNCTION void sincos(QuarterTurn r, int* sinv, int* cosv)
{
*sinv = sin(r);
*cosv = cos(r);
}
//!@}

//---------------------------------------------------------------------------//
Expand Down
1 change: 1 addition & 0 deletions src/orange/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ list(APPEND SOURCES
surf/detail/SurfaceClipperImpl.cc
surf/detail/SurfaceTranslator.cc
surf/detail/SurfaceTransformer.cc
transform/SignedPermutation.cc
transform/Transformation.cc
transform/VariantTransform.cc
)
Expand Down
34 changes: 7 additions & 27 deletions src/orange/MatrixUtils.cc
Original file line number Diff line number Diff line change
Expand Up @@ -144,44 +144,24 @@ Mat3 make_rotation(Axis ax, Turn theta)
{
CELER_EXPECT(ax < Axis::size_);

// Calculate sin and cosine with less precision loss
// Calculate sin and cosine with less precision loss using "turn" value
real_type cost;
real_type sint;
sincos(theta, &sint, &cost);

// Avoid signed zeros with this implementation
real_type msint = real_type(0) - sint;

// Fill result with zeros
Mat3 r;
r.fill(Real3{0, 0, 0});

// {i, i} gets 1
r[to_int(ax)][to_int(ax)] = 1;

switch (ax)
{
case Axis::x:
r[1][1] = cost;
r[1][2] = msint;
r[2][1] = sint;
r[2][2] = cost;
break;
case Axis::y:
r[0][0] = cost;
r[0][2] = sint;
r[2][0] = msint;
r[2][2] = cost;
break;
case Axis::z:
r[0][0] = cost;
r[0][1] = msint;
r[1][0] = sint;
r[1][1] = cost;
break;
default:
CELER_ASSERT_UNREACHABLE();
}
int uax = (to_int(ax) + 1) % 3;
int vax = (to_int(ax) + 2) % 3;
r[uax][uax] = cost;
r[uax][vax] = negate(sint); // avoid signed zeros
r[vax][uax] = sint;
r[vax][vax] = cost;
return r;
}

Expand Down
140 changes: 140 additions & 0 deletions src/orange/transform/SignedPermutation.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
//----------------------------------*-C++-*----------------------------------//
// Copyright 2024 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/SignedPermutation.cc
//---------------------------------------------------------------------------//
#include "SignedPermutation.hh"

#include "orange/MatrixUtils.hh"

namespace celeritas
{
//---------------------------------------------------------------------------//
/*!
* Construct with an identity permutation.
*/
SignedPermutation::SignedPermutation()
: SignedPermutation{
SignedAxes{{{'+', Axis::x}, {'+', Axis::y}, {'+', Axis::z}}}}
{
}

//---------------------------------------------------------------------------//
/*!
* Construct with a permutation vector.
*/
SignedPermutation::SignedPermutation(SignedAxes permutation) : compressed_{0}
{
EnumArray<Axis, bool> encountered_ax{false, false, false};

SquareMatrix<int, 3> explicit_mat;

for (auto ax : {Axis::z, Axis::y, Axis::x})
{
auto new_ax_sign = permutation[ax];
CELER_VALIDATE(new_ax_sign.first == '+' || new_ax_sign.first == '-',
<< "invalid permutation sign '" << new_ax_sign.first
<< "'");
CELER_VALIDATE(new_ax_sign.second < Axis::size_,
<< "invalid permutation axis");
CELER_VALIDATE(!encountered_ax[new_ax_sign.second],
<< "duplicate axis " << to_char(new_ax_sign.second)
<< " in permutation matrix input");
encountered_ax[new_ax_sign.second] = true;

// Push back "flip bit"
compressed_ <<= 1;
compressed_ |= (new_ax_sign.first == '-' ? 0b1 : 0b0);
// Push back "new axis"
compressed_ <<= 2;
compressed_ |= to_int(new_ax_sign.second);

// Explicitly construct row in error checking matrix
for (auto oax : {Axis::x, Axis::y, Axis::z})
{
explicit_mat[to_int(ax)][to_int(oax)] = [oax, new_ax_sign] {
if (oax != new_ax_sign.second)
return 0;
if (new_ax_sign.first == '-')
return -1;
return 1;
}();
}
}
int det = determinant(explicit_mat);
CELER_VALIDATE(
det == +1,
<< "invalid rotation matrix: determinant should be +1 but is " << det);
}

//---------------------------------------------------------------------------//
/*!
* Reconstruct the permutation.
*/
auto SignedPermutation::permutation() const -> SignedAxes
{
SignedAxes result;

UIntT temp = compressed_;
for (auto ax : range(Axis::size_))
{
// Copy "new axis"
result[ax].second = to_axis(temp & 0b11);
temp >>= 2;
// Push back "flip bit"
result[ax].first = temp & 0b1 ? '-' : '+';
temp >>= 1;
}
return result;
}

//---------------------------------------------------------------------------//
/*!
* Get a view to the data for type-deleted storage.
*/
auto SignedPermutation::data() const -> DataArray
{
static_assert(
max_value() == static_cast<UIntT>(static_cast<real_type>(max_value())),
"cannot round-trip real type and short int");

return {static_cast<real_type>(compressed_)};
}

//---------------------------------------------------------------------------//
/*!
* Make a permutation by rotating about the given axis.
*/
SignedPermutation make_permutation(Axis ax, QuarterTurn theta)
{
CELER_EXPECT(ax < Axis::size_);

auto to_sign = [](int v) { return v < 0 ? '-' : '+'; };

int const cost = cos(theta);
int const sint = sin(theta);
Axis const uax = to_axis((to_int(ax) + 1) % 3);
Axis const vax = to_axis((to_int(ax) + 2) % 3);

SignedPermutation::SignedAxes r;

// Axis of rotation is unchanged
r[ax] = {to_sign(1), ax};
if (cost != 0)
{
r[uax] = {to_sign(cost), uax};
r[vax] = {to_sign(cost), vax};
}
else
{
r[uax] = {to_sign(-sint), vax};
r[vax] = {to_sign(sint), uax};
}

return SignedPermutation{r};
}

//---------------------------------------------------------------------------//
} // namespace celeritas

0 comments on commit 33eea0f

Please sign in to comment.