Fix single-precision soft equivalence for plane (celeritas-project#1139)

* Uninline local surface inserter

* Fix relative instead of absolute tolerance

* Add test demonstrating failure for single-precision LSI

* Add surface soft equal documentation for plane

* Add additional plane test

* Fix plane comparison by adding extra machine precision

* Fix pedantic semicolons

* fixup! Fix relative instead of absolute tolerance

doc/CMakeLists.txt

@@ -58,6 +58,7 @@ set(DOXYGEN_QT_AUTOBRIEF YES)
 set(DOXYGEN_DISTRIBUTE_GROUP_DOC YES)
 set(DOXYGEN_USE_MDFILE_AS_MAINPAGE "${PROJECT_SOURCE_DIR}/README.md")
 set(DOXYGEN_FORMULA_MACROFILE "${CMAKE_CURRENT_SOURCE_DIR}/_static/macros.tex")
+set(DOXYGEN_IMAGE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/_static")
 set(DOXYGEN_EXTRA_PACKAGES "amsmath")
 set(DOXYGEN_PREDEFINED "__DOXYGEN__=${_doxy_version_hex}")
 set(DOXYGEN_USE_MATHJAX YES)

src/orange/orangeinp/ConvexSurfaceBuilder.cc

@@ -188,7 +188,7 @@ void visit(ConvexSurfaceBuilder& csb, Sense sense, VariantSurface const& surf)
 //---------------------------------------------------------------------------//
 //! \cond
 #define CSB_INSTANTIATE(SURF) \
-    template void ConvexSurfaceBuilder::operator()(Sense, SURF const&);
+    template void ConvexSurfaceBuilder::operator()(Sense, SURF const&)
 CSB_INSTANTIATE(PlaneAligned<Axis::x>);
 CSB_INSTANTIATE(PlaneAligned<Axis::y>);
 CSB_INSTANTIATE(PlaneAligned<Axis::z>);

src/orange/orangeinp/detail/LocalSurfaceInserter.cc

@@ -25,6 +25,57 @@ LocalSurfaceInserter::LocalSurfaceInserter(VecSurface* v,
     CELER_EXPECT(surfaces_->empty());
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Construct a surface with deduplication.
+ */
+template<class S>
+LocalSurfaceId LocalSurfaceInserter::operator()(S const& source)
+{
+    VecSurface& all_surf = *surfaces_;
+
+    // Test for soft equality against all existing surfaces
+    auto is_soft_equal = [this, &source](VariantSurface const& vtarget) {
+        if (S const* target = std::get_if<S>(&vtarget))
+        {
+            return soft_surface_equal_(source, *target);
+        }
+        return false;
+    };
+
+    // TODO: instead of linear search (making overall unit insertion
+    // quadratic!) accelerate by mapping surfaces into a hash and comparing all
+    // neighboring hash cells
+    auto iter = std::find_if(all_surf.begin(), all_surf.end(), is_soft_equal);
+    if (iter == all_surf.end())
+    {
+        // Surface is completely unique
+        LocalSurfaceId result(all_surf.size());
+        all_surf.emplace_back(std::in_place_type<S>, source);
+        return result;
+    }
+
+    // Surface is soft equivalent to an existing surface at this index
+    LocalSurfaceId target_id{static_cast<size_type>(iter - all_surf.begin())};
+
+    CELER_ASSUME(std::holds_alternative<S>(*iter));
+    if (exact_surface_equal_(source, std::get<S>(*iter)))
+    {
+        // Surfaces are *exactly* equal: don't insert
+        CELER_ENSURE(target_id < all_surf.size());
+        return target_id;
+    }
+
+    // Surface is a little bit different, so we still need to insert it
+    // to chain duplicates
+    LocalSurfaceId source_id(all_surf.size());
+    all_surf.emplace_back(std::in_place_type<S>, source);
+
+    // Store the equivalency relationship and potentially chain equivalent
+    // surfaces
+    return this->merge_impl(source_id, target_id);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Look for duplicate surfaces and store equivalency relationship.
@@ -46,6 +97,34 @@ LocalSurfaceInserter::merge_impl(LocalSurfaceId source, LocalSurfaceId target)
     return target;
 }
 
+//---------------------------------------------------------------------------//
+// EXPLICIT INSTANTIATIONS
+//---------------------------------------------------------------------------//
+//! \cond
+#define LSI_INSTANTIATE(SURF) \
+    template LocalSurfaceId LocalSurfaceInserter::operator()(SURF const&)
+
+LSI_INSTANTIATE(PlaneAligned<Axis::x>);
+LSI_INSTANTIATE(PlaneAligned<Axis::y>);
+LSI_INSTANTIATE(PlaneAligned<Axis::z>);
+LSI_INSTANTIATE(CylCentered<Axis::x>);
+LSI_INSTANTIATE(CylCentered<Axis::y>);
+LSI_INSTANTIATE(CylCentered<Axis::z>);
+LSI_INSTANTIATE(SphereCentered);
+LSI_INSTANTIATE(CylAligned<Axis::x>);
+LSI_INSTANTIATE(CylAligned<Axis::y>);
+LSI_INSTANTIATE(CylAligned<Axis::z>);
+LSI_INSTANTIATE(Plane);
+LSI_INSTANTIATE(Sphere);
+LSI_INSTANTIATE(ConeAligned<Axis::x>);
+LSI_INSTANTIATE(ConeAligned<Axis::y>);
+LSI_INSTANTIATE(ConeAligned<Axis::z>);
+LSI_INSTANTIATE(SimpleQuadric);
+LSI_INSTANTIATE(GeneralQuadric);
+
+#undef LSI_INSTANTIATE
+//! \endcond
+
 //---------------------------------------------------------------------------//
 }  // namespace detail
 }  // namespace orangeinp

src/orange/orangeinp/detail/LocalSurfaceInserter.hh

@@ -53,7 +53,7 @@ class LocalSurfaceInserter
 
     // Construct a surface with deduplication
     template<class S>
-    inline LocalSurfaceId operator()(S const& surface);
+    LocalSurfaceId operator()(S const& surface);
 
   private:
     //// TYPES ////
@@ -72,59 +72,6 @@ class LocalSurfaceInserter
     LocalSurfaceId merge_impl(LocalSurfaceId source, LocalSurfaceId target);
 };
 
