Implement second-order convex regions (celeritas-project#1125)

* Revert "REVERTME: move second-order convex regions to a follow-on PR"
* Fix signed zeros in quadric converters
* Address review feedback

src/orange/MatrixUtils.hh

@@ -4,6 +4,7 @@
 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
 //---------------------------------------------------------------------------//
 //! \file orange/MatrixUtils.hh
+// TODO: split into BLAS and host-only utils
 //---------------------------------------------------------------------------//
 #pragma once
 

src/orange/orangeinp/ConvexRegion.cc

@@ -75,6 +75,180 @@ void Box::build(ConvexSurfaceBuilder& insert_surface) const
     insert_surface(Sense::inside, PlaneZ{hw_[Z]});
 }
 
+//---------------------------------------------------------------------------//
+// CONE
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with Z half-height and lo, hi radii.
+ */
+Cone::Cone(Real2 const& radii, real_type halfheight)
+    : radii_{radii}, hh_{halfheight}
+{
+    for (auto i : range(2))
+    {
+        CELER_VALIDATE(radii_[i] >= 0, << "negative radius: " << radii_[i]);
+    }
+    CELER_VALIDATE(radii_[0] != radii_[1], << "radii cannot be equal");
+    CELER_VALIDATE(hh_ > 0, << "nonpositive halfheight: " << hh_);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Build surfaces.
+ *
+ * The inner bounding box of a cone is determined with the following procedure:
+ * - Represent a radial slice of the cone as a right triangle with base b
+ *   (aka the higher radius) and height h (translated vanishing point)
+ * - An interior bounding box (along the xy diagonal cut!) will satisfy
+ *   r = b - tangent * z
+ * - Maximize the area of that box to obtain r = b / 2, i.e. z = h / 2
+ * - Truncate z so that it's not outside of the half-height
+ * - Project that radial slice onto the xz plane by multiplying 1/sqrt(2)
+ */
+void Cone::build(ConvexSurfaceBuilder& insert_surface) const
+{
+    // Build the bottom and top planes
+    insert_surface(Sense::outside, PlaneZ{-hh_});
+    insert_surface(Sense::inside, PlaneZ{hh_});
+
+    // Calculate the cone using lo and hi radii
+    real_type const lo{radii_[0]};
+    real_type const hi{radii_[1]};
+
+    // Arctangent of the opening angle of the cone (opposite / adjacent)
+    real_type const tangent = std::fabs(lo - hi) / (2 * hh_);
+
+    // Calculate vanishing point (origin)
+    real_type vanish_z = 0;
+    if (lo > hi)
+    {
+        // Cone opens downward (base is on bottom)
+        vanish_z = -hh_ + lo / tangent;
+        CELER_ASSERT(vanish_z > 0);
+    }
+    else
+    {
+        // Cone opens upward
+        vanish_z = hh_ - hi / tangent;
+        CELER_ASSERT(vanish_z < 0);
+    }
+
+    // Build the cone surface along the given axis
+    ConeZ cone{Real3{0, 0, vanish_z}, tangent};
+    insert_surface(cone);
+
+    // Set radial extents of exterior bbox
+    insert_surface(Sense::inside, make_xyradial_bbox(std::fmax(lo, hi)));
+
+    // Calculate the interior bounding box:
+    real_type const b = std::fmax(lo, hi);
+    real_type const h = b / tangent;
+    real_type const z = std::fmin(h / 2, 2 * hh_);
+    real_type const r = b - tangent * z;
+
+    // Now convert from "triangle above z=0" to "cone centered on z=0"
+    real_type zmin = -hh_;
+    real_type zmax = zmin + z;
+    if (lo < hi)
+    {
+        // Base is on top
+        zmax = hh_;
+        zmin = zmax - z;
+    }
+    CELER_ASSERT(zmin < zmax);
+    real_type const rbox = (constants::sqrt_two / 2) * r;
+    BBox const interior_bbox{{-rbox, -rbox, zmin}, {rbox, rbox, zmax}};
+
+    // Check that the corners are actually inside the cone
+    CELER_ASSERT(cone.calc_sense(interior_bbox.lower() * real_type(1 - 1e-5))
+                 == SignedSense::inside);
+    CELER_ASSERT(cone.calc_sense(interior_bbox.upper() * real_type(1 - 1e-5))
+                 == SignedSense::inside);
+    insert_surface(Sense::outside, interior_bbox);
+}
+
+//---------------------------------------------------------------------------//
+// CYLINDER
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with radius.
+ */
+Cylinder::Cylinder(real_type radius, real_type halfheight)
+    : radius_{radius}, hh_{halfheight}
+{
+    CELER_VALIDATE(radius_ > 0, << "nonpositive radius: " << radius_);
+    CELER_VALIDATE(hh_ > 0, << "nonpositive half-height: " << hh_);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Build surfaces.
+ */
+void Cylinder::build(ConvexSurfaceBuilder& insert_surface) const
+{
+    insert_surface(Sense::outside, PlaneZ{-hh_});
+    insert_surface(Sense::inside, PlaneZ{hh_});
+    insert_surface(CCylZ{radius_});
+}
+
+//---------------------------------------------------------------------------//
+// ELLIPSOID
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with radii.
+ */
+Ellipsoid::Ellipsoid(Real3 const& radii) : radii_{radii}
+{
+    for (auto ax : range(Axis::size_))
+    {
+        CELER_VALIDATE(radii_[to_int(ax)] > 0,
+                       << "nonpositive radius " << to_char(ax)
+                       << " axis: " << radii_[to_int(ax)]);
+    }
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Build surfaces.
+ */
+void Ellipsoid::build(ConvexSurfaceBuilder& insert_surface) const
+{
+    // Second-order coefficients are product of the other two squared radii;
+    // Zeroth-order coefficient is the product of all three squared radii
+    Real3 rsq;
+    for (auto ax : range(to_int(Axis::size_)))
+    {
+        rsq[ax] = ipow<2>(radii_[ax]);
+    }
+
+    Real3 abc{1, 1, 1};
+    real_type g = -1;
+    for (auto ax : range(to_int(Axis::size_)))
+    {
+        g *= rsq[ax];
+        for (auto nax : range(to_int(Axis::size_)))
+        {
+            if (ax != nax)
+            {
+                abc[ax] *= rsq[nax];
+            }
+        }
+    }
+
+    insert_surface(SimpleQuadric{abc, Real3{0, 0, 0}, g});
+
+    // Set exterior bbox
+    insert_surface(Sense::inside, BBox{-radii_, radii_});
+
+    // Set an interior bbox with maximum volume: a scaled inscribed cube
+    Real3 inner_radii = radii_;
+    for (real_type& r : inner_radii)
+    {
+        r *= 1 / constants::sqrt_three;
+    }
+    insert_surface(Sense::outside, BBox{-inner_radii, inner_radii});
+}
+
 //---------------------------------------------------------------------------//
 // PRISM
 //---------------------------------------------------------------------------//
@@ -145,6 +319,26 @@ void Prism::build(ConvexSurfaceBuilder& insert_surface) const
     insert_surface(Sense::outside, interior_bbox);
 }
 
+//---------------------------------------------------------------------------//
+// SPHERE
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with radius.
+ */
+Sphere::Sphere(real_type radius) : radius_{radius}
+{
+    CELER_VALIDATE(radius_ > 0, << "nonpositive radius: " << radius_);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Build surfaces.
+ */
+void Sphere::build(ConvexSurfaceBuilder& insert_surface) const
+{
+    insert_surface(SphereCentered{radius_});
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace orangeinp
 }  // namespace celeritas

src/orange/orangeinp/ConvexRegion.hh

@@ -69,6 +69,77 @@ class Box final : public ConvexRegionInterface
     Real3 hw_;
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * A closed cone along the Z axis centered on the origin.
+ *
+ * A quadric cone technically defines two opposing cones that touch at a single
+ * vanishing point, but this cone is required to be truncated so that the
+ * vanishing point is on our outside the cone.
+ *
+ * The midpoint along the Z axis of the cone is the origin. A cone is \em not
+ * allowed to have equal radii: for that, use a cylinder. However, it \em may
+ * have a single radius of zero, which puts the vanishing point on one end of
+ * the cone.
+ *
+ * This convex region, along with the Cylinder, is a base component of the
+ * G4Polycone (PCON).
+ */
+class Cone final : public ConvexRegionInterface
+{
+  public:
+    //!@{
+    //! \name Type aliases
+    using Real2 = Array<real_type, 2>;
+    //!@}
+
+  public:
+    // Construct with Z half-height and lo, hi radii
+    Cone(Real2 const& radii, real_type halfheight);
+
+    // Build surfaces
+    void build(ConvexSurfaceBuilder&) const final;
+
+  private:
+    Real2 radii_;
+    real_type hh_;
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * A Z-aligned cylinder centered on the origin.
+ */
+class Cylinder final : public ConvexRegionInterface
+{
+  public:
+    // Construct with radius
+    Cylinder(real_type radius, real_type halfheight);
+
+    // Build surfaces
+    void build(ConvexSurfaceBuilder&) const final;
+
+  private:
+    real_type radius_;
+    real_type hh_;
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * An axis-alligned ellipsoid centered at the origin.
+ */
+class Ellipsoid final : public ConvexRegionInterface
+{
+  public:
+    // Construct with radius
+    explicit Ellipsoid(Real3 const& radii);
+
+    // Build surfaces
+    void build(ConvexSurfaceBuilder&) const final;
+
+  private:
+    Real3 radii_;
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * A regular, z-extruded polygon centered on the origin.
@@ -114,6 +185,23 @@ class Prism final : public ConvexRegionInterface
     real_type orientation_;
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * A sphere centered on the origin.
+ */
+class Sphere final : public ConvexRegionInterface
+{
+  public:
+    // Construct with radius
+    explicit Sphere(real_type radius);
+
+    // Build surfaces
+    void build(ConvexSurfaceBuilder&) const final;
+
+  private:
+    real_type radius_;
+};
+
 //---------------------------------------------------------------------------//
 }  // namespace orangeinp
 }  // namespace celeritas

src/orange/orangeinp/ConvexSurfaceBuilder.hh

@@ -36,6 +36,15 @@ struct ConvexSurfaceState;
  * - Metadata combining the surfaces names with the object name
  * - Bounding boxes (interior and exterior)
  *
+ * Internally, this class uses:
+ * - \c SurfaceClipper to update bounding boxes for closed quadric surfaces
+ *   (axis aligned cylinders, spheres, planes)
+ * - \c apply_transform (which calls \c detail::SurfaceTransformer and its
+ *   siblings) to generate new surfaces based on the local transformed
+ *   coordinate system
+ * - \c RecursiveSimplifier to take transformed or user-input surfaces and
+ *   reduce them to more compact quadric representations
+ *
  * \todo Should we require that the user implicitly guarantee that the result
  * is convex, e.g. prohibit quadrics outside "saddle" points? What about a
  * torus, which (unless degenerate) is never convex? Or should we just require

src/orange/surf/detail/QuadricConeConverter.hh

@@ -120,6 +120,8 @@ QuadricConeConverter::operator()(AxisTag<T>, SimpleQuadric const& sq) const
         return {};
     }
 
+    // Clear potential signed zeros before returning
+    origin += real_type{0};
     return ConeAligned<T>::from_tangent_sq(origin, tsq);
 }
 

src/orange/surf/detail/QuadricCylConverter.hh

@@ -120,6 +120,9 @@ QuadricCylConverter::operator()(AxisTag<T>, SimpleQuadric const& sq) const
         // No real solution
         return {};
     }
+
+    // Clear potential signed zeros before returning
+    origin += real_type{0};
     return CylAligned<T>::from_radius_sq(origin, radius_sq);
 }
 

src/orange/surf/detail/QuadricSphereConverter.hh

@@ -103,6 +103,9 @@ QuadricSphereConverter::operator()(SimpleQuadric const& sq) const
         // No real solution
         return {};
     }
+
+    // Clear potential signed zeros before returning
+    origin += real_type{0};
     return Sphere::from_radius_sq(origin, radius_sq);
 }