diff --git a/src/geocel/BoundingBox.hh b/src/geocel/BoundingBox.hh index 2b0bbfb397..5caf294135 100644 --- a/src/geocel/BoundingBox.hh +++ b/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. */ diff --git a/src/orange/orangeinp/ConvexRegion.cc b/src/orange/orangeinp/ConvexRegion.cc index 90ce78ecc9..8b24f0f95f 100644 --- a/src/orange/orangeinp/ConvexRegion.cc +++ b/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. + */ +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_ > -Turn{0.25} && alpha_ < Turn{0.25}, + << "invalid angle " << alpha_.value() + << " [turns]: must be in the range (-0.25, 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 //---------------------------------------------------------------------------// diff --git a/src/orange/orangeinp/ConvexRegion.hh b/src/orange/orangeinp/ConvexRegion.hh index 392738ad13..7824280b2d 100644 --- a/src/orange/orangeinp/ConvexRegion.hh +++ b/src/orange/orangeinp/ConvexRegion.hh @@ -230,6 +230,62 @@ 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. Note that all angle parameters are expressed in terms + * of fractions of a 360deg turn. + * + * The shape has the center in the origin and it is defined by: + * + * - `halfedges:` a 3-vector (dY, dY, dZ) with half-lengths of the + * projections of the edges on X, Y, Z. The lower Z face is positioned at + * `-dZ`, and the upper one at `+dZ`. + * - `alpha:` angle between the segment defined by the centers of the + * X-parallel edges and Y axis. Validity range is `(-1/4, 1/4)`; + * - `theta:` polar angle of the shape's main axis, e.g. the segment defined + * by the centers of the Z faces. Validity range is `[0, 1/4)`; + * - `phi:` azimuthal angle of the shape's main axis (as explained above). + *   Validity range is `[0, 1)`. + */ +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. diff --git a/src/orange/orangeinp/ObjectIO.json.cc b/src/orange/orangeinp/ObjectIO.json.cc index 93526f3690..36a1a876a6 100644 --- a/src/orange/orangeinp/ObjectIO.json.cc +++ b/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"}, diff --git a/src/orange/orangeinp/ObjectIO.json.hh b/src/orange/orangeinp/ObjectIO.json.hh index c6ef5a4620..4f846bcd90 100644 --- a/src/orange/orangeinp/ObjectIO.json.hh +++ b/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); diff --git a/src/orange/orangeinp/Shape.cc b/src/orange/orangeinp/Shape.cc index 33f26e4794..ac14bcb213 100644 --- a/src/orange/orangeinp/Shape.cc +++ b/src/orange/orangeinp/Shape.cc @@ -46,6 +46,7 @@ template class Shape; template class Shape; template class Shape; template class Shape; +template class Shape; template class Shape; template class Shape; diff --git a/src/orange/orangeinp/Shape.hh b/src/orange/orangeinp/Shape.hh index d48a61995d..32b49d5497 100644 --- a/src/orange/orangeinp/Shape.hh +++ b/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; using ConeShape = Shape; using CylinderShape = Shape; using EllipsoidShape = Shape; +using ParallelepipedShape = Shape; using PrismShape = Shape; using SphereShape = Shape; diff --git a/test/orange/orangeinp/ConvexRegion.test.cc b/test/orange/orangeinp/ConvexRegion.test.cc index a992c84a96..ffa207c14e 100644 --- a/test/orange/orangeinp/ConvexRegion.test.cc +++ b/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 //---------------------------------------------------------------------------//