Add Parallelepiped shape

src/geocel/BoundingBox.hh

@@ -34,7 +34,7 @@ namespace celeritas
  * it (center, surface area, volume) are prohibited.
  *
  * A "degenerate" bounding box is one that is well-defined but has zero volume
- * because at least one lower coorindate is equal to the corresponding upper
+ * because at least one lower coordinate is equal to the corresponding upper
  * coordinate. Any point on the surface of this bounding box is still "inside".
  * It may have nonzero surface area but will have zero volume.
  */

src/orange/orangeinp/ConvexRegion.cc

@@ -354,6 +354,91 @@ void InfWedge::output(JsonPimpl* j) const
     to_json_pimpl(j, *this);
 }
 
+//---------------------------------------------------------------------------//
+// PARALLELEPIPED
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with a 3-vector of half-edges and three angles (in radians).
+ */
+Parallelepiped::Parallelepiped(Real3 const& half_projs,
+                               Turn alpha,
+                               Turn theta,
+                               Turn phi)
+    : hpr_{half_projs}, alpha_{alpha}, theta_{theta}, phi_{phi}
+{
+    for (auto ax : range(Axis::size_))
+    {
+        CELER_VALIDATE(hpr_[to_int(ax)] > 0,
+                       << "nonpositive half-edge - roughly along "
+                       << to_char(ax) << " axis: " << hpr_[to_int(ax)]);
+    }
+
+    CELER_VALIDATE(alpha_ >= zero_quantity() && alpha_ < Turn{0.25},
+                   << "invalid angle " << alpha_.value()
+                   << " [turns]: must be in the range [0, 0.25]");
+
+    CELER_VALIDATE(theta_ >= zero_quantity() && theta_ < Turn{0.25},
+                   << "invalid angle " << theta_.value()
+                   << " [turns]: must be in the range [0, 0.25]");
+
+    CELER_VALIDATE(phi_ >= zero_quantity() && phi_ < Turn{1.},
+                   << "invalid angle " << phi_.value()
+                   << " [turns]: must be in the range [0, 1)");
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Build surfaces.
+ */
+void Parallelepiped::build(ConvexSurfaceBuilder& insert_surface) const
+{
+    constexpr auto X = to_int(Axis::x);
+    constexpr auto Y = to_int(Axis::y);
+    constexpr auto Z = to_int(Axis::z);
+
+    // 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
+    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
+    auto xnorm = make_unit_vector(cross_product(b, c));
+    auto ynorm = make_unit_vector(cross_product(c, a));
+    auto xoffset = dot_product(a, xnorm);
+    auto yoffset = dot_product(b, ynorm);
+
+    // Build top and bottom planes
+    insert_surface(Sense::outside, PlaneZ{-hpr_[Z]});
+    insert_surface(Sense::inside, PlaneZ{hpr_[Z]});
+
+    // Build the side planes roughly perpendicular to y-axis
+    insert_surface(Sense::outside, Plane{ynorm, -yoffset});
+    insert_surface(Sense::inside, Plane{ynorm, yoffset});
+
+    // Build the side planes roughly perpendicular to x-axis
+    insert_surface(Sense::inside, Plane{-xnorm, xoffset});
+    insert_surface(Sense::inside, Plane{xnorm, xoffset});
+
+    // Add an exterior bounding box
+    auto half_diagonal = a + b + c;
+    insert_surface(Sense::inside, BBox{-half_diagonal, half_diagonal});
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Write output to the given JSON object.
+ */
+void Parallelepiped::output(JsonPimpl* j) const
+{
+    to_json_pimpl(j, *this);
+}
+
 //---------------------------------------------------------------------------//
 // PRISM
 //---------------------------------------------------------------------------//

src/orange/orangeinp/ConvexRegion.hh

@@ -230,6 +230,60 @@ class InfWedge final : public ConvexRegionInterface
     Turn interior_;
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * A general parallelepiped centered on the origin.
+ *
+ * A parallelepiped is a shape having 3 pairs of parallel faces out of
+ * which one is parallel with the XY plane (Z faces). All faces are
+ * parallelograms in the general case. The Z faces have 2 edges parallel
+ * with the X-axis.
+ *
+ * The shape has the center in the origin and it is defined by:
+ *
+ *   - `dX, dY, dZ:` half-lengths of the projections of the edges on X, Y
+ *     and Z. The lower Z face is positioned at `-dZ`, while the upper at
+ * `+dZ`.
+ *   - `alpha:` angle between the segment defined by the centers of the
+ *     X-parallel edges and Y axis `[-90,90]` in degrees;
+ *   - `theta:` polar angle of the shape's main axis, e.g. the segment defined
+ *     by the centers of the Z faces;
+ *   - `phi:` azimuthal angle of the shape's main axis (explained above).
+ */
+class Parallelepiped final : public ConvexRegionInterface
+{
+  public:
+    // Construct with half widths and 3 angles
+    Parallelepiped(Real3 const& halfedges, Turn alpha, Turn theta, Turn phi);
+
+    // Build surfaces
+    void build(ConvexSurfaceBuilder&) const final;
+
+    // Output to JSON
+    void output(JsonPimpl*) const final;
+
+    //// ACCESSORS ////
+
+    //! Half-lengths of edge projections along each axis
+    Real3 const& half_projs() 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)
+    Turn theta() const { return theta_; }
+    //! Azimuthal angle of main axis (in turns)
+    Turn phi() const { return phi_; }
+
+  private:
+    // half-lengths
+    Real3 hpr_;
+    // angle between slanted y-edges and the y-axis
+    Turn alpha_;
+    // polar angle of main axis
+    Turn theta_;
+    // azimuthal angle of main axis
+    Turn phi_;
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * A regular, z-extruded polygon centered on the origin.

src/orange/orangeinp/ObjectIO.json.cc

@@ -191,6 +191,14 @@ void to_json(nlohmann::json& j, InfWedge const& cr)
          {"start", cr.start().value()},
          {"interior", cr.interior().value()}};
 }
