Add external bounding box to generic trapezoid (celeritas-project#1244)

* Add exterior bounding box for gentrap

* Add bounding box grower

* Just use fmin/fmax instead of algorithm min/max

* Add 'grow' function for convex region

* Address @mrguilima feedback from celeritas-project#1237

* Updated unit test for new boundaries

src/corecel/math/Algorithms.hh

@@ -329,23 +329,17 @@ CELER_FORCEINLINE_FUNCTION void sort(RandomAccessIt first, RandomAccessIt last)
  * This function is specialized when building CUDA device code, which has
  * special intrinsics for max.
  */
-#ifndef __CUDA_ARCH__
-template<class T>
-#else
-template<class T, std::enable_if_t<!std::is_arithmetic<T>::value, bool> = true>
-#endif
+template<class T, std::enable_if_t<!std::is_floating_point<T>::value, bool> = true>
 CELER_CONSTEXPR_FUNCTION T const& max(T const& a, T const& b) noexcept
 {
     return (b > a) ? b : a;
 }
 
-#ifdef __CUDA_ARCH__
-template<class T, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true>
+template<class T, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
 CELER_CONSTEXPR_FUNCTION T max(T a, T b) noexcept
 {
-    return ::max(a, b);
+    return std::fmax(a, b);
 }
-#endif
 
 //---------------------------------------------------------------------------//
 /*!
@@ -354,23 +348,17 @@ CELER_CONSTEXPR_FUNCTION T max(T a, T b) noexcept
  * This function is specialized when building CUDA device code, which has
  * special intrinsics for min.
  */
-#ifndef __CUDA_ARCH__
-template<class T>
-#else
-template<class T, std::enable_if_t<!std::is_arithmetic<T>::value, bool> = true>
-#endif
+template<class T, std::enable_if_t<!std::is_floating_point<T>::value, bool> = true>
 CELER_CONSTEXPR_FUNCTION T const& min(T const& a, T const& b) noexcept
 {
     return (b < a) ? b : a;
 }
 
-#ifdef __CUDA_ARCH__
-template<class T, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true>
+template<class T, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
 CELER_CONSTEXPR_FUNCTION T min(T a, T b) noexcept
 {
-    return ::min(a, b);
+    return std::fmin(a, b);
 }
-#endif
 
 //---------------------------------------------------------------------------//
 /*!

src/geocel/BoundingBox.hh

@@ -7,13 +7,13 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
+#include <cmath>
 #include <type_traits>
 
 #include "corecel/Assert.hh"
 #include "corecel/Macros.hh"
 #include "corecel/Types.hh"
 #include "corecel/cont/Array.hh"
-#include "corecel/math/Algorithms.hh"
 #include "corecel/math/NumericLimits.hh"
 
 #include "Types.hh"
@@ -95,6 +95,9 @@ class BoundingBox
     CELER_CONSTEXPR_FUNCTION void
     grow(Bound bnd, Axis axis, real_type position);
 
+    // Increase the bounding box's extent on both bounds
+    CELER_CONSTEXPR_FUNCTION void grow(Axis axis, real_type position);
+
   private:
     Array<Real3, 2> points_;  //!< lo/hi points
 
@@ -257,11 +260,11 @@ BoundingBox<T>::shrink(Bound bnd, Axis axis, real_type position)
     real_type p = points_[to_int(bnd)][to_int(axis)];
     if (bnd == Bound::lo)
     {
-        p = ::celeritas::max(p, position);
+        p = std::fmax(p, position);
     }
     else
     {
-        p = ::celeritas::min(p, position);
+        p = std::fmin(p, position);
     }
     points_[to_int(bnd)][to_int(axis)] = p;
 }
@@ -280,14 +283,29 @@ BoundingBox<T>::grow(Bound bnd, Axis axis, real_type position)
     real_type p = points_[to_int(bnd)][to_int(axis)];
     if (bnd == Bound::lo)
     {
-        p = ::celeritas::min(p, position);
+        p = std::fmin(p, position);
     }
     else
     {
-        p = ::celeritas::max(p, position);
+        p = std::fmax(p, position);
     }
     points_[to_int(bnd)][to_int(axis)] = p;
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Increase (expand) the bounding box's extent along an axis.
+ *
+ * If the point is outside the box, the box is expanded so the given boundary
+ * is on that point. Otherwise no change is made.
+ */
+template<class T>
+CELER_CONSTEXPR_FUNCTION void
+BoundingBox<T>::grow(Axis axis, real_type position)
+{
+    this->grow(Bound::lo, axis, position);
+    this->grow(Bound::hi, axis, position);
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/orange/orangeinp/ConvexRegion.cc

@@ -487,6 +487,22 @@ void GenTrap::build(ConvexSurfaceBuilder& insert_surface) const
                            GeneralQuadric{abc, def, ghi, offset});
         }
     }