-//---------------------------------------------------------------------------//
-// INLINE DEFINITIONS
-//---------------------------------------------------------------------------//
-/*!
- * Construct a surface with deduplication.
- */
-template<class S>
-LocalSurfaceId LocalSurfaceInserter::operator()(S const& source)
-{
-    VecSurface& all_surf = *surfaces_;
-
-    // Test for soft equality against all existing surfaces
-    auto is_soft_equal = [this, &source](VariantSurface const& vtarget) {
-        if (S const* target = std::get_if<S>(&vtarget))
-        {
-            return soft_surface_equal_(source, *target);
-        }
-        return false;
-    };
-
-    // TODO: instead of linear search (making overall unit insertion
-    // quadratic!) accelerate by mapping surfaces into a hash and comparing all
-    // neighboring hash cells
-    auto iter = std::find_if(all_surf.begin(), all_surf.end(), is_soft_equal);
-    if (iter == all_surf.end())
-    {
-        // Surface is completely unique
-        LocalSurfaceId result(all_surf.size());
-        all_surf.emplace_back(std::in_place_type<S>, source);
-        return result;
-    }
-
-    // Surface is soft equivalent to an existing surface at this index
-    LocalSurfaceId target_id{static_cast<size_type>(iter - all_surf.begin())};
-
-    CELER_ASSUME(std::holds_alternative<S>(*iter));
-    if (exact_surface_equal_(source, std::get<S>(*iter)))
-    {
-        // Surfaces are *exactly* equal: don't insert
-        CELER_ENSURE(target_id < all_surf.size());
-        return target_id;
-    }
-
-    // Surface is a little bit different, so we still need to insert it
-    // to chain duplicates
-    LocalSurfaceId source_id(all_surf.size());
-    all_surf.emplace_back(std::in_place_type<S>, source);
-
-    // Store the equivalency relationship and potentially chain equivalent
-    // surfaces
-    return this->merge_impl(source_id, target_id);
-}
-
 //---------------------------------------------------------------------------//
 }  // namespace detail
 }  // namespace orangeinp