+void to_json(nlohmann::json& j, Parallelepiped const& cr)
+{
+    j = {{"_type", "parallelepiped"},
+         SIO_ATTR_PAIR(cr, half_projs),
+         {"alpha", cr.alpha().value()},
+         {"theta", cr.theta().value()},
+         {"phi", cr.phi().value()}};
+}
 void to_json(nlohmann::json& j, Prism const& cr)
 {
     j = {{"_type", "prism"},

src/orange/orangeinp/ObjectIO.json.hh

@@ -33,6 +33,7 @@ class Cone;
 class Cylinder;
 class Ellipsoid;
 class InfWedge;
+class Parallelepiped;
 class Prism;
 class Sphere;
 
@@ -59,6 +60,7 @@ 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, Parallelepiped 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

@@ -46,6 +46,7 @@ template class Shape<Box>;
 template class Shape<Cone>;
 template class Shape<Cylinder>;
 template class Shape<Ellipsoid>;
+template class Shape<Parallelepiped>;
 template class Shape<Prism>;
 template class Shape<Sphere>;
 

src/orange/orangeinp/Shape.hh

@@ -30,6 +30,7 @@ namespace orangeinp
  * - \c ConeShape
  * - \c CylinderShape
  * - \c EllipsoidShape
+ * - \c ParallelepipedShape
  * - \c PrismShape
  * - \c SphereShape
  */
@@ -115,6 +116,7 @@ using BoxShape = Shape<Box>;
 using ConeShape = Shape<Cone>;
 using CylinderShape = Shape<Cylinder>;
 using EllipsoidShape = Shape<Ellipsoid>;
+using ParallelepipedShape = Shape<Parallelepiped>;
 using PrismShape = Shape<Prism>;
 using SphereShape = Shape<Sphere>;
 

test/orange/orangeinp/ConvexRegion.test.cc

@@ -474,6 +474,125 @@ TEST_F(InfWedgeTest, half_turn)
     }
 }
 