+
+    // Construct exterior bounding box
+    BBox exterior_bbox;
+    for (VecReal2 const* p : {&lo_, &hi_})
+    {
+        for (Real2 const& xy : *p)
+        {
+            for (auto ax : {Axis::x, Axis::y})
+            {
+                exterior_bbox.grow(ax, xy[to_int(ax)]);
+            }
+        }
+    }
+    exterior_bbox.grow(Bound::lo, Axis::z, -hz_);
+    exterior_bbox.grow(Bound::hi, Axis::z, hz_);
+    insert_surface(Sense::inside, exterior_bbox);
 }
 
 //---------------------------------------------------------------------------//

src/orange/surf/SurfaceSimplifier.hh

@@ -82,7 +82,7 @@ class SurfaceSimplifier
     // Quadric can be normalized or simplified
     Optional<SimpleQuadric, GeneralQuadric> operator()(GeneralQuadric const&);
 
-    //! Default: no simplifcation
+    //! Default: no simplification
     template<class S>
     std::variant<std::monostate> operator()(S const&) const
     {

test/orange/data/inputbuilder-incomplete-bb.org.json

@@ -210,13 +210,13 @@
 {
 "bbox": [
 [
--5.000050000000001,
--5.000050000000001,
+-2.0,
+-2.0,
 -3.0
 ],
 [
-5.000050000000001,
-5.000050000000001,
+2.0,
+2.0,
 3.0
 ]
 ],

test/orange/orangeinp/ConvexRegion.test.cc

@@ -476,8 +476,8 @@ TEST_F(GenTrapTest, trd1)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -inf, -3}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, inf, 3}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-2, -2, -3}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{2, 2, 3}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, trd2)
@@ -496,8 +496,8 @@ TEST_F(GenTrapTest, trd2)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -inf, -3}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, inf, 3}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-2, -2, -3}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{2, 2, 3}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, ppiped)
@@ -518,8 +518,8 @@ TEST_F(GenTrapTest, ppiped)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -inf, -4}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, inf, 4}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-2, -2, -4}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{2, 2, 4}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, triang_prism)
@@ -539,8 +539,8 @@ TEST_F(GenTrapTest, triang_prism)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-1, -inf, -3}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, inf, 3}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-1, -1, -3}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{2, 1, 3}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, trap_corners)
@@ -563,8 +563,8 @@ TEST_F(GenTrapTest, trap_corners)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -30, -40}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, 30, 40}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-21, -30, -40}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{21, 30, 40}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, trapezoid_trans)
@@ -588,8 +588,8 @@ TEST_F(GenTrapTest, trapezoid_trans)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -60, -40}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, 0, 40}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-51, -60, -40}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{-9, 0, 40}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, trapezoid_ccw)
@@ -611,8 +611,8 @@ TEST_F(GenTrapTest, trapezoid_ccw)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -30, -40}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, 30, 40}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-21, -30, -40}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{21, 30, 40}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, trap_theta)
@@ -632,8 +632,8 @@ TEST_F(GenTrapTest, trap_theta)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -20, -40}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, 20, 40}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-50, -20, -40}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{50, 20, 40}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, trap_thetaphi)
@@ -653,8 +653,8 @@ TEST_F(GenTrapTest, trap_thetaphi)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-10, -inf, -40}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{10, inf, 40}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-10, -60, -40}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{10, 60, 40}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, trap_g4)
@@ -679,8 +679,10 @@ TEST_F(GenTrapTest, trap_g4)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -inf, -4}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, inf, 4}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-1.95920952072934, -2.93923101204883, -4}),
+                       result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{2.64848563385739, 3.06076898795117, 4}),
+                       result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, trap_full)
@@ -700,8 +702,10 @@ TEST_F(GenTrapTest, trap_full)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -inf, -40}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, inf, 40}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-40.2842712474619, -48.2842712474619, -40}),
+                       result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{40.2842712474619, 48.2842712474619, 40}),
+                       result.exterior.upper());
 }
 
 // TODO: this should be valid
