Integrate involute surface into ORANGE construction

src/orange/OrangeTypes.cc

@@ -125,6 +125,7 @@ char const* to_cstring(SurfaceType value)
         "kz",
         "sq",
         "gq",
+        "inv",
         // clang-format on
     };
     return to_cstring_impl(value);

src/orange/OrangeTypes.hh

@@ -102,8 +102,13 @@ enum class Sense : bool
 /*!
  * Enumeration for mapping surface classes to integers.
  *
- * These are ordered by number of coefficients needed for their representation:
- * 1 for `p.|sc|c.c`, 3 for `c.`, 4 for `[ps]|k.`, 7 for `sq`, and 10 for `gq`.
+ * These are ordered roughly by complexity. The storage requirement for
+ * corresponding surfaces are:
+ * - 1 for `p.|sc|c.c`,
+ * - 3 for `c.`,
+ * - 4 for `[ps]|k.`,
+ * - 7 for `sq`, and
+ * - 10 for `gq`.
  *
  * See \c orange/surf/SurfaceTypeTraits.hh for how these map to classes.
  */
@@ -126,6 +131,7 @@ enum class SurfaceType : unsigned char
     kz,  //!< Cone parallel to Z axis
     sq,  //!< Simple quadric
     gq,  //!< General quadric
+    inv,  //!< Involute
     size_  //!< Sentinel value for number of surface types
 };
 
