Skip to content

Commit

Permalink
Fix gentrap orientation (celeritas-project#1206)
Browse files Browse the repository at this point in the history 
* Add translated test

* Add has_orientation function

* Add validation for orientation

* Avoid debug failures when interior and exterior bboxes are both null

* Fix polygon orientation
  • Loading branch information
@sethrj
sethrj committed Apr 26, 2024
1 parent be754c7 commit d6720d6
Show file tree
Hide file tree
Showing 5 changed files with 150 additions and 39 deletions.
58 changes: 37 additions & 21 deletions src/orange/orangeinp/ConvexRegion.cc
View file
Expand Up @@ -12,6 +12,7 @@
#include "corecel/Constants.hh"
#include "corecel/cont/Range.hh"
#include "corecel/io/JsonPimpl.hh"
#include "corecel/math/SoftEqual.hh"
#include "geocel/BoundingBox.hh"
#include "geocel/Types.hh"
#include "orange/orangeinp/detail/PolygonUtils.hh"
Expand Down Expand Up @@ -323,19 +324,32 @@ GenTrap::GenTrap(real_type halfz, VecReal2 const& lo, VecReal2 const& 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");
<< ") for -z polygon");
CELER_VALIDATE(hi_.size() == lo_.size(),
<< "incompatible number of vertices (" << hi_.size()
<< ") for +z polygon: expected " << lo_.size());

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

// Input vertices must be arranged in a CCW order
// Input vertices must be arranged in the same counter/clockwise order
// and be convex
using detail::calc_orientation;
constexpr auto cw = detail::Orientation::clockwise;
CELER_VALIDATE(detail::is_convex(make_span(lo_)),
<< "-Z polygon is not convex");
<< "-z polygon is not convex");
CELER_VALIDATE(detail::is_convex(make_span(hi_)),
<< "+Z polygon is not convex");
<< "+z polygon is not convex");
CELER_VALIDATE(calc_orientation(lo_[0], lo_[1], lo_[2])
== calc_orientation(hi_[0], hi_[1], hi_[2]),
<< "-z and +z polygons have different orientations");
if (calc_orientation(lo_[0], lo_[1], lo_[2]) == cw)
{
// Reverse point orders so it's counterclockwise, needed for vectors to
// point outward
std::reverse(lo_.begin(), lo_.end());
std::reverse(hi_.begin(), hi_.end());
}

// TODO: Temporarily ensure that all side faces are planar
for (auto i : range(lo_.size()))
Expand Down Expand Up @@ -367,22 +381,24 @@ void GenTrap::build(ConvexSurfaceBuilder& insert_surface) const
// Build the side planes
for (auto i : range(lo_.size()))
{
// Viewed from the outside of the trapezoid, the points on the polygon
// here are from the lower left counterclockwise to the upper right
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)));
Real3 const ilo{lo_[i][0], lo_[i][1], -hz_};
Real3 const jlo{lo_[j][0], lo_[j][1], -hz_};
Real3 const jhi{hi_[j][0], hi_[j][1], hz_};
Real3 const ihi{hi_[i][0], hi_[i][1], hz_};

// Calculate outward normal by taking the cross product of the edges
auto normal = make_unit_vector(cross_product(jlo - ilo, ihi - ilo));
// Assert that the surface is (for now!) not twisted
CELER_ASSERT(soft_equal(
dot_product(make_unit_vector(cross_product(ihi - jhi, jlo - jhi)),
normal),
real_type{1}));

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