@@ -737,8 +741,8 @@ TEST_F(GenTrapTest, full)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -2, -4}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, inf, 4}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-2, -2, -4}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{2, 2, 4}), result.exterior.upper());
 }
 
 TEST_F(GenTrapTest, full2)
@@ -758,33 +762,28 @@ TEST_F(GenTrapTest, full2)
     EXPECT_EQ(expected_node, result.node);
     EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
     EXPECT_FALSE(result.interior) << result.interior;
-    EXPECT_VEC_SOFT_EQ((Real3{-inf, -20, -40}), result.exterior.lower());
-    EXPECT_VEC_SOFT_EQ((Real3{inf, 20, 40}), result.exterior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-52, -20, -40}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{54, 20, 40}), result.exterior.upper());
 }
 
+/*!
+ * Test deduplication of two opposing quadric surfaces.
+ *
+ * \verbatim
+ * Lower polygons:      Upper polygons:
+ *
+ * x=-1      x=1           x=-0.5
+ * +----+----+ y=1      +--+------+ y=1
+ * |    |    |          |   \     |
+ * |    |  R |          |    \  R |
+ * |  L |    |          |  L  \   |
+ * |    |    |          |      \  |
+ * +----+----+ y=-1     +-------+-+ y=-1
+ *      x=0                     x=0.5
+ * \endverbatim
+ */
 TEST_F(GenTrapTest, adjacent_twisted)
 {
-    /* Lower polygons:
-     *
-     * x=-1      x=1
-     * +----+----+ y=1
-     * |    |    |
-     * |  L |  R |
-     * |    |    |
-     * |    |    |
-     * +----+----+ y=-1
-     *
-     *
-     * Upper polygons:
-     *    x=-0.5
-     * +--+------+ y=1
-     * |   \     |
-     * |    \  R |
-     * |  L  \   |
-     * |      \  |
-     * +-------+-+ y=-1
-     *         x=0.5
-     */
     {
         // Left
         auto result
@@ -796,7 +795,7 @@ TEST_F(GenTrapTest, adjacent_twisted)
 
         EXPECT_EQ(expected_node, result.node);
         EXPECT_VEC_SOFT_EQ((Real3{-1, -1, -1}), result.exterior.lower());
-        EXPECT_VEC_SOFT_EQ((Real3{inf, 1, 1}), result.exterior.upper());
+        EXPECT_VEC_SOFT_EQ((Real3{0.5, 1, 1}), result.exterior.upper());
     }
     {
         // Right
@@ -808,7 +807,7 @@ TEST_F(GenTrapTest, adjacent_twisted)
         static char const expected_node[] = "all(+0, -1, +2, +3, -4, -6)";
 
         EXPECT_EQ(expected_node, result.node);
-        EXPECT_VEC_SOFT_EQ((Real3{-inf, -1, -1}), result.exterior.lower());
+        EXPECT_VEC_SOFT_EQ((Real3{-0.5, -1, -1}), result.exterior.lower());
         EXPECT_VEC_SOFT_EQ((Real3{1, 1, 1}), result.exterior.upper());
     }
     {
@@ -821,7 +820,7 @@ TEST_F(GenTrapTest, adjacent_twisted)
         static char const expected_node[] = "all(+0, -1, +7, -8, -9, +10)";
 
         EXPECT_EQ(expected_node, result.node);
-        EXPECT_VEC_SOFT_EQ((Real3{-inf, -2, -1}), result.exterior.lower());
+        EXPECT_VEC_SOFT_EQ((Real3{-1, -2, -1}), result.exterior.lower());
         EXPECT_VEC_SOFT_EQ((Real3{2, 2, 1}), result.exterior.upper());
     }
 

test/orange/orangeinp/UnitProto.test.cc

@@ -466,7 +466,8 @@ class InputBuilderTest : public UnitProtoTest
 
         std::stringstream expected;
         expected << ref.rdbuf();
-        EXPECT_JSON_EQ(expected.str(), actual.str());
+        EXPECT_JSON_EQ(expected.str(), actual.str())
+            << "update the file at " << ref_path;
     }
 };