+//---------------------------------------------------------------------------//
+// PARALLELEPIPED
+//---------------------------------------------------------------------------//
+using ParallelepipedTest = ConvexRegionTest;
+
+TEST_F(ParallelepipedTest, errors)
+{
+    EXPECT_THROW(Parallelepiped({0, 1, 2}, Turn(0.1), Turn(0.1), Turn(0.1)),
+                 RuntimeError);  // bad x
+    EXPECT_THROW(Parallelepiped({2, 0, 1}, Turn(0.2), Turn(0.0), Turn(0.1)),
+                 RuntimeError);  // bad y
+    EXPECT_THROW(Parallelepiped({2, 1, 0}, Turn(0.1), Turn(0.1), Turn(0.1)),
+                 RuntimeError);  // bad z
+
+    Real3 sides{1, 2, 3};
+    EXPECT_THROW(Parallelepiped(sides, Turn(0.3), Turn(0.1), Turn(0.1)),
+                 RuntimeError);  // alpha
+    EXPECT_THROW(Parallelepiped(sides, Turn(0.1), Turn(0.3), Turn(0.1)),
+                 RuntimeError);  // theta
+    EXPECT_THROW(Parallelepiped(sides, Turn(0.1), Turn(0.1), Turn(1.0)),
+                 RuntimeError);  // phi
+}
+
+TEST_F(ParallelepipedTest, box)
+{
+    Real3 sides{1, 2, 3};
+    auto result
+        = this->test(Parallelepiped(sides, Turn(0.0), Turn(0.0), Turn(0.0)));
+
+    static char const expected_node[] = "all(+0, -1, +2, -3, +4, -5)";
+    static char const* const expected_surfaces[] = {"Plane: z=-3",
+                                                    "Plane: z=3",
+                                                    "Plane: y=-2",
+                                                    "Plane: y=2",
+                                                    "Plane: x=-1",
+                                                    "Plane: x=1"};
+
+    EXPECT_EQ(expected_node, result.node);
+    EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
+    EXPECT_VEC_SOFT_EQ((Real3{-1, -2, -3}), result.interior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{1, 2, 3}), result.interior.upper());
+    EXPECT_VEC_SOFT_EQ((Real3{-1, -2, -3}), result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{1, 2, 3}), result.exterior.upper());
+}
+
+TEST_F(ParallelepipedTest, alpha)
+{
+    Real3 sides{1, 2, 3};
+    auto result
+        = this->test(Parallelepiped(sides, Turn(0.1), Turn(0.0), Turn(0.0)));
+
+    static char const expected_node[] = "all(+0, -1, +2, -3, +4, -5)";
+    static char const* const expected_surfaces[]
+        = {"Plane: z=-3",
+           "Plane: z=3",
+           "Plane: y=-1.618",
+           "Plane: y=1.618",
+           "Plane: n={0.80902,-0.58779,0}, d=-0.80902",
+           "Plane: n={0.80902,-0.58779,0}, d=0.80902"};
+
+    EXPECT_EQ(expected_node, result.node);
+    EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
+    EXPECT_FALSE(result.interior) << result.interior;
+    EXPECT_VEC_SOFT_EQ((Real3{-2.1755705045849, -1.6180339887499, -3}),
+                       result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{2.1755705045849, 1.6180339887499, 3}),
+                       result.exterior.upper());
+}
+
+TEST_F(ParallelepipedTest, theta)
+{
+    Real3 sides{1, 2, 3};
+    auto result
+        = this->test(Parallelepiped(sides, Turn(0), Turn(0.1), Turn(0)));
+
+    static char const expected_node[] = "all(+0, -1, +2, -3, +4, -5)";
+    static char const* const expected_surfaces[]
+        = {"Plane: z=-3",
+           "Plane: z=3",
+           "Plane: y=-2",
+           "Plane: y=2",
+           "Plane: n={0.80902,0,-0.58779}, d=-0.80902",
+           "Plane: n={0.80902,0,-0.58779}, d=0.80902"};
+
+    EXPECT_EQ(expected_node, result.node);
+    EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
+    EXPECT_FALSE(result.interior) << result.interior;
+    EXPECT_VEC_SOFT_EQ((Real3{-2.7633557568774, -2, -2.4270509831248}),
+                       result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ((Real3{2.7633557568774, 2, 2.4270509831248}),
+                       result.exterior.upper());
+}
+
+TEST_F(ParallelepipedTest, full)
+{
+    Real3 sides{1, 2, 3};
+    auto result
+        = this->test(Parallelepiped(sides, Turn(0.1), Turn(0.05), Turn(0.15)));
+
+    static char const expected_node[] = "all(+0, -1, +2, -3, +4, -5)";
+    static char const* const expected_surfaces[]
+        = {"Plane: z=-3",
+           "Plane: z=3",
+           "Plane: n={0,0.96714,-0.25423}, d=-1.5649",
+           "Plane: n={0,0.96714,-0.25423}, d=1.5649",
+           "Plane: n={0.80902,-0.58779,0}, d=-0.80902",
+           "Plane: n={0.80902,-0.58779,0}, d=0.80902"};
+
+    EXPECT_EQ(expected_node, result.node);
+    EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
+    EXPECT_FALSE(result.interior) << result.interior;
+    EXPECT_VEC_SOFT_EQ(
+        (Real3{-2.720477400589, -2.3680339887499, -2.8531695488855}),
+        result.exterior.lower());
+    EXPECT_VEC_SOFT_EQ(
+        (Real3{2.720477400589, 2.3680339887499, 2.8531695488855}),
+        result.exterior.upper());
+}
+
 //---------------------------------------------------------------------------//
 // PRISM
 //---------------------------------------------------------------------------//