Expand Down
4 changes: 3 additions & 1 deletion src/orange/orangeinp/detail/ConvexSurfaceState.cc
View file
Expand Up @@ -25,7 +25,9 @@ namespace detail
BoundingZone calc_merged_bzone(ConvexSurfaceState const& css)
{
CELER_EXPECT(css.transform);
CELER_EXPECT(encloses(css.local_bzone.exterior, css.local_bzone.interior));
CELER_EXPECT(
(!css.local_bzone.exterior && !css.local_bzone.interior)
|| encloses(css.local_bzone.exterior, css.local_bzone.interior));
CELER_EXPECT(!css.local_bzone.negated);
CELER_EXPECT(!css.global_bzone.negated);

Expand Down
23 changes: 22 additions & 1 deletion src/orange/orangeinp/detail/PolygonUtils.hh
View file
Expand Up @@ -44,7 +44,7 @@ enum class Orientation
/*!
* Find orientation of ordered vertices in 2D coordinates.
*/
CELER_FORCEINLINE_FUNCTION Orientation orientation(Real2 const& a,
inline CELER_FUNCTION Orientation calc_orientation(Real2 const& a,
Real2 const& b,
Real2 const& c)
{
Expand All @@ -54,6 +54,26 @@ CELER_FORCEINLINE_FUNCTION Orientation orientation(Real2 const& a,
: Orientation::collinear;
}

//---------------------------------------------------------------------------//
/*!
* Test whether a 2D polygon has the given orientation.
*
* The list of input corners must have at least 3 points to be a polygon.
*/
inline CELER_FUNCTION bool
has_orientation(Span<Real2 const> const& corners, Orientation o)
{
CELER_EXPECT(corners.size() > 2);
for (auto i : range(corners.size()))
{
auto j = (i + 1) % corners.size();
auto k = (i + 2) % corners.size();
if (calc_orientation(corners[i], corners[j], corners[k]) != o)
return false;
}
return true;
}

//---------------------------------------------------------------------------//
/*!
* Check if a 2D polygon is convex.
Expand All @@ -65,6 +85,7 @@ 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) {
Expand Down
82 changes: 70 additions & 12 deletions test/orange/orangeinp/ConvexRegion.test.cc
View file
Expand Up @@ -83,8 +83,16 @@ auto ConvexRegionTest::test(ConvexRegionInterface const& r,

ConvexSurfaceBuilder insert_surface{&unit_builder_, &css};
r.build(insert_surface);
EXPECT_TRUE(encloses(css.local_bzone.exterior, css.local_bzone.interior));
EXPECT_TRUE(encloses(css.global_bzone.exterior, css.global_bzone.interior));
if (css.local_bzone.exterior || css.local_bzone.interior)
{
EXPECT_TRUE(
encloses(css.local_bzone.exterior, css.local_bzone.interior));
}
if (css.global_bzone.exterior || css.global_bzone.interior)
{
EXPECT_TRUE(
encloses(css.global_bzone.exterior, css.global_bzone.interior));
}

// Intersect the given surfaces
NodeId node_id
Expand Down Expand Up @@ -495,22 +503,72 @@ TEST_F(GenTrapTest, triang_prism)
auto result = this->test(
GenTrap(3, {{-1, -1}, {-1, 1}, {2, 0}}, {{-1, -1}, {-1, 1}, {2, 0}}));

static char const expected_node[] = "all(+0, -1, -2, +3, +4)";
static char const* const expected_surfaces[]
= {"Plane: z=-3",
"Plane: z=3",
"Plane: x=-1",
"Plane: n={0.31623,0.94868,0}, d=0.63246",
"Plane: n={0.31623,-0.94868,0}, d=0.63246"};
static char const expected_node[] = "all(+0, -1, -2, +3, -4)";
static char const* const expected_surfaces[] = {
"Plane: z=-3",
"Plane: z=3",
"Plane: n={0.31623,0.94868,-0}, d=0.63246",
"Plane: x=-1",
"Plane: n={0.31623,-0.94868,0}, d=0.63246",
};

EXPECT_EQ(expected_node, result.node);
EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
EXPECT_FALSE(result.interior) << result.interior;
EXPECT_VEC_SOFT_EQ((Real3{-inf, -inf, -3}), result.exterior.lower());
EXPECT_VEC_SOFT_EQ((Real3{-1, inf, 3}), result.exterior.upper());
EXPECT_VEC_SOFT_EQ((Real3{-1, -inf, -3}), result.exterior.lower());
EXPECT_VEC_SOFT_EQ((Real3{inf, inf, 3}), result.exterior.upper());
}

TEST_F(GenTrapTest, trapezoid)
{
auto result
= this->test(GenTrap(40,
{{-19, -30}, {-19, 30}, {21, 30}, {21, -30}},
{{-21, -30}, {-21, 30}, {19, 30}, {19, -30}}));

static char const expected_node[] = "all(+0, -1, -2, -3, +4, +5)";
static char const* const expected_surfaces[] = {
"Plane: z=-40",
"Plane: z=40",
"Plane: n={0.99969,-0,0.024992}, d=19.994",
"Plane: y=30",
"Plane: n={0.99969,0,0.024992}, d=-19.994",
"Plane: y=-30",
};

EXPECT_EQ(expected_node, result.node);
EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
EXPECT_FALSE(result.interior) << result.interior;
EXPECT_VEC_SOFT_EQ((Real3{-inf, -30, -40}), result.exterior.lower());
EXPECT_VEC_SOFT_EQ((Real3{inf, 30, 40}), result.exterior.upper());
}

TEST_F(GenTrapTest, trapezoid_trans)
{
// trapezoid but translated -30, -30
auto result
= this->test(GenTrap(40,
{{-49, -60}, {-49, 0}, {-9, 0}, {-9, -60}},
{{-51, -60}, {-51, 0}, {-11, 0}, {-11, -60}}));

static char const expected_node[] = "all(+0, -1, -2, -3, +4, +5)";
static char const* const expected_surfaces[] = {
"Plane: z=-40",
"Plane: z=40",
"Plane: n={0.99969,-0,0.024992}, d=-9.9969",
"Plane: y=0",
"Plane: n={0.99969,0,0.024992}, d=-49.984",
"Plane: y=-60",
};

EXPECT_EQ(expected_node, result.node);
EXPECT_VEC_EQ(expected_surfaces, result.surfaces);
EXPECT_FALSE(result.interior) << result.interior;
EXPECT_VEC_SOFT_EQ((Real3{-inf, -60, -40}), result.exterior.lower());
EXPECT_VEC_SOFT_EQ((Real3{inf, 0, 40}), result.exterior.upper());
}

TEST_F(GenTrapTest, CCWtrap)
TEST_F(GenTrapTest, trapezoid_ccw)
{
auto result
= this->test(GenTrap(40,
Expand Down
22 changes: 18 additions & 4 deletions test/orange/orangeinp/detail/PolygonUtils.test.cc
View file
Expand Up @@ -33,16 +33,30 @@ using constants::pi;
// TESTS
//---------------------------------------------------------------------------//

TEST(PolygonUtilsTest, orientation)
TEST(PolygonUtilsTest, calc_orientation)
{
EXPECT_TRUE(orientation(Real2{0, 0}, Real2{4, 4}, Real2{1, 2})
EXPECT_TRUE(calc_orientation(Real2{0, 0}, Real2{4, 4}, Real2{1, 2})
== Orientation::counterclockwise);
EXPECT_TRUE(orientation(Real2{0, 0}, Real2{4, 4}, Real2{2, 1})
EXPECT_TRUE(calc_orientation(Real2{0, 0}, Real2{4, 4}, Real2{2, 1})
== Orientation::clockwise);
EXPECT_TRUE(orientation(Real2{0, 0}, Real2{4, 4}, Real2{2, 2})
EXPECT_TRUE(calc_orientation(Real2{0, 0}, Real2{4, 4}, Real2{2, 2})
== Orientation::collinear);
}

TEST(PolygonUtilsTest, has_orientation)
{
EXPECT_TRUE(has_orientation(
make_span(VecReal2{{-19, -30}, {-19, 30}, {21, 30}, {21, -30}}),
Orientation::clockwise));
EXPECT_FALSE(has_orientation(
make_span(VecReal2{{-19, -30}, {-19, 30}, {21, 30}, {21, -30}}),
Orientation::counterclockwise));

EXPECT_TRUE(has_orientation(
make_span(VecReal2{{-2, -2}, {0, -2}, {0, 0}, {-2, 0}}),
Orientation::counterclockwise));
}

TEST(PolygonUtilsTest, convexity)
{
VecReal2 cw{{1, 1}, {1, 2}, {2, 2}, {2, 1}};
Expand Down

0 comments on commit d6720d6

Please sign in to comment.