src/orange/surf/RecursiveSimplifier.hh

@@ -22,8 +22,7 @@ namespace celeritas
 /*!
  * Recursively simplify, then call the given function on the final surface.
  *
- * The tolerance for this class should be "absolute" tolerance, i.e., the
- * relative tolerance for an O(1) sized object.
+ * The tolerance for this class should be relative.
  *
  * Example:
  * \code
@@ -49,7 +48,7 @@ class RecursiveSimplifier
 
     //! Construct with tolerance and function to apply
     RecursiveSimplifier(F&& func, Tolerance<> tol)
-        : RecursiveSimplifier(std::forward<F>(func), tol.abs)
+        : RecursiveSimplifier(std::forward<F>(func), tol.rel)
     {
     }
 

src/orange/surf/SoftSurfaceEqual.cc

@@ -8,6 +8,7 @@
 #include "SoftSurfaceEqual.hh"
 
 #include <cmath>
+#include <limits>
 
 #include "detail/AllSurfaces.hh"
 
@@ -88,14 +89,15 @@ ORANGE_INSTANTIATE_OP(CylAligned);
  * Compare two planes for near equality.
  *
  * Consider two planes with normals \em a and \em b , without loss of
- * generality where \em a is on the \em x axis and \em b is in the \em xy
+ * generality where \em a is at \f$x=0\f$ and \em b is in the \em xy
  * plane. Suppose that to be equal, we want a intercept error of no more than
  * \f$\delta\f$ at a unit distance from the normal (since we're assuming the
  * error is based on the length scale of the problem). Taking a ray
  * from (1, 1) along (-1, 0), the distance to the plane with normal \em a is 1,
- * and the distance to plane \em b is then \f$ 1 + \delta \f$. This results
- * in a right triangle with legs of 1 and delta and an angle \f$ \theta \f$
- * at the origin, which is equal to the angle between the normals of the two
+ * and the distance to plane \em b is then \f$ 1 \pm \delta \f$. This results
+ * in a right triangle with legs of 1 and \f$ \delta \f$ and opening angle from
+ * the origin of \f$ \theta \f$, which is equal to the angle between the
+ * normals of the two
  * planes. Thus we have:
  * \f[
    \tan \theta = \frac{\delta}{1}
@@ -117,17 +119,24 @@ ORANGE_INSTANTIATE_OP(CylAligned);
  * half-space by checking for \f$ \mu > 0 \f$.
  *
  * Since this derivation is based on an absolute length scale of 1, the
- * "absolute" (rather than relative) tolerance of the plane normals is used.
- * For problems with larger or smaller length scales, that value will be scaled
- * accordingly.
+ * relative tolerance should be used.
+ *
+ * Due to floating point arithmetic, \f$mu\f$ can be slightly greater than
+ * unity, and since epsilon is often smaller than
+ * \f$\sqrt{\epsilon_\mathrm{mach}}\f$ for single precision arithmetic, the
+ * comparison here adds an extra bump to account for the precision loss.
+ *
+ * \image html orange-surface-softeq-plane.png
  */
 bool SoftSurfaceEqual::operator()(Plane const& a, Plane const& b) const
 {
     if (!soft_eq_(a.displacement(), b.displacement()))
         return false;
 
     real_type const mu = dot_product(a.normal(), b.normal());
-    return mu > 0 && (1 / ipow<2>(mu) - 1) < ipow<2>(soft_eq_.abs());
+    constexpr real_type eps_mach = std::numeric_limits<real_type>::epsilon();
+    return mu > 0
+           && (1 / ipow<2>(mu) - 1) <= (ipow<2>(soft_eq_.rel()) + eps_mach);
 }
 
 //---------------------------------------------------------------------------//
