Add transform simplifier (celeritas-project#1131)

* Add trace utility for matrices
* Add arbitrary-direction rotation matrix helper
* Add transform simplifier
* Fix and document single precision issues
* Add deduplication test
* Deduplicate transforms on insertion
* Address review feedback

src/orange/CMakeLists.txt

@@ -58,6 +58,7 @@ list(APPEND SOURCES
   transform/SignedPermutation.cc
   transform/TransformHasher.cc
   transform/TransformIO.cc
+  transform/TransformSimplifier.cc
   transform/Transformation.cc
   transform/VariantTransform.cc
 )

src/orange/MatrixUtils.cc

@@ -23,11 +23,26 @@ namespace celeritas
 template<class T>
 T determinant(SquareMatrix<T, 3> const& mat)
 {
-    return mat[0][0] * mat[1][1] * mat[2][2] + mat[1][0] * mat[2][1] * mat[0][2]
+    // clang-format off
+    return   mat[0][0] * mat[1][1] * mat[2][2]
+           + mat[1][0] * mat[2][1] * mat[0][2]
            + mat[2][0] * mat[0][1] * mat[1][2]
            - mat[2][0] * mat[1][1] * mat[0][2]
            - mat[1][0] * mat[0][1] * mat[2][2]
            - mat[0][0] * mat[2][1] * mat[1][2];
+    // clang-format on
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate the trace of a 3x3 matrix.
+ *
+ * The trace is just the sum of the diagonal elements.
+ */
+template<class T>
+T trace(SquareMatrix<T, 3> const& mat)
+{
+    return mat[0][0] + mat[1][1] + mat[2][2];
 }
 
 //---------------------------------------------------------------------------//
@@ -136,6 +151,48 @@ void orthonormalize(SquareMatrix<T, N>* mat)
     CELER_ENSURE(soft_equal(std::fabs(determinant(*mat)), T{1}));
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Create a C-ordered rotation matrix from an arbitrary rotation.
+ *
+ * This is equation (38) in "Rotation Matrices in Two, Three, and Many
+ * Dimensions",  Physics 116A, UC Santa Cruz,
+ * http://scipp.ucsc.edu/~haber/ph116A/.
+ *
+ * \param ax Axis of rotation (unit vector)
+ * \param theta Rotation
+ */
+Mat3 make_rotation(Real3 const& ax, Turn theta)
+{
+    CELER_EXPECT(is_soft_unit_vector(ax));
+    CELER_EXPECT(theta >= Turn{0} && theta <= Turn{real_type(0.5)});
+
+    // Axis/direction enumeration
+    enum
+    {
+        X = 0,
+        Y = 1,
+        Z = 2
+    };
+
+    // Calculate sin and cosine with less precision loss using "turn" value
+    real_type cost;
+    real_type sint;
+    sincos(theta, &sint, &cost);
+
+    Mat3 r{Real3{cost + ipow<2>(ax[X]) * (1 - cost),
+                 ax[X] * ax[Y] * (1 - cost) - ax[Z] * sint,
+                 ax[X] * ax[Z] * (1 - cost) + ax[Y] * sint},
+           Real3{ax[X] * ax[Y] * (1 - cost) + ax[Z] * sint,
+                 cost + ipow<2>(ax[Y]) * (1 - cost),
+                 ax[Y] * ax[Z] * (1 - cost) - ax[X] * sint},
+           Real3{ax[X] * ax[Z] * (1 - cost) - ax[Y] * sint,
+                 ax[Y] * ax[Z] * (1 - cost) + ax[X] * sint,
+                 cost + ipow<2>(ax[Z]) * (1 - cost)}};
+    CELER_ENSURE(soft_equal(std::fabs(determinant(r)), real_type{1}));
+    return r;
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Create a C-ordered rotation matrix.
@@ -190,6 +247,10 @@ template int determinant(SquareMatrix<int, 3> const&);
 template float determinant(SquareMatrix<float, 3> const&);
 template double determinant(SquareMatrix<double, 3> const&);
 
+template int trace(SquareMatrix<int, 3> const&);
+template float trace(SquareMatrix<float, 3> const&);
+template double trace(SquareMatrix<double, 3> const&);
+
 template void orthonormalize(SquareMatrix<float, 3>*);
 template void orthonormalize(SquareMatrix<double, 3>*);
 

src/orange/MatrixUtils.hh

@@ -74,6 +74,10 @@ gemv(matrix::TransposePolicy, SquareMatrix<T, N> const& a, Array<T, N> const& x)
 template<class T>
 T determinant(SquareMatrix<T, 3> const& mat);
 
+// Calculate the trace of a matrix
+template<class T>
+T trace(SquareMatrix<T, 3> const& mat);
+
 // Perform a matrix-matrix multiply
 template<class T, size_type N>
 SquareMatrix<T, N>
@@ -89,7 +93,10 @@ SquareMatrix<T, N> gemm(matrix::TransposePolicy,
 template<class T, size_type N>
 void orthonormalize(SquareMatrix<T, N>* mat);
 
-// Create a C-ordered rotation matrix
+// Create a C-ordered rotation matrix about an arbitrary axis
+SquareMatrixReal3 make_rotation(Real3 const& ax, Turn rev);
+
+// Create a C-ordered rotation matrix about a cartesian axis
 SquareMatrixReal3 make_rotation(Axis ax, Turn rev);
 
 // Apply a rotation to an existing C-ordered rotation matrix

src/orange/orangeinp/detail/CsgUnitBuilder.cc

@@ -10,6 +10,7 @@
 #include "corecel/io/Logger.hh"
 #include "corecel/io/StreamableVariant.hh"
 #include "orange/transform/TransformIO.hh"
+#include "orange/transform/TransformSimplifier.hh"
 
 namespace celeritas
 {
@@ -49,6 +50,16 @@ BoundingZone const& CsgUnitBuilder::bounds(NodeId nid) const
     return iter->second.bounds;
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Insert transform with simplification and deduplication.
+ */
+TransformId CsgUnitBuilder::insert_transform(VariantTransform const& vt)
+{
+    auto simplified = std::visit(TransformSimplifier(tol_), vt);
+    return this->insert_transform_(std::move(simplified));
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Set a bounding zone and transform for a node.

src/orange/orangeinp/detail/CsgUnitBuilder.hh

@@ -79,7 +79,7 @@ class CsgUnitBuilder
     inline NodeInsertion insert_csg(Args&&... args);
 
     // Insert a transform
-    inline TransformId insert_transform(VariantTransform&& vt);
+    TransformId insert_transform(VariantTransform const& vt);
 
     // Insert node metadata
     inline void insert_md(NodeId node, Metadata&& md);
@@ -173,15 +173,6 @@ auto CsgUnitBuilder::insert_csg(Args&&... args) -> NodeInsertion
     return result;
 }
 
-//---------------------------------------------------------------------------//
-/*!
- * Insert transform with deduplication.
- */
-TransformId CsgUnitBuilder::insert_transform(VariantTransform&& vt)
-{
-    return this->insert_transform_(std::move(vt));
-}
-
 //---------------------------------------------------------------------------//
 /*!
  * Insert node metadata.

src/orange/surf/SoftSurfaceEqual.hh

@@ -38,7 +38,7 @@ class SoftSurfaceEqual
 {
   public:
     // Construct with tolerance
-    inline SoftSurfaceEqual(Tolerance<> const& tol);
+    explicit inline SoftSurfaceEqual(Tolerance<> const& tol);
 
     //! Construct with relative tolerance only
     explicit SoftSurfaceEqual(real_type rel)

src/orange/surf/SurfaceSimplifier.hh

@@ -29,6 +29,9 @@ namespace celeritas
  * was applied.
  *
  * The embedded sense may be flipped as part of the simplifications.
+ *
+ * \todo Use a \c Tolerance object instead of a single tolerance, and compare
+ * implementations with \c SoftSurfaceEqual for consistency.
  */
 class SurfaceSimplifier
 {

src/orange/transform/TransformSimplifier.cc

@@ -0,0 +1,47 @@
+//----------------------------------*-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/TransformSimplifier.cc
+//---------------------------------------------------------------------------//
+#include "TransformSimplifier.hh"
+
+#include "corecel/math/ArrayUtils.hh"
+#include "orange/MatrixUtils.hh"
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Translation may simplify to no transformation.
+ */
+VariantTransform TransformSimplifier::operator()(Translation const& t)
+{
+    if (norm(t.translation()) <= eps_)
+    {
+        // Distance to origin is less than tolerance
+        return NoTransformation{};
+    }
+    return t;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Simplify, possibly to translation or no transform.
+ *
+ * See the derivation in the class documentation.
+ */
+VariantTransform TransformSimplifier::operator()(Transformation const& t)
+{
+    real_type tr = trace(t.rotation());
+    if (tr >= 3 - ipow<2>(eps_))
+    {
+        // Rotation results in no more then epsilon movement
+        return (*this)(Translation{t.translation()});
+    }
+    return t;
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

src/orange/transform/TransformSimplifier.hh

@@ -0,0 +1,77 @@
+//----------------------------------*-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/TransformSimplifier.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "corecel/Macros.hh"
+#include "orange/OrangeTypes.hh"
+
+#include "VariantTransform.hh"
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Return a simplified version of a transformation.
+ *
+ * Like surface simplification, we want to consider whether two different
+ * transformations will result in a distance change of \f$\epsilon\f$ for a
+ * point that's at the length scale from the origin. Setting the length scale
+ * to unity (the default), we use the relative tolerance.
+ *
+ * A *translation* can be deleted if its magnitude is less than epsilon.
+ *
+ * For a translation, we use the fact that the trace (sum of diagonal elements)
+ * of any proper (non-reflecting) rotation matrix has an angle of rotation \f[
+   \mathrm{Tr}[R] = 2 \cos \theta + 1
+ * \f] about an axis of rotation. Applying the rotation to a point
+ * at a unit distance will yield an iscoceles triangle with sides 1 and inner
+ * angle \f$\theta\f$. For the displacement to be no more than \f$\epsilon\f$
+ * then the angle must be \f[
+  \sin \theta/2 \le \epsilon/2
+  \f]
+ * which with some manipulation means that a "soft zero" rotation has a trace
+ * \f[
+ \mathrm{Tr}[R] \ge 3 - \epsilon^2 \,.
+ \f]
+ *
+ * Note that this means no rotational simplifications may be performed when the
+ * geometry tolerance is less than the square root of machine precision.
+ */
+class TransformSimplifier
+{
+  public:
+    // Construct with tolerance
+    explicit inline TransformSimplifier(Tolerance<> const& tol);
+
+    //! No simplification can be applied to a null transformation
+    VariantTransform operator()(NoTransformation const& nt) { return nt; }
+
+    // Translation may simplify to no transformation
+    VariantTransform operator()(Translation const& nt);
+
+    // Simplify, possibly to translation or no transform
+    VariantTransform operator()(Transformation const& nt);
+
+  private:
+    real_type eps_;
+};
+
+//---------------------------------------------------------------------------//
+// INLINE DEFINITIONS
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with tolerance.
+ */
+TransformSimplifier::TransformSimplifier(Tolerance<> const& tol)
+    : eps_{tol.rel}
+{
+    CELER_EXPECT(tol);
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

test/corecel/math/Algorithms.test.cc

@@ -363,27 +363,38 @@ TEST(MathTest, sincospi)
     EXPECT_DOUBLE_EQ(std::sin(m_pi * 0.1), sinpi(0.1));
     EXPECT_DOUBLE_EQ(std::cos(m_pi * 0.1), cospi(0.1));
 
-    double s{0}, c{0};
-    sincospi(0.123, &s, &c);
-    EXPECT_DOUBLE_EQ(std::sin(m_pi * 0.123), s);
-    EXPECT_DOUBLE_EQ(std::cos(m_pi * 0.123), c);
-
-    // Test special cases
-    sincospi(0, &s, &c);
-    EXPECT_EQ(double(0.0), s);
-    EXPECT_EQ(double(1.0), c);
-
-    sincospi(0.5, &s, &c);
-    EXPECT_EQ(double(1.0), s);
-    EXPECT_EQ(double(0.0), c);
-
-    sincospi(1.0, &s, &c);
-    EXPECT_EQ(double(0.0), s);
-    EXPECT_EQ(double(-1.0), c);
-
-    sincospi(1.5, &s, &c);
-    EXPECT_EQ(double(-1.0), s);
-    EXPECT_EQ(double(0.0), c);
+    {
+        double s{0}, c{0};
+        sincospi(0.123, &s, &c);
+        EXPECT_DOUBLE_EQ(std::sin(m_pi * 0.123), s);
+        EXPECT_DOUBLE_EQ(std::cos(m_pi * 0.123), c);
+
+        // Test special cases
+        sincospi(0, &s, &c);
+        EXPECT_EQ(double(0.0), s);
+        EXPECT_EQ(double(1.0), c);
+
+        sincospi(0.5, &s, &c);
+        EXPECT_EQ(double(1.0), s);
+        EXPECT_EQ(double(0.0), c);
+
+        sincospi(1.0, &s, &c);
+        EXPECT_EQ(double(0.0), s);
+        EXPECT_EQ(double(-1.0), c);
+
+        sincospi(1.5, &s, &c);
+        EXPECT_EQ(double(-1.0), s);
+        EXPECT_EQ(double(0.0), c);
+    }
+    {
+        // This is about the threshold where sin(x) == 1.0f
+        float inp{0.000233115f};
+        float s{0}, c{0};
+        sincospi(inp, &s, &c);
+        EXPECT_FLOAT_EQ(1.0f, c);
+        EXPECT_FLOAT_EQ(std::sin(static_cast<float>(m_pi * inp)), s);
+        EXPECT_FLOAT_EQ(std::cos(static_cast<float>(m_pi * inp)), c);
+    }
 }
 
 //---------------------------------------------------------------------------//

test/orange/CMakeLists.txt

@@ -12,9 +12,9 @@ celeritas_add_library(testcel_orange
 )
 celeritas_target_link_libraries(testcel_orange
   PUBLIC
-    Celeritas::orange Celeritas::testcel_harness
+    Celeritas::orange Celeritas::testcel_harness Celeritas::testcel_core
   PRIVATE
-    Celeritas::testcel_geocel Celeritas::testcel_core Celeritas::orange ${nlohmann_json_LIBRARIES}
+    Celeritas::testcel_geocel Celeritas::orange ${nlohmann_json_LIBRARIES}
 )
 
 #-----------------------------------------------------------------------------#
@@ -69,6 +69,7 @@ celeritas_add_test(orangeinp/detail/TransformInserter.test.cc)
 # Transforms
 set(CELERITASTEST_PREFIX orange/transform)
 celeritas_add_test(transform/SignedPermutation.test.cc)
+celeritas_add_test(transform/TransformSimplifier.test.cc)
 celeritas_add_test(transform/Transformation.test.cc)
 celeritas_add_test(transform/Translation.test.cc)
 celeritas_add_test(transform/VariantTransform.test.cc)