Add additional ObjectInterface and orangeinp helpers (#1174)

* Add helper function for building solid or shape
* Add transform helper functions
* Update parallelepiped names for consistency
* Clean region store when destroying geometry
* Add factory helper function to Transformed
* Add factory helper to CsgObject

src/geocel/GeantGeoUtils.cc

@@ -17,6 +17,7 @@
 #include <G4LogicalVolumeStore.hh>
 #include <G4PhysicalVolumeStore.hh>
 #include <G4ReflectionFactory.hh>
+#include <G4RegionStore.hh>
 #include <G4SolidStore.hh>
 #include <G4Threading.hh>
 #include <G4TouchableHistory.hh>
@@ -185,6 +186,7 @@ void reset_geant_geometry()
 
         G4PhysicalVolumeStore::Clean();
         G4LogicalVolumeStore::Clean();
+        G4RegionStore::Clean();
         G4SolidStore::Clean();
 #if G4VERSION_NUMBER >= 1100
         G4ReflectionFactory::Instance()->Clean();

src/orange/MatrixUtils.cc

@@ -240,6 +240,27 @@ Mat3 make_rotation(Axis ax, Turn theta, Mat3 const& other)
     return gemm(make_rotation(ax, theta), other);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Construct a transposed matrix.
+ *
+ * This should only be used for preprocessing. Prefer methods that transpose on
+ * the fly.
+ */
+SquareMatrixReal3 make_transpose(SquareMatrixReal3 const& mat)
+{
+    SquareMatrixReal3 result;
+    for (size_type i = 0; i != 3; ++i)
+    {
+        for (size_type j = 0; j != 3; ++j)
+        {
+            result[i][j] = mat[j][i];
+        }
+    }
+
+    return result;
+}
+
 //---------------------------------------------------------------------------//
 // EXPLICIT INSTANTIATIONS
 //---------------------------------------------------------------------------//

src/orange/MatrixUtils.hh

@@ -102,6 +102,9 @@ SquareMatrixReal3 make_rotation(Axis ax, Turn rev);
 // Apply a rotation to an existing C-ordered rotation matrix
 SquareMatrixReal3 make_rotation(Axis ax, Turn rev, SquareMatrixReal3 const&);
 
+// Construct a transposed matrix
+SquareMatrixReal3 make_transpose(SquareMatrixReal3 const&);
+
 //---------------------------------------------------------------------------//
 // INLINE DEFINITIONS
 //---------------------------------------------------------------------------//

src/orange/orangeinp/ConvexRegion.hh

@@ -314,7 +314,7 @@ class Parallelepiped final : public ConvexRegionInterface
     //// ACCESSORS ////
 
     //! Half-lengths of edge projections along each axis
-    Real3 const& half_projs() const { return hpr_; }
+    Real3 const& halfedges() const { return hpr_; }
     //! Angle between slanted y-edges and the y-axis (in turns)
     Turn alpha() const { return alpha_; }
     //! Polar angle of main axis (in turns)

src/orange/orangeinp/CsgObject.cc

@@ -70,6 +70,28 @@ void NegatedObject::output(JsonPimpl* j) const
 template<OperatorToken Op>
 constexpr OperatorToken JoinObjects<Op>::op_token;
 
+//---------------------------------------------------------------------------//
+/*!
+ * Construct a joined object if nontrivial, or return the original.
+ *
+ * This should only be called if the label of the resulting object is not
+ * important.
+ */
+template<OperatorToken Op>
+auto JoinObjects<Op>::or_object(std::string&& label, VecObject&& objects)
+    -> SPConstObject
+{
+    CELER_EXPECT(!objects.empty());
+    if (objects.size() == 1)
+    {
+        SPConstObject result = std::move(objects.front());
+        objects.clear();
+        return result;
+    }
+    return std::make_shared<JoinObjects<Op>>(std::move(label),
+                                             std::move(objects));
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Construct with a name and a vector of objects.

src/orange/orangeinp/CsgObject.hh

@@ -66,6 +66,9 @@ class JoinObjects : public ObjectInterface
     static constexpr OperatorToken op_token = Op;
 
   public:
+    // Construct a joined object if nontrivial, or return the original
+    static SPConstObject or_object(std::string&& label, VecObject&& objects);
+
     // Construct with a label and vector of objects
     JoinObjects(std::string&& label, VecObject&& objects);
 

src/orange/orangeinp/ObjectIO.json.cc

@@ -201,7 +201,7 @@ void to_json(nlohmann::json& j, InfWedge const& cr)
 void to_json(nlohmann::json& j, Parallelepiped const& cr)
 {
     j = {{"_type", "parallelepiped"},
-         SIO_ATTR_PAIR(cr, half_projs),
+         SIO_ATTR_PAIR(cr, halfedges),
          {"alpha", cr.alpha().value()},
          {"theta", cr.theta().value()},
          {"phi", cr.phi().value()}};

src/orange/orangeinp/Solid.cc

@@ -22,6 +22,8 @@
 #    include "ObjectIO.json.hh"
 #endif
 
+#include "Shape.hh"
+
 namespace celeritas
 {
 namespace orangeinp
@@ -111,43 +113,52 @@ void SolidBase::output(JsonPimpl* j) const
 
 //---------------------------------------------------------------------------//
 /*!
- * Construct with an excluded interior and enclosed angle.
+ * Return a solid or shape given an optional interior or enclosed angle.
  */
 template<class T>
-Solid<T>::Solid(std::string&& label,
-                T&& interior,
-                T&& excluded,
-                SolidEnclosedAngle enclosed)
-    : label_{std::move(label)}
-    , interior_{std::move(interior)}
-    , exclusion_{std::move(excluded)}
-    , enclosed_(enclosed)
+auto Solid<T>::or_shape(std::string&& label,
+                        T&& interior,
+                        OptionalRegion&& excluded,
+                        SolidEnclosedAngle&& enclosed) -> SPConstObject
 {
-    CELER_VALIDATE(interior_.encloses(*exclusion_),
-                   << "solid '" << this->label()
-                   << "' was given an interior region that is not enclosed by "
-                      "its exterior");
+    if (!excluded && !enclosed)
+    {
+        // Just a shape
+        return std::make_shared<Shape<T>>(std::move(label),
+                                          std::move(interior));
+    }
+
+    return std::make_shared<Solid<T>>(std::move(label),
+                                      std::move(interior),
+                                      std::move(excluded),
+                                      std::move(enclosed));
 }
 
 //---------------------------------------------------------------------------//
 /*!
- * Construct with an enclosed angle.
+ * Construct with optional components.
  */
 template<class T>
-Solid<T>::Solid(std::string&& label, T&& interior, SolidEnclosedAngle enclosed)
+Solid<T>::Solid(std::string&& label,
+                T&& interior,
+                OptionalRegion&& excluded,
+                SolidEnclosedAngle&& enclosed)
     : label_{std::move(label)}
     , interior_{std::move(interior)}
-    , enclosed_(enclosed)
+    , exclusion_{std::move(excluded)}
+    , enclosed_{std::move(enclosed)}
 {
-    CELER_VALIDATE(enclosed_,
+    CELER_VALIDATE(exclusion_ || enclosed,
+                   << "solid requires either an excluded region or a shape");
+    CELER_VALIDATE(!exclusion_ || interior_.encloses(*exclusion_),
                    << "solid '" << this->label()
-                   << "' did not exclude an interior or a wedge (use a Shape "
-                      "instead)");
+                   << "' was given an interior region that is not enclosed by "
+                      "its exterior");
 }
 
 //---------------------------------------------------------------------------//
 /*!
- * Construct with only an excluded interior.
+ * Construct with an excluded interior.
  */
 template<class T>
 Solid<T>::Solid(std::string&& label, T&& interior, T&& excluded)
@@ -158,6 +169,21 @@ Solid<T>::Solid(std::string&& label, T&& interior, T&& excluded)
 {
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with an enclosed angle.
+ */
+template<class T>
+Solid<T>::Solid(std::string&& label, T&& interior, SolidEnclosedAngle&& enclosed)
+    : Solid{
+        std::move(label), std::move(interior), std::nullopt, std::move(enclosed)}
+{
+    CELER_VALIDATE(enclosed_,
+                   << "solid '" << this->label()
+                   << "' did not exclude an interior or a wedge (use a Shape "
+                      "instead)");
+}
+
 //---------------------------------------------------------------------------//
 // EXPLICIT INSTANTIATION
 //---------------------------------------------------------------------------//

src/orange/orangeinp/Solid.hh

@@ -125,14 +125,26 @@ class Solid final : public SolidBase
     static_assert(std::is_base_of_v<ConvexRegionInterface, T>);
 
   public:
-    // Construct with an excluded interior and enclosed angle
+    //!@{
+    //! \name Type aliases
+    using OptionalRegion = std::optional<T>;
+    //!@}
+
+  public:
+    // Return a solid *or* shape given an optional interior or enclosed angle
+    static SPConstObject or_shape(std::string&& label,
+                                  T&& interior,
+                                  OptionalRegion&& excluded,
+                                  SolidEnclosedAngle&& enclosed);
+
+    // Construct with an excluded interior *and/or* enclosed angle
     Solid(std::string&& label,
           T&& interior,
-          T&& excluded,
-          SolidEnclosedAngle enclosed);
+          OptionalRegion&& excluded,
+          SolidEnclosedAngle&& enclosed);
 
-    // Construct with an less-than-full enclosed angle
-    Solid(std::string&& label, T&& interior, SolidEnclosedAngle enclosed);
+    // Construct with only an enclosed angle
+    Solid(std::string&& label, T&& interior, SolidEnclosedAngle&& enclosed);
 
     // Construct with only an excluded interior
     Solid(std::string&& label, T&& interior, T&& excluded);
@@ -152,7 +164,7 @@ class Solid final : public SolidBase
   private:
     std::string label_;
     T interior_;
-    std::optional<T> exclusion_;
+    OptionalRegion exclusion_;
     SolidEnclosedAngle enclosed_;
 };
 

src/orange/orangeinp/Transformed.cc

@@ -20,13 +20,32 @@ namespace celeritas
 {
 namespace orangeinp
 {
+//---------------------------------------------------------------------------//
+/*!
+ * Construct a transformed object if nontrivial, or return the original.
+ */
+auto Transformed::or_object(SPConstObject obj,
+                            VariantTransform const& transform) -> SPConstObject
+{
+    if (std::holds_alternative<NoTransformation>(transform))
+    {
+        return obj;
+    }
+    return std::make_shared<Transformed>(std::move(obj), transform);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Construct with daughter object and transform.
+ *
+ * The input transform should *not* be "no transform". If you don't know
+ * whether it is or not, use the \c Transformed::or_object factory function.
  */
 Transformed::Transformed(SPConstObject obj, VariantTransform const& transform)
     : obj_{std::move(obj)}, transform_{transform}
 {
+    CELER_EXPECT(!transform_.valueless_by_exception());
+    CELER_EXPECT(!std::holds_alternative<NoTransformation>(transform_));
     CELER_EXPECT(obj_);
 }
 

src/orange/orangeinp/Transformed.hh

@@ -22,6 +22,10 @@ namespace orangeinp
 class Transformed final : public ObjectInterface
 {
   public:
+    // Construct a transformed object if nontrivial, or return the original
+    static SPConstObject
+    or_object(SPConstObject obj, VariantTransform const& transform);
+
     // Construct with daughter object and transform
     Transformed(SPConstObject obj, VariantTransform const& transform);
 

src/orange/orangeinp/UnitProto.cc

@@ -340,11 +340,11 @@ auto UnitProto::build(Tol const& tol, ExteriorBoundary ext) const -> Unit
 /*!
  * Construct the daughter's shape in this unit's reference frame.
  */
-std::shared_ptr<Transformed> UnitProto::DaughterInput::make_interior() const
+auto UnitProto::DaughterInput::make_interior() const -> SPConstObject
 {
     CELER_EXPECT(*this);
-    return std::make_shared<Transformed>(this->fill->interior(),
-                                         VariantTransform{this->transform});
+    return Transformed::or_object(this->fill->interior(),
+                                  VariantTransform{this->transform});
 }
 
 //---------------------------------------------------------------------------//

src/orange/orangeinp/UnitProto.hh

@@ -98,7 +98,7 @@ class UnitProto : public ProtoInterface
         explicit inline operator bool() const;
 
         // Construct the daughter's shape in this unit's reference frame
-        std::shared_ptr<Transformed> make_interior() const;
+        SPConstObject make_interior() const;
     };
 
     //! Boundary conditions for the unit

src/orange/transform/Transformation.cc

@@ -16,6 +16,20 @@
 
 namespace celeritas
 {
+//---------------------------------------------------------------------------//
+/*!
+ * Construct by inverting another transformation.
+ */
+Transformation Transformation::from_inverse(Mat3 const& rot, Real3 const& trans)
+{
+    // Transpose the rotation
+    Mat3 const rinv = make_transpose(rot);
+
+    // Calculate the updated position
+    Real3 tinv = gemv(real_type{-1}, rinv, trans, real_type{0}, {});
+    return Transformation{rinv, tinv};
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Construct as an identity transform.

src/orange/transform/Transformation.hh

@@ -63,6 +63,9 @@ class Transformation
     }
 
   public:
+    // Construct by inverting a parent-to-daughter transformation
+    static Transformation from_inverse(Mat3 const& rot, Real3 const& trans);
+
     //// CONSTRUCTORS ////
 
     // Construct and check the input

test/orange/orangeinp/Solid.test.cc

@@ -7,6 +7,8 @@
 //---------------------------------------------------------------------------//
 #include "orange/orangeinp/Solid.hh"
 
+#include "orange/orangeinp/Shape.hh"
+
 #include "CsgTestUtils.hh"
 #include "ObjectTestBase.hh"
 #include "celeritas_test.hh"
@@ -75,6 +77,10 @@ TEST_F(SolidTest, errors)
     // No exclusion, no wedge
     EXPECT_THROW(ConeSolid("cone", Cone{{1, 2}, 10.0}, SolidEnclosedAngle{}),
                  RuntimeError);
+    EXPECT_THROW(
+        ConeSolid(
+            "cone", Cone{{1, 2}, 10.0}, std::nullopt, SolidEnclosedAngle{}),
+        RuntimeError);
 }
 
 TEST_F(SolidTest, inner)
@@ -292,6 +298,24 @@ TEST_F(SolidTest, cyl)
     EXPECT_VEC_EQ(expected_volume_strings, volume_strings(u));
 }
 
+TEST_F(SolidTest, or_shape)
+{
+    {
+        auto shape = ConeSolid::or_shape(
+            "cone", Cone{{1, 2}, 10.0}, std::nullopt, SolidEnclosedAngle{});
+        EXPECT_TRUE(shape);
+        EXPECT_TRUE(dynamic_cast<ConeShape const*>(shape.get()));
+    }
+    {
+        auto solid = ConeSolid::or_shape("cone",
+                                         Cone{{1.1, 2}, 10.0},
+                                         Cone{{0.9, 1.9}, 10.0},
+                                         SolidEnclosedAngle{});
+        EXPECT_TRUE(solid);
+        EXPECT_TRUE(dynamic_cast<ConeSolid const*>(solid.get()));
+    }
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace test
 }  // namespace orangeinp