@@ -209,7 +218,7 @@ bool SoftSurfaceEqual::soft_eq_sq(real_type a, real_type b) const
 bool SoftSurfaceEqual::soft_eq_distance(Real3 const& a, Real3 const& b) const
 {
     // This is soft equal formula but using vector distance.
-    real_type rel = soft_eq_.rel() * std::fmax(norm(a), norm(b));
+    real_type rel = soft_eq_.abs() * std::fmax(norm(a), norm(b));
     return distance(a, b) < std::fmax(soft_eq_.abs(), rel);
 }
 

test/orange/orangeinp/detail/LocalSurfaceInserter.test.cc

@@ -106,6 +106,27 @@ TEST_F(LocalSurfaceInserterTest, chained_duplicates)
     EXPECT_EQ(5, surfaces.size());
 }
 
+// Replicates InfWedge.quarter_turn from convex region test
+TEST_F(LocalSurfaceInserterTest, infwedge_quadrant)
+{
+    auto tol = Tolerance<>::from_relative(1e-4);
+    LocalSurfaceInserter insert(&surfaces, tol);
+
+    constexpr real_type sqrt_half{0.70710678118655};
+    EXPECT_EQ(0, insert(PlaneY(0)).unchecked_get());
+    EXPECT_EQ(1, insert(PlaneX(0)).unchecked_get());
+    EXPECT_EQ(1, insert(PlaneX(0)).unchecked_get());
+    EXPECT_EQ(0, insert(PlaneY(0)).unchecked_get());
+    EXPECT_EQ(1, insert(PlaneX(0)).unchecked_get());
+    EXPECT_EQ(0, insert(PlaneY(0)).unchecked_get());
+    EXPECT_EQ(2, insert(Plane({sqrt_half, -sqrt_half, 0}, 0)).unchecked_get());
+    EXPECT_EQ(3, insert(Plane({sqrt_half, sqrt_half, 0}, 0)).unchecked_get());
+    EXPECT_EQ(3, insert(Plane({sqrt_half, sqrt_half, 0}, 0)).unchecked_get());
+    EXPECT_EQ(2, insert(Plane({sqrt_half, -sqrt_half, 0}, 0)).unchecked_get());
+    EXPECT_EQ(3, insert(Plane({sqrt_half, sqrt_half, 0}, 0)).unchecked_get());
+    EXPECT_EQ(2, insert(Plane({sqrt_half, -sqrt_half, 0}, 0)).unchecked_get());
+}
+
 TEST_F(LocalSurfaceInserterTest, DISABLED_performance_test)
 {
     std::mt19937 rng;

test/orange/surf/SoftSurfaceEqual.test.cc

@@ -87,10 +87,8 @@ TEST_F(SoftSurfaceEqualTest, plane)
     Real3 const n = make_unit_vector(Real3{1, 1, 0});
     Plane const ref{n, p};
 
-    if (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE)
-    {
-        EXPECT_TRUE(softeq_(ref, Plane{n, p + Real3{small, 0, 0}}));
-    }
+    EXPECT_TRUE(softeq_(ref, ref));
+    EXPECT_TRUE(softeq_(ref, Plane{n, p + Real3{small, 0, 0}}));
     EXPECT_FALSE(softeq_(ref, Plane{n, p + Real3{large, 0, 0}}));
 
     Real3 const npert = make_unit_vector(n + Real3{small, 0, 0});
@@ -102,6 +100,15 @@ TEST_F(SoftSurfaceEqualTest, plane)
     EXPECT_FALSE(softeq_(ref, Plane{make_unit_vector(Real3{1, -1, 0}), p}));
 }
 
+TEST_F(SoftSurfaceEqualTest, infwedge_quadrant)
+{
+    constexpr real_type sqrt_half{0.70710678118655};
+    Plane const p1({sqrt_half, sqrt_half, 0}, 0);
+    Plane const p2({sqrt_half, -sqrt_half, 0}, 0);
+    EXPECT_TRUE(softeq_(p1, p1));
+    EXPECT_TRUE(softeq_(p2, p2));
+}
+
 TEST_F(SoftSurfaceEqualTest, sphere)
 {
     this->check_equality_s<Sphere>({0, 1, 2}, 1);