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

@@ -126,6 +126,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
 };
 

src/orange/orangeinp/IntersectRegion.cc

@@ -23,6 +23,7 @@
 #include "orange/surf/PlaneAligned.hh"
 #include "orange/surf/SimpleQuadric.hh"
 #include "orange/surf/SphereCentered.hh"
+#include "orange/surf/Involute.hh"
 
 #include "IntersectSurfaceBuilder.hh"
 #include "ObjectIO.json.hh"
@@ -90,6 +91,83 @@ void Box::output(JsonPimpl* j) const
     to_json_pimpl(j, *this);
 }
 
+//---------------------------------------------------------------------------//
+// Involute
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with prarameters and half height.
+ */
+Involute::Involute(Real3 const& radii, Real2 const& displacement, 
+                   Sign sign, real_type halfheight) : 
+radii_(radii), a_(displacement), tmin_(), tmax_(), sign_(sign), hh_{halfheight}
+{
+    for (auto i : range(3))
+    {
+        CELER_VALIDATE(radii_[i] >= 0, << "negative radius: " << radii_[i]);
+    }
+
+    tmin_ = std::sqrt(clamp_to_nonneg(ipow<2>(radii_[1]/radii_[0]) - 1));
+    tmax_ = std::sqrt(clamp_to_nonneg(ipow<2>(radii_[2]/radii_[0]) - 1));
+
+    for (auto i : range(2))
+    {
+        CELER_VALIDATE(a_[i] >= 0, << "negative displacment angle: " << a_[i]);
+    }
+    CELER_VALIDATE(a_[1] > a_[0], 
+                   << "nonpositive delta displacment: " << a_[1]-a_[0]);
+    CELER_VALIDATE(hh_ > 0, << "nonpositive half-height: " << hh_); 
+
+    CELER_VALIDATE(tmin_ > 0, << "nonpositive inner radius: " << tmin_); 
+    CELER_VALIDATE(tmax_ > tmin_, << "outer radius smaller than inner radius: " 
+                                  << tmax_); 
+    CELER_VALIDATE(tmax_ <  tmin_ + 2 * constants::pi - (a_[1]-a_[0]),
+                   << "radial bounds result in angular overlaped: " 
+                   << tmin_ + 2 * constants::pi - (a_[1]-a_[0]) - tmax_);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Build surfaces.
+ */
+void Involute::build(IntersectSurfaceBuilder& insert_surface) const
+{
+    if (sign_ == Sign::clockwise)
+    {
+        insert_surface(Sense::outside, celeritas::Involute{{0,0}, radii_[0], 
+                                                eumod(a_[0], 2*constants::pi), 
+                                                        sign_, tmin_, 
+                                                        tmax_+a_[1]-a_[0]});
+        insert_surface(Sense::inside, celeritas::Involute{{0,0}, radii_[0], 
+                                                eumod(a_[1], 2*constants::pi),
+                                                        sign_, tmin_, 
+                                                        tmax_+a_[1]-a_[0]});
+    } 
+    else
+    {
+        insert_surface(Sense::inside, celeritas::Involute{{0,0}, radii_[0], 
+                                                eumod(a_[0], 2*constants::pi), 
+                                                        sign_, tmin_, 
+                                                        tmax_+a_[1]-a_[0]});
+        insert_surface(Sense::outside, celeritas::Involute{{0,0}, radii_[0], 
+                                                eumod(a_[1], 2*constants::pi),
+                                                        sign_, tmin_, 
+                                                        tmax_+a_[1]-a_[0]});
+    }
+    insert_surface(Sense::outside, PlaneZ{-hh_});
+    insert_surface(Sense::inside, PlaneZ{hh_});
+    insert_surface(Sense::inside, CCylZ{radii_[2]});
+    insert_surface(Sense::outside, CCylZ{radii_[1]});
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Write output to the given JSON object.
+ */
+void Involute::output(JsonPimpl* j) const
+{
+    to_json_pimpl(j, *this);
+}
+
 //---------------------------------------------------------------------------//
 // CONE
 //---------------------------------------------------------------------------//

src/orange/orangeinp/IntersectRegion.hh

@@ -13,6 +13,8 @@
 #include "corecel/math/Turn.hh"
 #include "orange/OrangeTypes.hh"
 
+#include "../surf/detail/InvoluteSolver.hh"
+
 namespace celeritas
 {
 struct JsonPimpl;
@@ -84,6 +86,49 @@ class Box final : public IntersectRegionInterface
     Real3 hw_;
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * An involute centered on the origin.
+ *
+ * \note Be aware there's also an involute *surface* at orange/surf/Involute.hh 
+ * in a different namespace.
+ */
+class Involute final : public IntersectRegionInterface
+{
+  public:
+    //! Enum defining chirality of involute
+    using Sign = celeritas::detail::InvoluteSolver::Sign;
+
+
+    // Construct with radius
+    explicit Involute(Real3 const& radii, Real2 const& displacement, 
+                      Sign sign, real_type halfheight);
+
+    // Build surfaces
+    void build(IntersectSurfaceBuilder&) const final;
+
+    // Output to JSON
+    void output(JsonPimpl*) const final;
+
+    //// ACCESSORS ////
+
+    //! Radius
+    Real3 radii() const { return radii_; }
+    Real2 a() const { return a_; }
+    real_type tmin() const { return tmin_; }
+    real_type tmax() const { return tmax_; }
+    Sign sign() const { return sign_; }
+    real_type hh() const { return hh_; }
+
+  private:
+    Real3 radii_;
+    Real2 a_;
+    real_type tmin_;
+    real_type tmax_;
+    Sign sign_;
+    real_type hh_;
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * A closed cone along the Z axis 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

@@ -212,6 +212,16 @@ 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, a),
+         SIO_ATTR_PAIR(cr, tmin),
+         SIO_ATTR_PAIR(cr, tmax),
+         SIO_ATTR_PAIR(cr, sign),
+         SIO_ATTR_PAIR(cr, hh)};
+}
 //!@}
 
 //---------------------------------------------------------------------------//

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

@@ -47,6 +47,8 @@ struct SurfaceHashPoint
         return std::sqrt(s.radius_sq());
     }
 