@@ -211,6 +217,16 @@ enum class BoundaryResult : bool
     exiting = true
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * Chirality of a twirly object (currently only Involute).
+ */
+enum class Chirality : bool
+{
+    left,  //!< Sinistral, spiraling counterclockwise
+    right,  //!< Dextral, spiraling clockwise
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * Volume logic encoding.

src/orange/detail/SurfacesRecordBuilder.cc

@@ -39,6 +39,13 @@ auto SurfacesRecordBuilder::operator()(VecSurface const& surfaces) -> result_typ
 
     // Functor to save the surface type and data, and the data offset
     auto emplace_surface = [this](auto&& s) {
+        if constexpr (std::remove_reference_t<decltype(s)>::surface_type()
+                      == SurfaceType::inv)
+        {
+            // See discussion on
+            // https://github.com/celeritas-project/celeritas/pull/1342
+            CELER_NOT_IMPLEMENTED("runtime involute support");
+        }
         types_.push_back(s.surface_type());
         auto data = s.data();
         auto real_range = reals_.insert_back(data.begin(), data.end());

src/orange/orangeinp/IntersectRegion.cc

@@ -17,9 +17,11 @@
 #include "corecel/math/SoftEqual.hh"
 #include "geocel/BoundingBox.hh"
 #include "geocel/Types.hh"
+#include "orange/OrangeTypes.hh"
 #include "orange/orangeinp/detail/PolygonUtils.hh"
 #include "orange/surf/ConeAligned.hh"
 #include "orange/surf/CylCentered.hh"
+#include "orange/surf/Involute.hh"
 #include "orange/surf/PlaneAligned.hh"
 #include "orange/surf/SimpleQuadric.hh"
 #include "orange/surf/SphereCentered.hh"
@@ -647,6 +649,82 @@ void InfWedge::output(JsonPimpl* j) const
     to_json_pimpl(j, *this);
 }
 
+//---------------------------------------------------------------------------//
+// Involute
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with prarameters and half height.
+ */
+Involute::Involute(Real3 const& radii,
+                   Real2 const& displacement,
+                   Chirality sign,
+                   real_type halfheight)
+    : radii_(radii), a_(displacement), t_bounds_(), sign_(sign), hh_{halfheight}
+{
+    CELER_VALIDATE(radii_[0] > 0,
+                   << "nonpositive involute radius: " << radii_[0]);
+    CELER_VALIDATE(radii_[1] > radii_[0],
+                   << "inner cylinder radius " << radii_[1]
+                   << " is not greater than involute radius " << radii_[0]);
+    CELER_VALIDATE(radii_[2] > radii_[1],
+                   << "outer cylinder radius " << radii_[2]
+                   << " is not greater than inner cyl radius " << radii_[1]);
+
+    CELER_VALIDATE(a_[1] > a_[0],
+                   << "nonpositive delta displacment: " << a_[1] - a_[0]);
+    CELER_VALIDATE(hh_ > 0, << "nonpositive half-height: " << hh_);
+
+    for (auto i : range(2))
+    {
+        t_bounds_[i] = std::sqrt(
+            clamp_to_nonneg(ipow<2>(radii_[i + 1] / radii_[0]) - 1));
+    }
+    auto outer_isect = t_bounds_[0] + 2 * constants::pi - (a_[1] - a_[0]);
+    CELER_VALIDATE(t_bounds_[1] < outer_isect,
+                   << "radial bounds result in angular overlap: "
+                   << outer_isect - t_bounds_[1]);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Build surfaces.
+ */
+void Involute::build(IntersectSurfaceBuilder& insert_surface) const
+{
+    using InvSurf = ::celeritas::Involute;
+
+    insert_surface(Sense::outside, PlaneZ{-hh_});
+    insert_surface(Sense::inside, PlaneZ{hh_});
+    insert_surface(Sense::outside, CCylZ{radii_[1]});
+    insert_surface(Sense::inside, CCylZ{radii_[2]});
+    // Make an inside and outside involute
+    Real2 const xy{0, 0};
+    auto sense = (sign_ == Chirality::right ? Sense::outside : Sense::inside);
+    static char const* names[] = {"invl", "invr"};
+
+    for (auto i : range(2))
+    {
+        insert_surface(sense,
+                       InvSurf{xy,
+                               radii_[0],
+                               eumod(a_[i], 2 * constants::pi),
+                               sign_,
+                               t_bounds_[0],
+                               t_bounds_[1] + a_[1] - a_[0]},
+                       std::string{names[i]});
+        sense = flip_sense(sense);
+    }
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Write output to the given JSON object.
+ */
+void Involute::output(JsonPimpl* j) const
+{
+    to_json_pimpl(j, *this);
+}
+
 //---------------------------------------------------------------------------//
 // PARALLELEPIPED
 //---------------------------------------------------------------------------//
@@ -687,18 +765,18 @@ void Parallelepiped::build(IntersectSurfaceBuilder& insert_surface) const
     constexpr auto Y = to_int(Axis::y);
     constexpr auto Z = to_int(Axis::z);
 
-    // cache trigonometric values
+    // Cache trigonometric values
     real_type sinth, costh, sinphi, cosphi, sinal, cosal;
     sincos(theta_, &sinth, &costh);
     sincos(phi_, &sinphi, &cosphi);
     sincos(alpha_, &sinal, &cosal);
 
-    // base vectors
+    // Base vectors
     auto a = hpr_[X] * Real3{1, 0, 0};
     auto b = hpr_[Y] * Real3{sinal, cosal, 0};
     auto c = hpr_[Z] * Real3{sinth * cosphi, sinth * sinphi, costh};
 
-    // positioning the planes
+    // Position the planes
     auto xnorm = make_unit_vector(cross_product(b, c));
     auto ynorm = make_unit_vector(cross_product(c, a));
     auto xoffset = dot_product(a, xnorm);

src/orange/orangeinp/IntersectRegion.hh

@@ -317,6 +317,53 @@ class InfWedge final : public IntersectRegionInterface
     Turn interior_;
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * An involute "blade" centered on the origin.
+ *
+ * This is the intersection of two parallel involutes with a cylindrical shell.
+ * The three radii, which must be in ascending order, are that of the involute,
+ * the inner cylinder, and the outer cylinder.
+ *
+ * The "chirality" of the involute is viewed from the \em +z axis looking down:
+ * whether it spirals to the right or left.
+ */
+class Involute final : public IntersectRegionInterface
+{
+  public:
+    // Construct with radius
+    explicit Involute(Real3 const& radii,
+                      Real2 const& displacement,
+                      Chirality chirality,
+                      real_type halfheight);
+
+    // Build surfaces
+    void build(IntersectSurfaceBuilder&) const final;
+
+    // Output to JSON
+    void output(JsonPimpl*) const final;
+
+    //// ACCESSORS ////
+
+    //! Radii: Rdius of involute, minimum radius, maximum radius
+    Real3 const& radii() const { return radii_; }
+    //! Displacement angle
+    Real2 const& displacement_angle() const { return a_; }
+    //!  Angular bounds of involute
+    Real2 const& t_bounds() const { return t_bounds_; }
+    //! Chirality of involute: turning left or right
+    Chirality chirality() const { return sign_; }
+    //! Halfheight
+    real_type halfheight() const { return hh_; }
+
+  private:
+    Real3 radii_;
+    Real2 a_;
+    Real2 t_bounds_;
+    Chirality sign_;
+    real_type hh_;
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * A general parallelepiped centered on the origin.

src/orange/orangeinp/IntersectSurfaceBuilder.cc

@@ -245,6 +245,7 @@ CSB_INSTANTIATE(ConeAligned<Axis::y>);
 CSB_INSTANTIATE(ConeAligned<Axis::z>);
 CSB_INSTANTIATE(SimpleQuadric);
 CSB_INSTANTIATE(GeneralQuadric);
+CSB_INSTANTIATE(Involute);
 #undef CSB_INSTANTIATE
 //! \endcond
 

src/orange/orangeinp/ObjectIO.json.cc

@@ -209,6 +209,15 @@ void to_json(nlohmann::json& j, Sphere const& cr)
 {
     j = {{"_type", "sphere"}, SIO_ATTR_PAIR(cr, radius)};
 }
+void to_json(nlohmann::json& j, Involute const& cr)
+{
+    j = {{"_type", "involute"},
+         SIO_ATTR_PAIR(cr, radii),
+         SIO_ATTR_PAIR(cr, displacement_angle),
+         SIO_ATTR_PAIR(cr, t_bounds),
+         SIO_ATTR_PAIR(cr, chirality),
+         SIO_ATTR_PAIR(cr, halfheight)};
+}
 //!@}
 
 //---------------------------------------------------------------------------//

src/orange/orangeinp/ObjectIO.json.hh

@@ -43,6 +43,7 @@ class InfWedge;
 class Parallelepiped;
 class Prism;
 class Sphere;
+class Involute;
 
 //---------------------------------------------------------------------------//
 
@@ -74,6 +75,7 @@ void to_json(nlohmann::json& j, InfWedge const& cr);
 void to_json(nlohmann::json& j, Parallelepiped const& cr);
 void to_json(nlohmann::json& j, Prism const& cr);
 void to_json(nlohmann::json& j, Sphere const& cr);
+void to_json(nlohmann::json& j, Involute const& cr);
 
 //---------------------------------------------------------------------------//
 }  // namespace orangeinp

src/orange/orangeinp/detail/LocalSurfaceInserter.cc

@@ -203,6 +203,7 @@ LSI_INSTANTIATE(ConeAligned<Axis::y>);
 LSI_INSTANTIATE(ConeAligned<Axis::z>);
 LSI_INSTANTIATE(SimpleQuadric);
 LSI_INSTANTIATE(GeneralQuadric);
+LSI_INSTANTIATE(Involute);
 
 #undef LSI_INSTANTIATE
 //! \endcond

src/orange/orangeinp/detail/SurfaceHashPoint.hh

@@ -66,6 +66,11 @@ struct SurfaceHashPoint
         return norm(s.origin());
     }
 
+    real_type operator()(Involute const& s) const
+    {
+        return s.displacement_angle();
+    }
+
     real_type operator()(SimpleQuadric const& s) const
     {
         return std::sqrt(s.zeroth());

src/orange/surf/FaceNamer.cc

@@ -111,6 +111,15 @@ std::string FaceNamer::Impl::operator()(Sphere const&) const
     return "s";
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Construct a name for an involute.
+ */
+std::string FaceNamer::Impl::operator()(Involute const&) const
+{
+    return "inv";
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Construct a name for a cone.

src/orange/surf/FaceNamer.hh

@@ -87,6 +87,8 @@ class FaceNamer
         std::string operator()(SimpleQuadric const&) const;
 
         std::string operator()(GeneralQuadric const&) const;
+
+        std::string operator()(Involute const&) const;
     };
 };
 

src/orange/surf/Involute.cc

@@ -23,7 +23,7 @@ namespace celeritas
 Involute::Involute(Real2 const& origin,
                    real_type radius,
                    real_type displacement,
-                   Sign sign,
+                   Chirality sign,
                    real_type tmin,
                    real_type tmax)
     : origin_(origin), r_b_(radius), a_(displacement), tmin_(tmin), tmax_(tmax)
@@ -33,7 +33,7 @@ Involute::Involute(Real2 const& origin,
     CELER_EXPECT(tmin_ >= 0);
     CELER_EXPECT(tmax_ > tmin_ && tmax_ < 2 * constants::pi + tmin_);
 
-    if (sign)
+    if (sign == Chirality::right)
     {
         a_ = constants::pi - a_;
         r_b_ = -r_b_;