Skip to content

Commit

Permalink
Add minimal version of GenTrap (Arb8) shape (celeritas-project#1171)
Browse files Browse the repository at this point in the history
  • Loading branch information
mrguilima committed Apr 3, 2024
1 parent 3a31f2e commit 5edf749
Show file tree
Hide file tree
Showing 10 changed files with 527 additions and 0 deletions.
86 changes: 86 additions & 0 deletions src/orange/orangeinp/ConvexRegion.cc
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "corecel/io/JsonPimpl.hh"
#include "geocel/BoundingBox.hh"
#include "geocel/Types.hh"
#include "orange/orangeinp/detail/PolygonUtils.hh"
#include "orange/surf/ConeAligned.hh"
#include "orange/surf/CylCentered.hh"
#include "orange/surf/PlaneAligned.hh"
Expand All @@ -32,6 +33,7 @@ namespace orangeinp
{
namespace
{
using Real2 = Array<real_type, 2>;
//---------------------------------------------------------------------------//
/*!
* Create a z-aligned bounding box infinite along z and symmetric in r.
Expand Down Expand Up @@ -309,6 +311,90 @@ void Ellipsoid::output(JsonPimpl* j) const
to_json_pimpl(j, *this);
}

//---------------------------------------------------------------------------//
// GENTRAP
//---------------------------------------------------------------------------//
/*!
* Construct from half Z height and 1-4 vertices for top and bottom planes.
*/
GenTrap::GenTrap(real_type halfz, VecReal2 const& lo, VecReal2 const& hi)
: hz_{halfz}, lo_{std::move(lo)}, hi_{std::move(hi)}
{
CELER_VALIDATE(hz_ > 0, << "nonpositive halfheight: " << hz_);
CELER_VALIDATE(lo_.size() >= 3 && lo_.size() <= 4,
<< "invalid number of vertices (" << lo_.size()
<< ") for -Z polygon");
CELER_VALIDATE(hi_.size() >= 3 && hi_.size() <= 4,
<< "invalid number of vertices <" << hi_.size()
<< ") for +Z polygon");

CELER_VALIDATE(lo_.size() >= 3 || hi_.size() >= 3,
<< "not enough vertices for both of the +Z/-Z polygons.");

// Input vertices must be arranged in a CCW order
CELER_VALIDATE(detail::is_convex(make_span(lo_)),
<< "-Z polygon is not convex");
CELER_VALIDATE(detail::is_convex(make_span(hi_)),
<< "+Z polygon is not convex");

// TODO: Temporarily ensure that all side faces are planar
for (auto i : range(lo_.size()))
{
auto j = (i + 1) % lo_.size();
Real3 const a{lo_[i][0], lo_[i][1], -hz_};
Real3 const b{lo_[j][0], lo_[j][1], -hz_};
Real3 const c{hi_[j][0], hi_[j][1], hz_};
Real3 const d{hi_[i][0], hi_[i][1], hz_};

// *Temporarily* throws if a side face is not planar
if (!detail::is_planar(a, b, c, d))
{
CELER_NOT_IMPLEMENTED("non-planar side faces");
}
}
}

//---------------------------------------------------------------------------//
/*!
* Build surfaces.
*/
void GenTrap::build(ConvexSurfaceBuilder& insert_surface) const
{
// Build the bottom and top planes
insert_surface(Sense::outside, PlaneZ{-hz_});
insert_surface(Sense::inside, PlaneZ{hz_});

// Build the side planes
for (auto i : range(lo_.size()))
{
auto j = (i + 1) % lo_.size();
Real3 const a{lo_[i][0], lo_[i][1], -hz_};
Real3 const b{lo_[j][0], lo_[j][1], -hz_};
Real3 const c{hi_[j][0], hi_[j][1], hz_};
Real3 const d{hi_[i][0], hi_[i][1], hz_};

// Calculate plane parameters
auto normal = make_unit_vector(cross_product(b - a, c - b));
auto offset = dot_product(d, normal);

// *Temporarily* throws if a side face is not planar
CELER_ASSERT(celeritas::orangeinp::detail::is_planar(a, b, c, d));
CELER_ASSERT(SoftZero<real_type>{}(dot_product(d - a, normal)));

// Insert the plane
insert_surface(Sense::inside, Plane{normal, offset});
}
}

//---------------------------------------------------------------------------//
/*!
* Write output to the given JSON object.
*/
void GenTrap::output(JsonPimpl* j) const
{
to_json_pimpl(j, *this);
}

//---------------------------------------------------------------------------//
// INFWEDGE
//---------------------------------------------------------------------------//
Expand Down
47 changes: 47 additions & 0 deletions src/orange/orangeinp/ConvexRegion.hh
Original file line number Diff line number Diff line change
Expand Up @@ -196,6 +196,53 @@ class Ellipsoid final : public ConvexRegionInterface
Real3 radii_;
};

//---------------------------------------------------------------------------//
/*!
* A generalized trapezoid, inspired by VecGeom's GenTrap and also ROOT's Arb8.
*
* A GenTrap represents a general trapezoidal volume with up to eight vertices,
* or two 4-point sitting on two parallel planes perpendicular to Z axis.
* Both sets of up to four points provided need to be in counter-clockwise
* ordering. This ensures that the -z face will have an outward normal, and
* that the +z face points will correspond to their -z face counterparts for
* each side face.
* TODO: Add a check for this.
* TODO: Add proper treatment for degenerate cases.
*/
class GenTrap final : public ConvexRegionInterface
{
//!@{
//! \name Type aliases
using Real2 = Array<real_type, 2>;
using VecReal2 = std::vector<Real2>;
//!@}

public:
// Construct from half Z height and 1-4 vertices for top and bottom planes
GenTrap(real_type halfz, VecReal2 const& lo, VecReal2 const& hi);

// Build surfaces
void build(ConvexSurfaceBuilder&) const final;

// Output to JSON
void output(JsonPimpl*) const final;

//// ACCESSORS ////

//! Half-height along Z
real_type halfheight() const { return hz_; }

//! Polygon on -z face
VecReal2 const& lower() const { return lo_; }
//! Polygon on +z face
VecReal2 const& upper() const { return hi_; }

private:
real_type hz_; //!< half-height
VecReal2 lo_; //!< corners of the -z face
VecReal2 hi_; //!< corners of the +z face
};

//---------------------------------------------------------------------------//
/*!
* An open wedge shape from the Z axis.
Expand Down
7 changes: 7 additions & 0 deletions src/orange/orangeinp/ObjectIO.json.cc
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,13 @@ void to_json(nlohmann::json& j, Ellipsoid const& cr)
{
j = {{"_type", "ellipsoid"}, SIO_ATTR_PAIR(cr, radii)};
}
void to_json(nlohmann::json& j, GenTrap const& cr)
{
j = {{"_type", "gentrap"},
SIO_ATTR_PAIR(cr, halfheight),
SIO_ATTR_PAIR(cr, lower),
SIO_ATTR_PAIR(cr, upper)};
}
void to_json(nlohmann::json& j, InfWedge const& cr)
{
j = {{"_type", "infwedge"},
Expand Down
2 changes: 2 additions & 0 deletions src/orange/orangeinp/ObjectIO.json.hh
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ class Box;
class Cone;
class Cylinder;
class Ellipsoid;
class GenTrap;
class InfWedge;
class Parallelepiped;
class Prism;
Expand Down Expand Up @@ -59,6 +60,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, GenTrap 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);
Expand Down
1 change: 1 addition & 0 deletions src/orange/orangeinp/Shape.cc
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ template class Shape<Box>;
template class Shape<Cone>;
template class Shape<Cylinder>;
template class Shape<Ellipsoid>;
template class Shape<GenTrap>;
template class Shape<Parallelepiped>;
template class Shape<Prism>;
template class Shape<Sphere>;
Expand Down
2 changes: 2 additions & 0 deletions src/orange/orangeinp/Shape.hh
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ namespace orangeinp
* - \c ConeShape
* - \c CylinderShape
* - \c EllipsoidShape
* - \c GenTrapShape
* - \c ParallelepipedShape
* - \c PrismShape
* - \c SphereShape
Expand Down Expand Up @@ -116,6 +117,7 @@ using BoxShape = Shape<Box>;
using ConeShape = Shape<Cone>;
using CylinderShape = Shape<Cylinder>;
using EllipsoidShape = Shape<Ellipsoid>;
using GenTrapShape = Shape<GenTrap>;
using ParallelepipedShape = Shape<Parallelepiped>;
using PrismShape = Shape<Prism>;
using SphereShape = Shape<Sphere>;
Expand Down
111 changes: 111 additions & 0 deletions src/orange/orangeinp/detail/PolygonUtils.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
//----------------------------------*-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/PolygonUtils.hh
//! \brief Utility standalone functions for polygons in 2D or 3D space.
//---------------------------------------------------------------------------//
#pragma once

#include <iostream>
#include <vector>

#include "corecel/cont/Array.hh"
#include "corecel/cont/Range.hh"
#include "corecel/cont/Span.hh"
#include "corecel/math/ArrayOperators.hh"
#include "corecel/math/ArrayUtils.hh"
#include "orange/OrangeTypes.hh"

namespace celeritas
{
namespace orangeinp
{
namespace detail
{
//---------------------------------------------------------------------------//
using Real2 = Array<real_type, 2>;

//---------------------------------------------------------------------------//
/*!
* Polygon orientation based on ordering of vertices.
*/
enum class Orientation
{
clockwise = -1,
collinear = 0,
counterclockwise = 1,
};

//---------------------------------------------------------------------------//
// FREE FUNCTIONS
//---------------------------------------------------------------------------//
/*!
* Find orientation of ordered vertices in 2D coordinates.
*/
CELER_FORCEINLINE_FUNCTION Orientation orientation(Real2 const& a,
Real2 const& b,
Real2 const& c)
{
auto crossp = (b[0] - a[0]) * (c[1] - b[1]) - (b[1] - a[1]) * (c[0] - b[0]);
return crossp < 0 ? Orientation::clockwise
: crossp > 0 ? Orientation::counterclockwise
: Orientation::collinear;
}

//---------------------------------------------------------------------------//
/*!
* Check if a 2D polygon is convex.
*
* \param corners the vertices of the polygon
* \param degen_ok allow consecutive degenerate points
*/
CELER_FUNCTION bool
is_convex(Span<const Real2> const& corners, bool degen_ok = false)
{
CELER_EXPECT(!corners.empty());
// The cross product of all vector pairs corresponding to ordered
// consecutive segments has to be positive.
auto crossp = [&](Real2 const& v1, Real2 const& v2) {
return v1[0] * v2[1] - v1[1] * v2[0];
};

// Use the cross_product(last, first) as sign reference
auto num_corners = corners.size();
Real2 vec1 = corners[0] - corners[num_corners - 1];
Real2 vec2 = corners[1] - corners[0];
real_type const ref = crossp(vec1, vec2);
for (auto i : range(num_corners - 1))
{
vec1 = vec2;
vec2 = corners[(i + 2) % num_corners] - corners[(i + 1) % num_corners];
auto val = crossp(vec1, vec2);
// Make sure the sign is the same as the reference
if (val * ref < 0 || (!degen_ok && val == 0))
{
return false;
}
}
return true;
}

//---------------------------------------------------------------------------//
/*!
* Check for coplanarity of four 3D polygon vertices.
*/
CELER_FUNCTION bool
is_planar(Real3 const& a, Real3 const& b, Real3 const& c, Real3 const& d)
{
// Use the cross_product(last, first) as sign reference
auto norm = make_unit_vector(cross_product(b - a, c - a));
auto val = dot_product(norm, d - a);

// FIXME: SoftEqual and SoftZero should have rel = abs
return SoftZero{SoftEqual<>{}.rel()}(val);
}

//---------------------------------------------------------------------------//
} // namespace detail
} // namespace orangeinp
} // namespace celeritas
1 change: 1 addition & 0 deletions test/orange/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@ celeritas_add_test(orangeinp/Transformed.test.cc)
celeritas_add_test(orangeinp/UnitProto.test.cc)

celeritas_add_test(orangeinp/detail/BoundingZone.test.cc)
celeritas_add_test(orangeinp/detail/PolygonUtils.test.cc)
celeritas_add_test(orangeinp/detail/ProtoMap.test.cc)
celeritas_add_test(orangeinp/detail/LocalSurfaceInserter.test.cc)
celeritas_add_test(orangeinp/detail/TransformInserter.test.cc)
Expand Down

0 comments on commit 5edf749

Please sign in to comment.