+    real_type operator()(Involute const& s) const { return s.r_b(); }
+
     template<Axis T>
     real_type operator()(CylAligned<T> const& s) const
     {

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

@@ -34,7 +34,7 @@ Involute::Involute(Real2 const& origin,
     CELER_EXPECT(tmin_ >= 0);
     CELER_EXPECT(tmax_ > tmin_ && tmax_ < 2 * constants::pi + tmin_);
 
-    if (sign)
+    if (sign == Sign::clockwise)
     {
         a_ = constants::pi - a_;
         r_b_ = -r_b_;

src/orange/surf/Involute.hh

@@ -67,9 +67,9 @@ class Involute
     //// CLASS ATTRIBUTES ////
 
     // Surface type identifier
-    static SurfaceType surface_type()
+    static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type()
     {
-        CELER_NOT_IMPLEMENTED("runtime involute");
+        return SurfaceType::inv;
     }
 
     // Safety

src/orange/surf/SoftSurfaceEqual.cc

@@ -68,6 +68,21 @@ bool SoftSurfaceEqual::operator()(SphereCentered const& a,
     return this->soft_eq_sq(a.radius_sq(), b.radius_sq());
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Compare two centered involutes for near equality.
+ */
+bool SoftSurfaceEqual::operator()(Involute const& a, Involute const& b) const
+{
+    return this->soft_eq_sq(a.r_b(), b.r_b()) && this->soft_eq_sq(a.a(), b.a())
+           && this->soft_eq_sq(a.sign() == Involute::Sign::clockwise, 
+                               b.sign() == Involute::Sign::clockwise)
+           && this->soft_eq_sq(a.tmin(), b.tmin())
+           && this->soft_eq_sq(a.tmax(), b.tmax())
+           && this->soft_eq_sq(a.origin()[0], b.origin()[0])
+           && this->soft_eq_sq(a.origin()[1], b.origin()[1]);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Compare two aligned cylinders for near equality.

src/orange/surf/SoftSurfaceEqual.hh

@@ -66,6 +66,8 @@ class SoftSurfaceEqual
 
     bool operator()(Sphere const&, Sphere const&) const;
 
+    bool operator()(Involute const&, Involute const&) const;
+
     template<Axis T>
     bool operator()(ConeAligned<T> const&, ConeAligned<T> const&) const;
 

src/orange/surf/SurfaceClipper.cc

@@ -167,6 +167,16 @@ void SurfaceClipper::operator()(GeneralQuadric const&) const
     *int_ = BoundingBox{};
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Clip the bounding boxes to an involute.
+ */
+void SurfaceClipper::operator()(Involute const&) const
+{
+    // We no longer can guarantee any point being inside the shape; reset it
+    *int_ = BoundingBox{};
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Clip a variant surface.

src/orange/surf/SurfaceClipper.hh

@@ -70,6 +70,8 @@ class SurfaceClipper
 
     void operator()(GeneralQuadric const&) const;
 
+    void operator()(Involute const&) const;
+
     // Apply to a surface with unknown type
     void operator()(VariantSurface const& surf) const;
 

src/orange/surf/SurfaceIO.cc

@@ -113,7 +113,8 @@ std::ostream& operator<<(std::ostream& os, SphereCentered const& s)
 std::ostream& operator<<(std::ostream& os, Involute const& s)
 {
     os << "Involute: r, a, sign, tmin, tmax =" << s.r_b() << ' ' << s.a()
-       << ' ' << s.sign() << ' ' << s.tmin() << ' ' << s.tmax() << " at "
+       << ' ' << (s.sign() == Involute::Sign::clockwise) << ' ' << s.tmin() 
+       << ' ' << s.tmax() << " at "
        << s.origin();
     return os;
 }

src/orange/surf/SurfaceTypeTraits.hh

@@ -51,6 +51,7 @@ ORANGE_SURFACE_TRAITS(ky,  ConeAligned<Axis::y>);
 ORANGE_SURFACE_TRAITS(kz,  ConeAligned<Axis::z>);
 ORANGE_SURFACE_TRAITS(sq,  SimpleQuadric);
 ORANGE_SURFACE_TRAITS(gq,  GeneralQuadric);
+ORANGE_SURFACE_TRAITS(inv, Involute);
 // clang-format on
 
 #undef ORANGE_SURFACE_TRAITS
@@ -90,6 +91,7 @@ visit_surface_type(F&& func, SurfaceType st)
         ORANGE_ST_VISIT_CASE(kz);
         ORANGE_ST_VISIT_CASE(sq);
         ORANGE_ST_VISIT_CASE(gq);
+        ORANGE_ST_VISIT_CASE(inv);
         default:
             CELER_ASSERT_UNREACHABLE();
     }

src/orange/surf/detail/AllSurfaces.hh

@@ -12,6 +12,7 @@
 #include "../CylAligned.hh"
 #include "../CylCentered.hh"
 #include "../GeneralQuadric.hh"
+#include "../Involute.hh"
 #include "../Plane.hh"
 #include "../PlaneAligned.hh"
 #include "../SimpleQuadric.hh"

src/orange/surf/detail/InvoluteSolver.hh

@@ -50,7 +50,7 @@ class InvoluteSolver
     using Real2 = Array<real_type, 2>;
 
     //! Enum defining chirality of involute
-    enum Sign : bool
+    enum class Sign : bool
     {
         counterclockwise = 0,
         clockwise,  //!< Apply symmetry when solving
@@ -152,7 +152,7 @@ InvoluteSolver::operator()(Real3 const& pos,
     real_type v = dir[1];
 
     // Mirror systemm for counterclockwise involutes
-    if (sign_)
+    if (sign_==Sign::clockwise)
     {
         x = -x;
         u = -u;