Add CSG wedge and helper functions (#1142)

* Add  eumod helper function
* Add infinite wedge and test
* Move shape construction into helper function
* Insert period after prefix for face builder
* Restrict start angle to [0,1)

src/corecel/math/Algorithms.hh

@@ -497,16 +497,45 @@ template<class T>
 /*!
  * Calculate the difference of squares \f$ a^2 - b^2 \f$.
  *
- * The expression \f$ a^2 - b^2 \f$ exhibits catastrophic cancellation when \f$
- * a \f$ and \f$ b \f$ are close in magnitude. This can be avoided by
- * rearranging the formula as \f$ (a - b)(a + b) \f$.
+ * This calculation exchanges one multiplication for one addition, but it does
+ * not increase the accuracy of the computed result. It is used
+ * occasionally in Geant4 but is likely a premature optimization... see
+ * https://github.com/celeritas-project/celeritas/pull/1082
  */
 template<class T>
 CELER_CONSTEXPR_FUNCTION T diffsq(T a, T b)
 {
     return (a - b) * (a + b);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate the Euclidian modulus of two numbers.
+ *
+ * If both numbers are positive, this should be the same as fmod. If the
+ * sign of the remainder and denominator don't match, the remainder will be
+ * remapped so that it is between zero and the denominator.
+ *
+ * This function is useful for normalizing user-provided angles.
+ */
+template<class T, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
+CELER_CONSTEXPR_FUNCTION T eumod(T numer, T denom)
+{
+    T r = std::fmod(numer, denom);
+    if (r < 0)
+    {
+        if (denom >= 0)
+        {
+            r += denom;
+        }
+        else
+        {
+            r -= denom;
+        }
+    }
+    return r;
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Double-precision math constant (POSIX derivative).

src/orange/CMakeLists.txt

@@ -34,6 +34,7 @@ list(APPEND SOURCES
   orangeinp/Shape.cc
   orangeinp/Transformed.cc
   orangeinp/detail/BoundingZone.cc
+  orangeinp/detail/BuildConvexRegion.cc
   orangeinp/detail/ConvexSurfaceState.cc
   orangeinp/detail/CsgUnitBuilder.cc
   orangeinp/detail/LocalSurfaceInserter.cc

src/orange/orangeinp/ConvexRegion.cc

@@ -290,6 +290,51 @@ void Ellipsoid::output(JsonPimpl* j) const
     to_json_pimpl(j, *this);
 }
 
+//---------------------------------------------------------------------------//
+// INFWEDGE
+//---------------------------------------------------------------------------//
+/*!
+ * Construct from a starting angle and interior angle.
+ */
+InfWedge::InfWedge(Turn start, Turn interior)
+    : start_{start}, interior_{interior}
+{
+    CELER_VALIDATE(start_ >= zero_quantity() && start_ < Turn{1},
+                   << "invalid start angle " << start_.value()
+                   << " [turns]: must be in the range [0, 1)");
+    CELER_VALIDATE(interior_ > zero_quantity() && interior_ <= Turn{0.5},
+                   << "invalid interior wedge angle " << interior.value()
+                   << " [turns]: must be in the range (0, 0.5]");
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Build surfaces.
+ *
+ * Both planes should point "outward" to the wedge. In the degenerate case of
+ * interior = 0.5 we rely on CSG object deduplication.
+ */
+void InfWedge::build(ConvexSurfaceBuilder& insert_surface) const
+{
+    real_type sinstart, cosstart, sinend, cosend;
+    sincos(start_, &sinstart, &cosstart);
+    sincos(start_ + interior_, &sinend, &cosend);
+
+    insert_surface(Sense::inside, Plane{Real3{sinstart, -cosstart, 0}, 0.0});
+    insert_surface(Sense::outside, Plane{Real3{sinend, -cosend, 0}, 0.0});
+
+    // TODO: restrict bounding boxes, at least eliminating two quadrants...
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Write output to the given JSON object.
+ */
+void InfWedge::output(JsonPimpl* j) const
+{
+    to_json_pimpl(j, *this);
+}
+
 //---------------------------------------------------------------------------//
 // PRISM
 //---------------------------------------------------------------------------//

src/orange/orangeinp/ConvexRegion.hh

@@ -10,6 +10,7 @@
 
 #include "corecel/Macros.hh"
 #include "corecel/cont/Array.hh"
+#include "corecel/math/Turn.hh"
 #include "orange/OrangeTypes.hh"
 
 namespace celeritas
@@ -183,6 +184,40 @@ class Ellipsoid final : public ConvexRegionInterface
     Real3 radii_;
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * An open wedge shape from the Z axis.
+ *
+ * The wedge is defined by an interior angle that *must* be less than or equal
+ * to 180 degrees (half a turn) and *must* be more than zero. It can be
+ * subtracted, or its negation can be subtracted. The start angle is mapped
+ * onto [0, 1) on construction.
+ */
+class InfWedge final : public ConvexRegionInterface
+{
+  public:
+    // Construct from a starting angle and interior angle
+    InfWedge(Turn start, Turn interior);
+
+    // Build surfaces
+    void build(ConvexSurfaceBuilder&) const final;
+
+    // Output to JSON
+    void output(JsonPimpl*) const final;
+
+    //// ACCESSORS ////
+
+    //! Starting angle
+    Turn start() const { return start_; }
+
+    //! Interior angle
+    Turn interior() const { return interior_; }
+
+  private:
+    Turn start_;
+    Turn interior_;
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * A regular, z-extruded polygon centered on the origin.

src/orange/orangeinp/ObjectIO.json.cc

@@ -161,6 +161,12 @@ void to_json(nlohmann::json& j, Ellipsoid const& cr)
 {
     j = {{"_type", "ellipsoid"}, SIO_ATTR_PAIR(cr, radii)};
 }
+void to_json(nlohmann::json& j, InfWedge const& cr)
+{
+    j = {{"_type", "infwedge"},
+         {"start", cr.start().value()},
+         {"interior", cr.interior().value()}};
+}
 void to_json(nlohmann::json& j, Prism const& cr)
 {
     j = {{"_type", "prism"},

src/orange/orangeinp/ObjectIO.json.hh

@@ -29,6 +29,7 @@ class Box;
 class Cone;
 class Cylinder;
 class Ellipsoid;
+class InfWedge;
 class Prism;
 class Sphere;
 
@@ -50,6 +51,7 @@ void to_json(nlohmann::json& j, Box const& cr);
 void to_json(nlohmann::json& j, Cone const& cr);
 void to_json(nlohmann::json& j, Cylinder const& cr);
 void to_json(nlohmann::json& j, Ellipsoid const& cr);
+void to_json(nlohmann::json& j, InfWedge const& cr);
 void to_json(nlohmann::json& j, Prism const& cr);
 void to_json(nlohmann::json& j, Sphere const& cr);
 

src/orange/orangeinp/Shape.cc

@@ -9,12 +9,7 @@
 
 #include "corecel/io/JsonPimpl.hh"
 
-#include "ConvexSurfaceBuilder.hh"
-#include "CsgTreeUtils.hh"
-
-#include "detail/ConvexSurfaceState.hh"
-#include "detail/CsgUnitBuilder.hh"
-#include "detail/VolumeBuilder.hh"
+#include "detail/BuildConvexRegion.hh"
 
 #if CELERITAS_USE_JSON
 #    include "ObjectIO.json.hh"
@@ -30,20 +25,8 @@ namespace orangeinp
  */
 NodeId ShapeBase::build(VolumeBuilder& vb) const
 {
-    // Set input attributes for surface state
-    detail::ConvexSurfaceState css;
-    css.transform = &vb.local_transform();
-    css.object_name = this->label();
-    css.make_face_name = {};  // No prefix for standalone shapes
-
-    // Construct surfaces
-    auto sb = ConvexSurfaceBuilder(&vb.unit_builder(), &css);
-    this->interior().build(sb);
-
-    // Intersect the given surfaces to create a new CSG node
-    return vb.insert_region(Label{std::move(css.object_name)},
-                            Joined{op_and, std::move(css.nodes)},
-                            calc_merged_bzone(css));
+    return detail::build_convex_region(
+        vb, std::string{this->label()}, {}, this->interior());
 }
 
 //---------------------------------------------------------------------------//

src/orange/orangeinp/Shape.hh

@@ -58,9 +58,12 @@ class ShapeBase : public ObjectInterface
 /*!
  * Shape that holds a convex region and forwards construction args to it.
  *
- * Construct as:
- *
- *    BoxShape s{"mybox", Real3{1, 2, 3}};
+ * Construct as: \code
+    BoxShape s{"mybox", Real3{1, 2, 3}};
+ * \endcode
+ * or \code
+ *  Shape s{"mybox", Box{{1, 2, 3}}};
+ * \endcode
  *
  * See ConvexRegion.hh for a list of the regions and their construction
  * arguments.

src/orange/orangeinp/detail/BuildConvexRegion.cc

@@ -0,0 +1,52 @@
+//----------------------------------*-C++-*----------------------------------//
+// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file orange/orangeinp/detail/BuildConvexRegion.cc
+//---------------------------------------------------------------------------//
+#include "BuildConvexRegion.hh"
+
+#include "orange/orangeinp/ConvexRegion.hh"
+#include "orange/orangeinp/ConvexSurfaceBuilder.hh"
+#include "orange/surf/FaceNamer.hh"
+
+#include "ConvexSurfaceState.hh"
+#include "CsgUnitBuilder.hh"
+#include "VolumeBuilder.hh"
+
+namespace celeritas
+{
+namespace orangeinp
+{
+namespace detail
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Build a convex region.
+ */
+NodeId build_convex_region(VolumeBuilder& vb,
+                           std::string&& label,
+                           std::string&& face_prefix,
+                           ConvexRegionInterface const& region)
+{
+    // Set input attributes for surface state
+    ConvexSurfaceState css;
+    css.transform = &vb.local_transform();
+    css.object_name = std::move(label);
+    css.make_face_name = FaceNamer{std::move(face_prefix)};
+
+    // Construct surfaces
+    auto sb = ConvexSurfaceBuilder(&vb.unit_builder(), &css);
+    region.build(sb);
+
+    // Intersect the given surfaces to create a new CSG node
+    return vb.insert_region(Label{std::move(css.object_name)},
+                            Joined{op_and, std::move(css.nodes)},
+                            calc_merged_bzone(css));
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace detail
+}  // namespace orangeinp
+}  // namespace celeritas

src/orange/orangeinp/detail/BuildConvexRegion.hh

@@ -0,0 +1,41 @@
+//----------------------------------*-C++-*----------------------------------//
+// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file orange/orangeinp/detail/BuildConvexRegion.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include <string>
+
+#include "orange/orangeinp/CsgTypes.hh"
+
+namespace celeritas
+{
+namespace orangeinp
+{
+class ConvexRegionInterface;
+namespace detail
+{
+class VolumeBuilder;
+//---------------------------------------------------------------------------//
+
+// Build a convex region
+NodeId build_convex_region(VolumeBuilder& vb,
+                           std::string&& label,
+                           std::string&& face_prefix,
+                           ConvexRegionInterface const& region);
+
+//! Build a convex region with no face prefix
+inline NodeId build_convex_region(VolumeBuilder& vb,
+                                  std::string&& label,
+                                  ConvexRegionInterface const& region)
+{
+    return build_convex_region(vb, std::move(label), {}, region);
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace detail
+}  // namespace orangeinp
+}  // namespace celeritas

src/orange/surf/FaceNamer.hh

@@ -32,8 +32,8 @@ class FaceNamer
     // Construct with no prefix
     FaceNamer() = default;
 
-    //! Construct with prefix
-    explicit FaceNamer(std::string&& prefix) : prefix_{std::move(prefix)} {}
+    // Construct with prefix
+    explicit inline FaceNamer(std::string&& prefix);
 
     // Apply to a surface with known type
     template<class S>
@@ -86,6 +86,18 @@ class FaceNamer
     };
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with prefix.
+ */
+FaceNamer::FaceNamer(std::string&& prefix) : prefix_{std::move(prefix)}
+{
+    if (!prefix_.empty() && prefix_.back() != '.')
+    {
+        prefix_.push_back('.');
+    }
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Apply to a surface with known type.

test/corecel/math/Algorithms.test.cc

@@ -340,6 +340,20 @@ TEST(MathTest, diffsq)
 
 //---------------------------------------------------------------------------//
 
+TEST(MathTest, eumod)
+{
+    // Wrap numbers to between [0, 360)
+    EXPECT_DOUBLE_EQ(270.0, eumod(-90.0 - 360.0, 360.0));
+    EXPECT_DOUBLE_EQ(270.0, eumod(-90.0, 360.0));
+    EXPECT_DOUBLE_EQ(0.0, eumod(0.0, 360.0));
+    EXPECT_DOUBLE_EQ(45.0, eumod(45.0, 360.0));
+    EXPECT_DOUBLE_EQ(0.0, eumod(360.0, 360.0));
+    EXPECT_DOUBLE_EQ(15.0, eumod(375.0, 360.0));
+    EXPECT_DOUBLE_EQ(30.0, eumod(720.0 + 30, 360.0));
+}
+
+//---------------------------------------------------------------------------//
+
 TEST(MathTest, sincos)
 {
     {