From 746badbc4c4aa1d8939967ca3a2a1cfc15257c12 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 19 Oct 2021 16:02:49 -0400 Subject: [PATCH 01/22] Add types and quadratic solver --- src/CMakeLists.txt | 1 + src/geometry/GeoParams.hh | 1 + src/geometry/Types.hh | 3 +- src/orange/Types.cc | 15 ++ src/orange/Types.hh | 254 ++++++++++++++++++ src/orange/surfaces/detail/QuadraticSolver.hh | 196 ++++++++++++++ test/CMakeLists.txt | 7 + .../surfaces/detail/QuadraticSolver.test.cc | 163 +++++++++++ 8 files changed, 638 insertions(+), 2 deletions(-) create mode 100644 src/orange/Types.cc create mode 100644 src/orange/Types.hh create mode 100644 src/orange/surfaces/detail/QuadraticSolver.hh create mode 100644 test/orange/surfaces/detail/QuadraticSolver.test.cc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9330cd4a11..5d839902b1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -61,6 +61,7 @@ list(APPEND SOURCES comm/ScopedMpiInit.cc comm/detail/LoggerMessage.cc geometry/detail/ScopedTimeAndRedirect.cc + orange/Types.cc io/ImportProcess.cc io/ImportPhysicsTable.cc io/ImportPhysicsVector.cc diff --git a/src/geometry/GeoParams.hh b/src/geometry/GeoParams.hh index 5ede13bcbf..58518b98af 100644 --- a/src/geometry/GeoParams.hh +++ b/src/geometry/GeoParams.hh @@ -43,6 +43,7 @@ class GeoParams const std::string& id_to_label(VolumeId vol_id) const; // Get the ID corresponding to a label + // TODO: rename to find ?? see MaterialParams etc. VolumeId label_to_id(const std::string& label) const; //! Number of volumes diff --git a/src/geometry/Types.hh b/src/geometry/Types.hh index f17bd0f187..fc7705dd5d 100644 --- a/src/geometry/Types.hh +++ b/src/geometry/Types.hh @@ -13,9 +13,8 @@ namespace celeritas { -class Geometry; //---------------------------------------------------------------------------// -//! Opaque numeric identifier for a geometry cell +//! Identifier for a geometry volume using VolumeId = OpaqueId; //---------------------------------------------------------------------------// diff --git a/src/orange/Types.cc b/src/orange/Types.cc new file mode 100644 index 0000000000..ef172b9b3a --- /dev/null +++ b/src/orange/Types.cc @@ -0,0 +1,15 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file Types.cc +//---------------------------------------------------------------------------// +#include "Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/Types.hh b/src/orange/Types.hh new file mode 100644 index 0000000000..82e76c0299 --- /dev/null +++ b/src/orange/Types.hh @@ -0,0 +1,254 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file Types.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +#include "base/OpaqueId.hh" +#include "base/NumericLimits.hh" +#include "geometry/Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +// TYPE ALIASES +//---------------------------------------------------------------------------// + +//! Integer type for volume CSG tree representation +using logic_int = unsigned short int; + +//! Identifier for a surface in a universe +using SurfaceId = OpaqueId; + +//! Identifier for a relocatable set of volumes +using UniverseId = OpaqueId; + +//---------------------------------------------------------------------------// +// ENUMERATIONS +//---------------------------------------------------------------------------// +/*! + * Whether a position is logically "inside" or "outside" a surface. + * + * For a plane, "POS" (outside/true) is equivalent to + * \f[ + \vec x \cdot \vec n >= 0 + * \f] + * and "inside" (NEG) is to the left of the plane's normal. Likewise, for a + * sphere, "inside" is where the dot product of the position and outward normal + * is negative. These are *opposite* the signs from KENO, where `-` has the + * sense of negating a closed region of space (i.e. the interior of a sphere). + */ +enum class Sense : bool +{ + inside, //!< Quadric expression is less than zero + outside, //!< Expression is greater than zero +}; + +//---------------------------------------------------------------------------// +/*! + * Enumeration for cartesian axes. + */ +enum class Axis +{ + x, + y, + z, + size_ +}; + +//---------------------------------------------------------------------------// +/*! + * Enumeration for mapping surface classes to integers. + * + * These are ordered by number of coefficients needed for their representation: + * 1 for `[ps].|c.o`, 3 for `c.`, 4 for `[ps]|k.`, 7 for `sq`, and 10 for `gq`. + */ +enum class SurfaceType : unsigned char +{ + px, //!< Plane aligned with X axis + py, //!< Plane aligned with Y axis + pz, //!< Plane aligned with Z axis + cxc, //!< Cylinder centered on X axis + cyc, //!< Cylinder centered on Y axis + czc, //!< Cylinder centered on Z axis +#if 0 + sc, //!< Sphere centered at the origin + cx, //!< Cylinder parallel to X axis + cy, //!< Cylinder parallel to Y axis + cz, //!< Cylinder parallel to Z axis + p, //!< General plane + s, //!< Sphere + kx, //!< Cone parallel to X axis + ky, //!< Cone parallel to Y axis + kz, //!< Cone parallel to Z axis + sq, //!< Simple quadric +#endif + gq, //!< General quadric + size_ +}; + +//---------------------------------------------------------------------------// +/*! + * . + * + * For a plane, "outside" is equivalent to + * \f[ + \vec x \cdot \vec n > 0 + * \f] + * and "inside" is to the left of the plane's normal (a negative dot product). + * The exact equality to zero is literally an "edge case" but it can happen + * with inter-universe coincident surfaces as well as carefully placed + * particle sources and ray tracing. + * + * As an implementataion detail, the "on" case is currently *exact*, but future + * changes might increase the width of "on" to a finite but small range + * ("fuzziness"). + */ +enum class SignedSense +{ + inside = -1, + on = 0, + outside = 1 +}; + +//---------------------------------------------------------------------------// +// CLASSES +//---------------------------------------------------------------------------// +/*! + * \struct Initialization + * \brief Volume ID and surface ID after initialization + * + * Possible configurations for the initialization result ('X' means 'has + * a valid ID', i.e. evaluates to true): + * + * Vol | Surface | Description + * :----: | :-----: | :------------------------------- + * | | Failed to find new vol + * X | | Initialized + * X | X | Crossed surface into new vol + * | X | Initialized on a surface (reject) + */ +struct Initialization +{ + VolumeId volume; + SurfaceId surface; + Sense sense = Sense::inside; // Sense if on a surface + + //! Whether initialization succeeded + explicit CELER_FUNCTION operator bool() const + { + return static_cast(volume); + } +}; + +//---------------------------------------------------------------------------// +// HELPER FUNCTIONS (HOST/DEVICE) +//---------------------------------------------------------------------------// +/*! + * Convert a boolean value to a Sense enum. + */ +CELER_CONSTEXPR_FUNCTION Sense to_sense(bool s) +{ + return static_cast(s); +} + +//---------------------------------------------------------------------------// +/*! + * Change the sense across a surface. + */ +CELER_CONSTEXPR_FUNCTION Sense flip_sense(Sense orig) +{ + return static_cast(!static_cast(orig)); +} + +//---------------------------------------------------------------------------// +/*! + * Evaluate the sense based on the LHS expression of the quadric equation. + * + * This is an optimized jump-free version of: + * \code + return quadric == 0 ? SignedSense::on + : quadric < 0 ? SignedSense::inside + : SignedSense::outside; + * \endcode + * as + * \code + int gz = !(quadric <= 0) ? 1 : 0; + int lz = quadric < 0 ? 1 : 0; + return static_cast(gz - lz); + * \endcode + * and compressed into a single line. + * + * NaN values are treated as "outside". + */ +CELER_CONSTEXPR_FUNCTION SignedSense real_to_sense(real_type quadric) +{ + return static_cast(!(quadric <= 0) - (quadric < 0)); +} + +//---------------------------------------------------------------------------// +/*! + * Convert a signed sense to a Sense enum. + */ +CELER_CONSTEXPR_FUNCTION Sense to_sense(SignedSense s) +{ + return Sense(static_cast(s) >= 0); +} + +//---------------------------------------------------------------------------// +/*! + * Sentinel value indicating "no intersection". + * + * \todo There is probably a better place to put this since it's not a "type". + */ +CELER_CONSTEXPR_FUNCTION real_type no_intersection() +{ + return numeric_limits::infinity(); +} + +//---------------------------------------------------------------------------// +// HELPER FUNCTIONS (HOST) +//---------------------------------------------------------------------------// + +//! Get a printable character corresponding to a sense. +inline static constexpr char to_char(Sense s) +{ + return s == Sense::inside ? '-' : '+'; +} + +//! Get the lowercase name of the axis. +inline static constexpr char to_char(Axis ax) +{ + return "xyz\a"[static_cast(ax)]; +} + +// Get a string corresponding to a surface type +const char* to_cstring(SurfaceType); + +//---------------------------------------------------------------------------// +} // namespace celeritas + +//---------------------------------------------------------------------------// +// STD::HASH SPECIALIZATION FOR HOST CODE +//---------------------------------------------------------------------------// +//! \cond +namespace std +{ +//! Specialization for std::hash for unordered storage. +template<> +struct hash +{ + using argument_type = celeritas::Sense; + using result_type = std::size_t; + result_type operator()(const argument_type& sense) const noexcept + { + return std::hash()(static_cast(sense)); + } +}; +} // namespace std +//! \endcond diff --git a/src/orange/surfaces/detail/QuadraticSolver.hh b/src/orange/surfaces/detail/QuadraticSolver.hh new file mode 100644 index 0000000000..649e16c081 --- /dev/null +++ b/src/orange/surfaces/detail/QuadraticSolver.hh @@ -0,0 +1,196 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file QuadraticSolver.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +#include "base/Algorithms.hh" +#include "base/Array.hh" +#include "base/Types.hh" +#include "orange/Types.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Find positive, real, nonzero roots for quadratic functions. + * + * These are for quadratic functions \f[ + a x^2 + b^2 + c = 0 + * \f] + * where a is nonzero (and not close to zero). + * + * This is used for all quadrics with potentially two roots (anything but + * planes). + * + * Each item in the Intersections result will be a positive valid intersection + * or the sentinel result \c no_intersection() . + */ +class QuadraticSolver +{ + public: + //!@{ + //! Type aliases + using Intersections = Array; + //!@} + + //! Fuzziness for "along surface" + static CELER_CONSTEXPR_FUNCTION real_type min_a() { return 1e-10; } + + // Solve when possibly along a surface (zeroish a) + static inline CELER_FUNCTION Intersections solve_general(real_type a, + real_type half_b, + real_type c, + bool on_surface); + + public: + // Construct with nonzero a, and b/2 + inline CELER_FUNCTION QuadraticSolver(real_type a, real_type half_b); + + // Solve fully general case + inline CELER_FUNCTION Intersections operator()(real_type c) const; + + // Solve degenerate case (known to be on surface) + inline CELER_FUNCTION Intersections operator()() const; + + private: + //// DATA //// + real_type a_inv_; + real_type hba_; // (b/2)/a +}; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Find all nonnegative roots for general quadric surfaces. + * + * This is used for cones, simple quadrics, and general quadrics. + */ +CELER_FUNCTION auto QuadraticSolver::solve_general(real_type a, + real_type half_b, + real_type c, + bool on_surface) + -> Intersections +{ + if (std::fabs(a) >= min_a()) + { + // Not along the surface + QuadraticSolver solve(a, half_b); + return on_surface ? solve() : solve(c); + } + else + { + // Travelling parallel to the quadric's surface + if (!on_surface) + { + QuadraticSolver solve(min_a(), half_b); + return solve(c); + } + else + { + // On and along surface: no intersections + return {no_intersection(), no_intersection()}; + } + } +} + +//---------------------------------------------------------------------------// +/*! + * Construct with b/2. + */ +CELER_FUNCTION QuadraticSolver::QuadraticSolver(real_type a, real_type half_b) + : a_inv_(1 / a), hba_(half_b * a_inv_) +{ + CELER_EXPECT(std::fabs(a) >= min_a()); +} + +//---------------------------------------------------------------------------// +/*! + * Find all nonnegative roots of x^2 + (b/a)*x + (c/a) = 0. + * + * In this case, the quadratic formula can be written as: \f[ + x = -b/2 \pm \sqrt{(b/2)^2 - c}. + * \f] + * + * Callers: + * - General quadratic solve: not on nor along surface + * - Sphere when not on surface + * - Cylinder when not on surface + */ +CELER_FUNCTION auto QuadraticSolver::operator()(real_type c) const + -> Intersections +{ + // Scale c by 1/a in accordance with scaling of b + c *= a_inv_; + real_type b2_4 = ipow<2>(hba_); // (b/2)^2 + + Intersections result; + + if (b2_4 > c) + { + // Two real roots, r1 and r2 + real_type t2 = std::sqrt(b2_4 - c); // (b^2 - 4ac) / 4 + result[0] = -hba_ - t2; + result[1] = -hba_ + t2; + + if (result[1] < 0) + { + // Both are negative + result[0] = no_intersection(); + result[1] = no_intersection(); + } + else if (result[0] < 0) + { + // Only first is negative + result[0] = no_intersection(); + } + } + else if (b2_4 == c) + { + // One real root, r1 + result[0] = -hba_; + result[1] = no_intersection(); + + if (result[0] < 0) + { + result[0] = no_intersection(); + } + } + else + { + // No real roots + result = {no_intersection(), no_intersection()}; + } + + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Solve degenerate case (known to be on surface). + * + * Since x = 0 is a root, then c = 0 and x = -b is the other root. + */ +CELER_FUNCTION auto QuadraticSolver::operator()() const -> Intersections +{ + Intersections result{-2 * hba_, no_intersection()}; + + if (result[0] <= 0) + { + result[0] = no_intersection(); + } + + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ecabd09e8a..4fe4a4721d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -178,6 +178,13 @@ if(CELERITAS_USE_VecGeom) celeritas_add_test(geometry/LinearPropagator.test.cc GPU) endif() +#-----------------------------------------------------------------------------# +# ORANGE geometry + +celeritas_setup_tests(SERIAL PREFIX orange) + +celeritas_add_test(orange/surfaces/detail/QuadraticSolver.test.cc) + #-----------------------------------------------------------------------------# # I/O (ROOT) diff --git a/test/orange/surfaces/detail/QuadraticSolver.test.cc b/test/orange/surfaces/detail/QuadraticSolver.test.cc new file mode 100644 index 0000000000..14ff7863f0 --- /dev/null +++ b/test/orange/surfaces/detail/QuadraticSolver.test.cc @@ -0,0 +1,163 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file QuadraticSolver.test.cc +//---------------------------------------------------------------------------// +#include "orange/surfaces/detail/QuadraticSolver.hh" + +#include "celeritas_test.hh" + +using celeritas::no_intersection; +using celeritas::detail::QuadraticSolver; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +TEST(SolveNonsurface, no_roots) +{ + // x^2 + 2*x + 1 = 0 -> one real root (x = -1) + // x^2 + 2*x + 1.001 = 0 -> two complex roots + { + double b_2 = 1.0; // b/2 + double c = 1.001; + + QuadraticSolver solve_quadratic(1, b_2); + auto x = solve_quadratic(c); + + // Verify that x was not changed + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_EQ(no_intersection(), x[1]); + } + + // x^2 - 4*x + 4 = 0 -> one real root (x = 2) + // x^2 - 4*x + 4.001 = 0 -> two complex roots + { + double b_2 = -2.0; + double c = 4.001; + + QuadraticSolver solve_quadratic(1, b_2); + auto x = solve_quadratic(c); + + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_EQ(no_intersection(), x[1]); + } +} + +TEST(SolveNonsurface, one_root) +{ + // x^2 - 2*x + 1 = 0 -> x = 1 + double b_2 = -1.0; + double c = 1.0; + + // For solve_quadratic() to detect the single root case, it must calculate + // that b/2 * b/2 - c == 0. It is assumed here that the floating point + // operations -1.*-1. - 1. will reliably yield 0. + + QuadraticSolver solve_quadratic(1, b_2); + auto x = solve_quadratic(c); + + EXPECT_SOFT_EQ(1.0, x[0]); + EXPECT_SOFT_EQ(no_intersection(), x[1]); +} + +TEST(SolveNonsurface, two_roots) +{ + // x^2 - 3*x + 2 = 0 -> x = 1, 2 + { + double b_2 = -1.5; + double c = 2.0; + + auto x = QuadraticSolver(1, b_2)(c); + + EXPECT_SOFT_EQ(1.0, x[0]); + EXPECT_SOFT_EQ(2.0, x[1]); + } + + // x^2 + x - 20 = 0 -> x = -5, 4 + { + double b_2 = 0.5; + double c = -20.0; + + auto x = QuadraticSolver(1, b_2)(c); + + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_EQ(4.0, x[1]); + } + + // x^2 + 3*x + 2 = 0 -> x = -1, -2 + { + double b_2 = 3.0 / 2; + double c = 2.0; + + auto x = QuadraticSolver(1, b_2)(c); + + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_EQ(no_intersection(), x[1]); + } + + // x^2 - 99999.999*x - 100 = 0 -> x = -0.001, 100000 + { + double b_2 = -99999.999 / 2.0; + double c = -100.0; + + auto x = QuadraticSolver(1, b_2)(c); + + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_EQ(100000.0, x[1]); + } +} + +TEST(SolveSurface, one_root) +{ + // x^2 + 2*x + 0 = 0, x -> -2 + { + double b_2 = 2.0 / 2; + + auto x = QuadraticSolver(1, b_2)(); + + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_EQ(no_intersection(), x[1]); + } + + // x^2 - 3.45*x + 0 = 0, x -> 3.45 + { + double b_2 = -3.45 / 2; + + auto x = QuadraticSolver(1, b_2)(); + + EXPECT_SOFT_EQ(3.45, x[0]); + EXPECT_SOFT_EQ(no_intersection(), x[1]); + } +} + +TEST(SolveGeneral, no_roots) +{ + // For a*x^2 + b*x + c = 0, as both a and b -> 0, |x| -> infinity. In this + // case, solve_degenerate_quadratic will return no positive roots. + // + // -1.0e-15*x^2 + 10000 = 0 -> x ~= +/- 3e9 + double a = -1e-15; + double b_2 = 0.; // (b/a)/2 + double c = 10000; // c/a + + auto x = QuadraticSolver::solve_general(a, b_2, c, false); + + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_EQ(no_intersection(), x[1]); +} + +TEST(SolveGeneral, one_root) +{ + // 1.0e-15*x^2 + 2*x - 2000 = 0 -> x = 1e3 + double a = 1e-15; + double b_2 = 2.0 / 2; + double c = -2000; + + auto x = QuadraticSolver::solve_general(a, b_2, c, false); + + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_NEAR(1.0e3, x[1], 1e-7); +} From fadcf73c693a6dbc639ab3c073048804283c3e93 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 19 Oct 2021 17:04:40 -0400 Subject: [PATCH 02/22] Add CylCentered --- src/orange/Types.hh | 18 +- src/orange/surfaces/CylCentered.hh | 230 +++++++++ src/orange/surfaces/detail/QuadraticSolver.hh | 18 +- test/CMakeLists.txt | 1 + test/orange/surfaces/CylCentered.test.cc | 455 ++++++++++++++++++ .../surfaces/detail/QuadraticSolver.test.cc | 5 +- 6 files changed, 712 insertions(+), 15 deletions(-) create mode 100644 src/orange/surfaces/CylCentered.hh create mode 100644 test/orange/surfaces/CylCentered.test.cc diff --git a/src/orange/Types.hh b/src/orange/Types.hh index 82e76c0299..95c2fa69ec 100644 --- a/src/orange/Types.hh +++ b/src/orange/Types.hh @@ -94,7 +94,7 @@ enum class SurfaceType : unsigned char //---------------------------------------------------------------------------// /*! - * . + * Evaluated quadric expression allowing for distinct 'on surface' state. * * For a plane, "outside" is equivalent to * \f[ @@ -116,6 +116,18 @@ enum class SignedSense outside = 1 }; +//---------------------------------------------------------------------------// +/*! + * When evaluating an intersection, whether the point is on the surface. + * + * This helps eliminate roundoff errors and other arithmetic issues. + */ +enum class SurfaceState : bool +{ + off = false, + on = true +}; + //---------------------------------------------------------------------------// // CLASSES //---------------------------------------------------------------------------// @@ -128,9 +140,9 @@ enum class SignedSense * * Vol | Surface | Description * :----: | :-----: | :------------------------------- - * | | Failed to find new vol + * | | Failed to find new volume * X | | Initialized - * X | X | Crossed surface into new vol + * X | X | Crossed surface into new volume * | X | Initialized on a surface (reject) */ struct Initialization diff --git a/src/orange/surfaces/CylCentered.hh b/src/orange/surfaces/CylCentered.hh new file mode 100644 index 0000000000..69ba265f9c --- /dev/null +++ b/src/orange/surfaces/CylCentered.hh @@ -0,0 +1,230 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file CylCentered.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Array.hh" +#include "base/ArrayUtils.hh" +#include "base/Assert.hh" +#include "base/Macros.hh" +#include "base/Span.hh" +#include "../Types.hh" +#include "detail/QuadraticSolver.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Axis-aligned cylinder centered about the origin. + * + * The cylinder is centered along an Axis template parameter. + * + * For a cylinder along the x axis: + * \f[ + y^2 + z^2 - R^2 = 0 + \f] + * + * This is an optimization of the Cyl. The motivations are: + * - Many geometries have units with concentric cylinders centered about the + * origin, so having this as a special case reduces the memory usage of those + * units (improved memory localization). + * - Cylindrical mesh geometries have lots of these cylinders, so efficient + * tracking through its cells should make this optimization worthwhile. + */ +template +class CylCentered +{ + public: + //@{ + //! Type aliases + using Intersections = Array; + using SpanConstRealN = Span; + //@} + + //// CLASS ATTRIBUTES //// + + // Surface type identifier + static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type(); + + //! Storage requirements + static CELER_CONSTEXPR_FUNCTION size_type size() + { + return SpanConstRealN::extent; + } + + public: + //// CONSTRUCTORS //// + + // Construct with radius + explicit inline CELER_FUNCTION CylCentered(real_type radius); + + // Construct from raw data + explicit inline CELER_FUNCTION CylCentered(SpanConstRealN); + + //// ACCESSORS //// + + //! Get the square of the radius + CELER_FUNCTION real_type radius_sq() const { return radius_sq_; } + + //! Get a view to the data for type-deleted storage + CELER_FUNCTION SpanConstRealN data() const { return {&radius_sq_, 1}; } + + //// CALCULATION //// + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION SignedSense calc_sense(const Real3& pos) const; + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION Intersections calc_intersections( + const Real3& pos, const Real3& dir, SurfaceState on_surface) const; + + // Calculate outward normal at a position + inline CELER_FUNCTION Real3 calc_normal(const Real3& pos) const; + + private: + //! Square of cylinder radius + real_type radius_sq_; + + static CELER_CONSTEXPR_FUNCTION int t_index(); + static CELER_CONSTEXPR_FUNCTION int u_index(); + static CELER_CONSTEXPR_FUNCTION int v_index(); +}; + +//---------------------------------------------------------------------------// +// TYPE ALIASES +//---------------------------------------------------------------------------// + +using CCylX = CylCentered; +using CCylY = CylCentered; +using CCylZ = CylCentered; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Surface type identifier. + */ +template +CELER_CONSTEXPR_FUNCTION SurfaceType CylCentered::surface_type() +{ + return (T == Axis::x + ? SurfaceType::cxc + : (T == Axis::y ? SurfaceType::cyc : SurfaceType::czc)); +} + +//---------------------------------------------------------------------------// +/*! + * Construct with radius. + */ +template +CELER_FUNCTION CylCentered::CylCentered(real_type radius) + : radius_sq_(ipow<2>(radius)) +{ + CELER_EXPECT(radius > 0); +} + +//---------------------------------------------------------------------------// +/*! + * Construct from raw data. + */ +template +CELER_FUNCTION CylCentered::CylCentered(SpanConstRealN data) + : radius_sq_(data[0]) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Determine the sense of the position relative to this surface. + */ +template +CELER_FUNCTION SignedSense CylCentered::calc_sense(const Real3& pos) const +{ + const real_type u = pos[u_index()]; + const real_type v = pos[v_index()]; + + return real_to_sense(ipow<2>(u) + ipow<2>(v) - radius_sq_); +} + +//---------------------------------------------------------------------------// +/*! + * Determine the sense of the position relative to this surface. + */ +template +CELER_FUNCTION auto +CylCentered::calc_intersections(const Real3& pos, + const Real3& dir, + SurfaceState on_surface) const + -> Intersections +{ + // 1 - \omega \dot e + const real_type a = 1 - ipow<2>(dir[t_index()]); + + if (a != 0) + { + const real_type u = pos[u_index()]; + const real_type v = pos[v_index()]; + + // b/2 = \omega \dot (x - x_0) + detail::QuadraticSolver solve_quadric( + a, dir[u_index()] * u + dir[v_index()] * v); + if (on_surface == SurfaceState::off) + { + // c = (x - x_0) \dot (x - x_0) - R * R + return solve_quadric(ipow<2>(u) + ipow<2>(v) - radius_sq_); + } + else + { + // Solve degenerate case (c=0) + return solve_quadric(); + } + } + else + { + // No intersection if we're traveling along the cylinder axis + return {no_intersection(), no_intersection()}; + } +} + +//---------------------------------------------------------------------------// +/*! + * Calculate outward normal at a position. + */ +template +CELER_FUNCTION Real3 CylCentered::calc_normal(const Real3& pos) const +{ + Real3 norm{0, 0, 0}; + + norm[u_index()] = pos[u_index()]; + norm[v_index()] = pos[v_index()]; + + normalize_direction(&norm); + return norm; +} + +//---------------------------------------------------------------------------// +//!@{ +//! Integer index values for primary and orthogonal axes. +template +CELER_CONSTEXPR_FUNCTION int CylCentered::t_index() +{ + return static_cast(T); +} +template +CELER_CONSTEXPR_FUNCTION int CylCentered::u_index() +{ + return static_cast(T == Axis::x ? Axis::y : Axis::x); +} +template +CELER_CONSTEXPR_FUNCTION int CylCentered::v_index() +{ + return static_cast(T == Axis::z ? Axis::y : Axis::z); +} +//!@} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/surfaces/detail/QuadraticSolver.hh b/src/orange/surfaces/detail/QuadraticSolver.hh index 649e16c081..fdbd9086c5 100644 --- a/src/orange/surfaces/detail/QuadraticSolver.hh +++ b/src/orange/surfaces/detail/QuadraticSolver.hh @@ -45,10 +45,8 @@ class QuadraticSolver static CELER_CONSTEXPR_FUNCTION real_type min_a() { return 1e-10; } // Solve when possibly along a surface (zeroish a) - static inline CELER_FUNCTION Intersections solve_general(real_type a, - real_type half_b, - real_type c, - bool on_surface); + static inline CELER_FUNCTION Intersections solve_general( + real_type a, real_type half_b, real_type c, SurfaceState on_surface); public: // Construct with nonzero a, and b/2 @@ -74,22 +72,22 @@ class QuadraticSolver * * This is used for cones, simple quadrics, and general quadrics. */ -CELER_FUNCTION auto QuadraticSolver::solve_general(real_type a, - real_type half_b, - real_type c, - bool on_surface) +CELER_FUNCTION auto QuadraticSolver::solve_general(real_type a, + real_type half_b, + real_type c, + SurfaceState on_surface) -> Intersections { if (std::fabs(a) >= min_a()) { // Not along the surface QuadraticSolver solve(a, half_b); - return on_surface ? solve() : solve(c); + return on_surface == SurfaceState::on ? solve() : solve(c); } else { // Travelling parallel to the quadric's surface - if (!on_surface) + if (on_surface == SurfaceState::off) { QuadraticSolver solve(min_a(), half_b); return solve(c); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4fe4a4721d..d837c43dec 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -184,6 +184,7 @@ endif() celeritas_setup_tests(SERIAL PREFIX orange) celeritas_add_test(orange/surfaces/detail/QuadraticSolver.test.cc) +celeritas_add_test(orange/surfaces/CylCentered.test.cc) #-----------------------------------------------------------------------------# # I/O (ROOT) diff --git a/test/orange/surfaces/CylCentered.test.cc b/test/orange/surfaces/CylCentered.test.cc new file mode 100644 index 0000000000..542d40cda8 --- /dev/null +++ b/test/orange/surfaces/CylCentered.test.cc @@ -0,0 +1,455 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file CylCentered.test.cc +//---------------------------------------------------------------------------// +#include "orange/surfaces/CylCentered.hh" + +#include +#include +#include +#include "base/Algorithms.hh" +#include "celeritas_test.hh" + +using celeritas::CCylX; +using celeritas::CCylY; +using celeritas::CCylZ; +using celeritas::no_intersection; +using celeritas::SignedSense; +using celeritas::SurfaceState; + +using celeritas::ipow; +using celeritas::Real3; +using celeritas::real_type; + +using Intersections = CCylX::Intersections; +using VecReal = std::vector; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// +TEST(TestCCylX, construction) +{ + EXPECT_EQ(1, CCylX::size()); + EXPECT_EQ(2, CCylX::Intersections{}.size()); + + CCylX c(4.0); + + const real_type expected_data[] = {ipow<2>(4)}; + + EXPECT_VEC_SOFT_EQ(expected_data, c.data()); +} + +TEST(TestCCylX, sense) +{ + CCylX cyl(4.0); + + EXPECT_EQ(SignedSense::inside, cyl.calc_sense(Real3{0, 3, 0})); + EXPECT_EQ(SignedSense::outside, cyl.calc_sense(Real3{0, 0, 9})); +} + +TEST(TestCCylX, normal) +{ + CCylX cyl(3.45); + + EXPECT_VEC_SOFT_EQ((Real3{0, 1, 0}), cyl.calc_normal(Real3{1.23, 3.45, 0})); + EXPECT_VEC_SOFT_EQ((Real3{0, 0, 1}), + cyl.calc_normal(Real3{0.0, 0.0, 4.45})); + + // Normal from center of cylinder is ill-defined but shouldn't raise an + // error + auto norm = cyl.calc_normal(Real3{10.0, 0.0, 0.0}); + EXPECT_TRUE(std::isnan(norm[0])); + EXPECT_TRUE(std::isnan(norm[1])); + EXPECT_TRUE(std::isnan(norm[2])); +} + +TEST(TestCCylX, intersect) +{ + Intersections distances{-1, -1}; + + // From inside + CCylX cyl(3.0); + distances = cyl.calc_intersections( + Real3{0, 0, 1.5}, Real3{0, 1, 0}, SurfaceState::off); + + EXPECT_SOFT_EQ(2.598076211353316, distances[1]); + EXPECT_EQ(no_intersection(), distances[0]); + + // From outside, hitting both + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{0, 0, -5.0}, Real3{0, 0, 1}, SurfaceState::off); + + EXPECT_SOFT_EQ(2.0, distances[0]); + EXPECT_SOFT_EQ(8.0, distances[1]); + + // From outside, tangent + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{0, -3, -5.0}, Real3{0, 0, 1}, SurfaceState::off); + + EXPECT_SOFT_EQ(5.0, distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + // From outside, hitting neither + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{0, 0, -5.0}, Real3{0, 1, 0}, SurfaceState::off); + + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); +} + +TEST(TestCCylX, intersect_from_surface) +{ + Intersections distances; + + CCylX cyl(3.45); + + // One intercept + + distances = cyl.calc_intersections( + Real3{1.23, 3.45, 0.0}, Real3{0, -1, 0}, SurfaceState::on); + EXPECT_SOFT_EQ(6.9, distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + // No intercepts + distances = cyl.calc_intersections( + Real3{1.23, 4.68, 2.34}, Real3{0, 1, 0}, SurfaceState::on); + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); +} + +//---------------------------------------------------------------------------// + +TEST(TestCCylY, sense) +{ + CCylY cyl(3.0); + + EXPECT_EQ(SignedSense::inside, cyl.calc_sense(Real3{1.5, 0, 0})); + EXPECT_EQ(SignedSense::outside, cyl.calc_sense(Real3{3.01, 0, 0})); +} + +TEST(TestCCylY, intersect) +{ + CCylY::Intersections distances{-1, -1}; + + // From inside + CCylY cyl(3.0); + distances = cyl.calc_intersections( + Real3{1.5, 0, 0}, Real3{0, 0, 1}, SurfaceState::off); + + EXPECT_SOFT_EQ(2.598076211353316, distances[1]); + EXPECT_EQ(no_intersection(), distances[0]); + + // From outside, hitting both + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{-5.0, 0, 0}, Real3{1, 0, 0}, SurfaceState::off); + + EXPECT_SOFT_EQ(2.0, distances[0]); + EXPECT_SOFT_EQ(8.0, distances[1]); + + // From outside, tangent + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{-5.0, 0, -3}, Real3{1, 0, 0}, SurfaceState::off); + + EXPECT_SOFT_EQ(5.0, distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + // From outside, hitting neither + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{-5.0, 0, 0}, Real3{0, 0, 1}, SurfaceState::off); + + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); +} + +TEST(TestCCylY, intersect_from_surface) +{ + CCylY::Intersections distances; + + CCylY cyl(3.45); + + // One intercept + + distances = cyl.calc_intersections( + Real3{3.45, 1.23, 0}, Real3{-1, 0, 0}, SurfaceState::on); + EXPECT_SOFT_EQ(6.9, distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + distances = cyl.calc_intersections( + Real3{1.5528869044748548, 12345., 3.08577380894971}, + Real3{0, 0.6, -0.8}, + SurfaceState::on); + EXPECT_SOFT_EQ(7.714434522374273, distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + // No intercepts + + distances = cyl.calc_intersections( + Real3{3.45, 1.23, 0.0}, Real3{1, 0, 0}, SurfaceState::on); + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + distances = cyl.calc_intersections( + Real3{1.5528869044748548, 12345., 3.08577380894971}, + Real3{0, 0.6, 0.8}, + SurfaceState::on); + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); +} + +//---------------------------------------------------------------------------// + +TEST(TestCCylZ, sense) +{ + CCylZ cyl(3.0); + + EXPECT_EQ(SignedSense::inside, cyl.calc_sense(Real3{1.5, 0, 0})); + EXPECT_EQ(SignedSense::outside, cyl.calc_sense(Real3{3.01, 0, 0})); +} + +TEST(TestCCylZ, calc_intersections) +{ + Intersections distances{-1, -1}; + + // From inside + CCylZ cyl(3.0); + distances = cyl.calc_intersections( + Real3{1.5, 0, 0}, Real3{0, 1, 0}, SurfaceState::off); + + EXPECT_SOFT_EQ(2.598076211353316, distances[1]); + EXPECT_EQ(no_intersection(), distances[0]); + + // From outside, hitting both + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{-5.0, 0, 0}, Real3{1, 0, 0}, SurfaceState::off); + + EXPECT_SOFT_EQ(2.0, distances[0]); + EXPECT_EQ(8.0, distances[1]); + + // From outside, tangent + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{-5.0, -3, 0}, Real3{1, 0, 0}, SurfaceState::off); + + EXPECT_SOFT_EQ(5.0, distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + // From outside, hitting neither + distances[0] = distances[1] = -1; + distances = cyl.calc_intersections( + Real3{-5.0, 0, 0}, Real3{0, 1, 0}, SurfaceState::off); + + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); +} + +TEST(TestCCylZ, calc_intersections_on_surface) +{ + CCylZ::Intersections distances; + const real_type eps = 1.e-4; + + { + CCylZ cyl(1.0); + + // Heading toward, slightly inside + distances = cyl.calc_intersections( + Real3{-1 + eps, 0, 0}, Real3{1, 0, 0}, SurfaceState::on); + EXPECT_SOFT_NEAR(2.0 - eps, distances[0], eps); + EXPECT_EQ(no_intersection(), distances[1]); + + // Heading away, slightly inside + distances = cyl.calc_intersections( + Real3{-1 + eps, 0, 0}, Real3{-1, 0, 0}, SurfaceState::on); + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + // Tangent + distances = cyl.calc_intersections( + Real3{-1 + eps, 0, 0}, Real3{0, 1, 0}, SurfaceState::on); + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + // Heading in, slightly inside + distances = cyl.calc_intersections( + Real3{-1 + eps, 0, 0}, Real3{1, 0, 0}, SurfaceState::on); + EXPECT_SOFT_NEAR(2.0 - eps, distances[0], eps); + EXPECT_EQ(no_intersection(), distances[1]); + + // Heading away, slightly outside + distances = cyl.calc_intersections( + Real3{1.0 - eps, 0, 0}, Real3{1, 0, 0}, SurfaceState::on); + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + + // Tangent, slightly inside + distances = cyl.calc_intersections( + Real3{-1 + eps, 0, 0}, Real3{0, 1, 0}, SurfaceState::on); + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); + } +} + +// Fused multiply-add on some CPUs in opt mode can cause the results of +// nearly-tangent cylinder checking to change. +TEST(TestCCylZ, multi_intersect) +{ + constexpr int Y = static_cast(celeritas::Axis::y); + + CCylZ cyl(10.0); + + VecReal all_first_distances; + VecReal all_distances; + VecReal all_y; + for (real_type x : {-10 + 1e-7, -1e-7, -0.0, 0.0, 1e-7, 10 - 1e-7}) + { + for (real_type v : {-1.0, 1.0}) + { + Real3 pos{x, -10.0001 * v, 0}; + Real3 dir{0, v, 0}; + + Intersections d; + + // Transport to inside of cylinder + d = cyl.calc_intersections(pos, dir, SurfaceState::off); + ASSERT_NE(no_intersection(), d[0]); + all_first_distances.push_back(d[0]); + pos[Y] += d[0] * dir[Y]; + all_y.push_back(pos[Y]); + + // Transport to other side of cylinder + d = cyl.calc_intersections(pos, dir, SurfaceState::on); + all_distances.push_back(d[0]); + ASSERT_NE(no_intersection(), d[0]); + pos[Y] += d[0] * dir[Y]; + + // We're done + d = cyl.calc_intersections(pos, dir, SurfaceState::on); + EXPECT_EQ(no_intersection(), d[0]); + } + } + + const real_type expected_all_first_distances[] = {9.99869, + 9.99869, + 0.0001, + 0.0001, + 0.0001, + 0.0001, + 0.0001, + 0.0001, + 0.0001, + 0.0001, + 9.99869, + 9.99869}; + EXPECT_VEC_NEAR(expected_all_first_distances, all_first_distances, 1e-5); + + const real_type expected_all_y[] = {0.00141421, + -0.00141421, + 10, + -10, + 10, + -10, + 10, + -10, + 10, + -10, + 0.00141421, + -0.00141421}; + EXPECT_VEC_NEAR(expected_all_y, all_y, 1e-5); + + const real_type expected_all_distances[] = {0.00282843, + 0.00282843, + 20, + 20, + 20, + 20, + 20, + 20, + 20, + 20, + 0.00282843, + 0.00282843}; + EXPECT_VEC_NEAR(expected_all_distances, all_distances, 1e-5); +} + +//---------------------------------------------------------------------------// +/*! + * Test initialization on or near boundary + */ +class DegenerateBoundaryTest : public celeritas::Test +{ + protected: + void run(real_type xdir) const; + + void run_all() + { + for (real_type r : {0.93, 1.0}) + { + radius = r; + for (real_type dir : {1, -1}) + { + std::ostringstream msg; + msg << "r=" << radius << ", dir_x=" << dir; + SCOPED_TRACE(msg.str()); + + run(1); + run(-1); + } + } + } + + protected: + real_type radius = -1; + real_type eps = -1; +}; + +void DegenerateBoundaryTest::run(real_type xdir) const +{ + CCylZ cyl(radius); + CCylZ::Intersections distances = {-1, -1}; + const real_type tol = std::max(1.e-14, 2 * std::fabs(eps)); + + // Distance across the cylinder + const real_type diameter = 2 * radius; + + Real3 pos = {0, 0, 0}; + Real3 dir = {xdir, 0, 0}; + + //// Inward boundary //// + pos[0] = -xdir * (diameter / 2 + eps); + distances = cyl.calc_intersections(pos, dir, SurfaceState::on); + EXPECT_SOFT_NEAR(diameter + eps, distances[0], tol); + EXPECT_EQ(no_intersection(), distances[1]); + + //// Outward boundary //// + pos[0] = xdir * (diameter / 2 + eps); + distances = cyl.calc_intersections(pos, dir, SurfaceState::on); + EXPECT_EQ(no_intersection(), distances[0]); + EXPECT_EQ(no_intersection(), distances[1]); +} + +TEST_F(DegenerateBoundaryTest, neg) +{ + eps = -1.0e-8; + run_all(); +} + +TEST_F(DegenerateBoundaryTest, DISABLED_zero) +{ + eps = 0.0; + run_all(); +} + +TEST_F(DegenerateBoundaryTest, pos) +{ + eps = 1.e-8; + run_all(); +} diff --git a/test/orange/surfaces/detail/QuadraticSolver.test.cc b/test/orange/surfaces/detail/QuadraticSolver.test.cc index 14ff7863f0..e56c69fdf6 100644 --- a/test/orange/surfaces/detail/QuadraticSolver.test.cc +++ b/test/orange/surfaces/detail/QuadraticSolver.test.cc @@ -10,6 +10,7 @@ #include "celeritas_test.hh" using celeritas::no_intersection; +using celeritas::SurfaceState; using celeritas::detail::QuadraticSolver; //---------------------------------------------------------------------------// @@ -143,7 +144,7 @@ TEST(SolveGeneral, no_roots) double b_2 = 0.; // (b/a)/2 double c = 10000; // c/a - auto x = QuadraticSolver::solve_general(a, b_2, c, false); + auto x = QuadraticSolver::solve_general(a, b_2, c, SurfaceState::off); EXPECT_SOFT_EQ(no_intersection(), x[0]); EXPECT_SOFT_EQ(no_intersection(), x[1]); @@ -156,7 +157,7 @@ TEST(SolveGeneral, one_root) double b_2 = 2.0 / 2; double c = -2000; - auto x = QuadraticSolver::solve_general(a, b_2, c, false); + auto x = QuadraticSolver::solve_general(a, b_2, c, SurfaceState::off); EXPECT_SOFT_EQ(no_intersection(), x[0]); EXPECT_SOFT_NEAR(1.0e3, x[1], 1e-7); From fda713e1169ee366c95aea3d0563aede92b2cf75 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 19 Oct 2021 17:12:45 -0400 Subject: [PATCH 03/22] Add aligned plane --- src/orange/surfaces/PlaneAligned.hh | 169 ++++++++++++++++++++++ test/CMakeLists.txt | 1 + test/orange/surfaces/PlaneAligned.test.cc | 133 +++++++++++++++++ 3 files changed, 303 insertions(+) create mode 100644 src/orange/surfaces/PlaneAligned.hh create mode 100644 test/orange/surfaces/PlaneAligned.test.cc diff --git a/src/orange/surfaces/PlaneAligned.hh b/src/orange/surfaces/PlaneAligned.hh new file mode 100644 index 0000000000..64cdb957f6 --- /dev/null +++ b/src/orange/surfaces/PlaneAligned.hh @@ -0,0 +1,169 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file PlaneAligned.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Array.hh" +#include "base/Assert.hh" +#include "base/Macros.hh" +#include "base/Span.hh" +#include "../Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Axis-aligned plane with positive-facing normal. + */ +template +class PlaneAligned +{ + public: + //@{ + //! Type aliases + using Intersections = Array; + using SpanConstRealN = Span; + //@} + + //// CLASS ATTRIBUTES //// + + // Surface type identifier + static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type(); + + //! Storage requirements + static CELER_CONSTEXPR_FUNCTION size_type size() + { + return SpanConstRealN::extent; + } + + public: + //// CONSTRUCTORS //// + + // Construct with radius + explicit inline CELER_FUNCTION PlaneAligned(real_type position); + + // Construct from raw data + explicit inline CELER_FUNCTION PlaneAligned(SpanConstRealN); + + //// ACCESSORS //// + + //! Get the square of the radius + CELER_FUNCTION real_type position() const { return position_; } + + //! Get a view to the data for type-deleted storage + CELER_FUNCTION SpanConstRealN data() const { return {&position_, 1}; } + + //// CALCULATION //// + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION SignedSense calc_sense(const Real3& pos) const; + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION Intersections calc_intersections( + const Real3& pos, const Real3& dir, SurfaceState on_surface) const; + + // Calculate outward normal at a position + inline CELER_FUNCTION Real3 calc_normal(const Real3& pos) const; + + private: + //! Intersection with the axis + real_type position_; + + static CELER_CONSTEXPR_FUNCTION int t_index() + { + return static_cast(T); + } +}; + +//---------------------------------------------------------------------------// +// TYPE ALIASES +//---------------------------------------------------------------------------// + +using PlaneX = PlaneAligned; +using PlaneY = PlaneAligned; +using PlaneZ = PlaneAligned; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Surface type identifier. + */ +template +CELER_CONSTEXPR_FUNCTION SurfaceType PlaneAligned::surface_type() +{ + return (T == Axis::x ? SurfaceType::px + : (T == Axis::y ? SurfaceType::py : SurfaceType::pz)); +} + +//---------------------------------------------------------------------------// +/*! + * Construct with radius. + */ +template +CELER_FUNCTION PlaneAligned::PlaneAligned(real_type position) + : position_(position) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Construct from raw data. + */ +template +CELER_FUNCTION PlaneAligned::PlaneAligned(SpanConstRealN data) + : position_(data[0]) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Determine the sense of the position relative to this surface. + */ +template +CELER_FUNCTION SignedSense PlaneAligned::calc_sense(const Real3& pos) const +{ + return real_to_sense(pos[t_index()] - position_); +} + +//---------------------------------------------------------------------------// +/*! + * Determine the sense of the position relative to this surface. + */ +template +CELER_FUNCTION auto +PlaneAligned::calc_intersections(const Real3& pos, + const Real3& dir, + SurfaceState on_surface) const + -> Intersections +{ + if (on_surface == SurfaceState::off && dir[t_index()] != 0) + { + real_type dist = (position_ - pos[t_index()]) / dir[t_index()]; + if (dist > 0) + { + return {dist}; + } + } + return {no_intersection()}; +} + +//---------------------------------------------------------------------------// +/*! + * Calculate outward normal at a position. + */ +template +CELER_FUNCTION Real3 PlaneAligned::calc_normal(const Real3& pos) const +{ + Real3 norm{0, 0, 0}; + + norm[t_index()] = 1.; + return norm; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d837c43dec..a690ac2802 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -185,6 +185,7 @@ celeritas_setup_tests(SERIAL PREFIX orange) celeritas_add_test(orange/surfaces/detail/QuadraticSolver.test.cc) celeritas_add_test(orange/surfaces/CylCentered.test.cc) +celeritas_add_test(orange/surfaces/PlaneAligned.test.cc) #-----------------------------------------------------------------------------# # I/O (ROOT) diff --git a/test/orange/surfaces/PlaneAligned.test.cc b/test/orange/surfaces/PlaneAligned.test.cc new file mode 100644 index 0000000000..c9b6c5f082 --- /dev/null +++ b/test/orange/surfaces/PlaneAligned.test.cc @@ -0,0 +1,133 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file PlaneAligned.test.cc +//---------------------------------------------------------------------------// +#include "orange/surfaces/PlaneAligned.hh" + +#include +#include +#include +#include "base/Algorithms.hh" +#include "celeritas_test.hh" + +using celeritas::no_intersection; +using celeritas::PlaneX; +using celeritas::PlaneY; +using celeritas::PlaneZ; +using celeritas::SignedSense; +using celeritas::SurfaceState; + +using celeritas::ipow; +using celeritas::Real3; +using celeritas::real_type; + +using VecReal = std::vector; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class PlaneAlignedTest : public celeritas::Test +{ + protected: + template + real_type calc_intersection(const S& surf, + Real3 pos, + Real3 dir, + SurfaceState s = SurfaceState::off) + { + static_assert(sizeof(typename S::Intersections) == sizeof(double), + "Expected plane to have a single intercept"); + return surf.calc_intersections(pos, dir, s)[0]; + } +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(PlaneAlignedTest, x_plane) +{ + PlaneX p(1.0); + + EXPECT_EQ((Real3{1.0, 0.0, 0.0}), p.calc_normal(Real3{1.0, 0.1, -4.2})); + + // Test sense + EXPECT_EQ(SignedSense::outside, p.calc_sense(Real3{1.01, 0.1, 0.0})); + EXPECT_EQ(SignedSense::inside, p.calc_sense(Real3{0.99, 0.1, 0.0})); + + // Test intersections + Real3 px{1.0, 0.0, 0.0}; + Real3 mx{-1.0, 0.0, 0.0}; + Real3 py{0.0, 1.0, 0.0}; + Real3 pz{0.0, 0.0, 1.0}; + + EXPECT_EQ(no_intersection(), + calc_intersection(p, {0.9999, 0.0, 0.0}, px, SurfaceState::on)); + EXPECT_EQ(no_intersection(), + calc_intersection(p, {1.0001, 0.0, 0.0}, mx, SurfaceState::on)); + EXPECT_SOFT_EQ(0.0, calc_intersection(p, {1.0, 0.0, 0.0}, px)); + EXPECT_SOFT_EQ(0.0, calc_intersection(p, {1.0, 0.0, 0.0}, mx)); + EXPECT_SOFT_EQ(0.01, calc_intersection(p, {0.99, 0.0, 0.0}, px)); + EXPECT_SOFT_EQ(0.01, calc_intersection(p, {1.01, 0.0, 0.0}, mx)); + EXPECT_SOFT_EQ(no_intersection(), + calc_intersection(p, {0.99, 0.0, 0.0}, mx)); + EXPECT_SOFT_EQ(no_intersection(), + calc_intersection(p, {1.01, 0.0, 0.0}, px)); + EXPECT_EQ(no_intersection(), calc_intersection(p, {1.01, 0.0, 0.0}, py)); + EXPECT_EQ(no_intersection(), calc_intersection(p, {0.99, 0.0, 0.0}, pz)); +} + +//---------------------------------------------------------------------------// + +TEST_F(PlaneAlignedTest, y_plane) +{ + PlaneY p(-1.0); + + EXPECT_EQ((Real3{0.0, 1.0, 0.0}), p.calc_normal(Real3{1.0, -1.0, -4.2})); + + // Test sense + EXPECT_EQ(SignedSense::outside, p.calc_sense(Real3{1.01, -0.99, 0.0})); + EXPECT_EQ(SignedSense::inside, p.calc_sense(Real3{0.99, -1.01, 0.0})); + + // Test intersections + Real3 py{0.0, 1.0, 0.0}; + Real3 my{0.0, -1.0, 0.0}; + Real3 px{1.0, 0.0, 0.0}; + + EXPECT_EQ(no_intersection(), + calc_intersection(p, {0, 1 - 1e-8, 0.0}, py, SurfaceState::on)); + EXPECT_SOFT_EQ(0.01, calc_intersection(p, {0.0, -1.01, 0.0}, py)); + EXPECT_SOFT_EQ(no_intersection(), + calc_intersection(p, {0.0, 0.99, 0.0}, my)); + EXPECT_EQ(no_intersection(), calc_intersection(p, {-1.01, 1.0, 0.0}, px)); +} + +//---------------------------------------------------------------------------// + +TEST_F(PlaneAlignedTest, plane_z) +{ + PlaneZ p(0.0); + + EXPECT_EQ((Real3{0.0, 0.0, 1.0}), p.calc_normal(Real3{1.0, 0.1, 0.0})); + + // Test sense + EXPECT_EQ(SignedSense::outside, p.calc_sense(Real3{1.01, 0.1, 0.01})); + EXPECT_EQ(SignedSense::inside, p.calc_sense(Real3{0.99, 0.1, -0.01})); + + // Test intersections + Real3 pz({0.0, 0.0, 1.0}); + Real3 mz({0.0, 0.0, -1.0}); + Real3 px({1.0, 0.0, 0.0}); + + EXPECT_EQ(no_intersection(), + calc_intersection(p, {0.0, 0.0, -1e-8}, pz, SurfaceState::on)); + EXPECT_SOFT_EQ(0.0, calc_intersection(p, {0.0, 0.0, 0.0}, pz)); + EXPECT_SOFT_EQ(0.0, calc_intersection(p, {0.0, 0.0, 0.0}, mz)); + EXPECT_SOFT_EQ(0.01, calc_intersection(p, {0.0, 0.0, -0.01}, pz)); + EXPECT_SOFT_EQ(0.01, calc_intersection(p, {0.0, 0.0, 0.01}, mz)); + EXPECT_EQ(no_intersection(), calc_intersection(p, {-1.01, 0.0, 0.0}, px)); +} From f721f00765264e4370360fc601058c9d8f175ea1 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 19 Oct 2021 17:22:20 -0400 Subject: [PATCH 04/22] Add general quadric and test placeholder --- src/orange/surfaces/CylCentered.hh | 1 - src/orange/surfaces/GeneralQuadric.hh | 213 ++++++++++++++++++++ src/orange/surfaces/PlaneAligned.hh | 1 - test/CMakeLists.txt | 1 + test/orange/surfaces/GeneralQuadric.test.cc | 91 +++++++++ 5 files changed, 305 insertions(+), 2 deletions(-) create mode 100644 src/orange/surfaces/GeneralQuadric.hh create mode 100644 test/orange/surfaces/GeneralQuadric.test.cc diff --git a/src/orange/surfaces/CylCentered.hh b/src/orange/surfaces/CylCentered.hh index 69ba265f9c..d1d41156fe 100644 --- a/src/orange/surfaces/CylCentered.hh +++ b/src/orange/surfaces/CylCentered.hh @@ -9,7 +9,6 @@ #include "base/Array.hh" #include "base/ArrayUtils.hh" -#include "base/Assert.hh" #include "base/Macros.hh" #include "base/Span.hh" #include "../Types.hh" diff --git a/src/orange/surfaces/GeneralQuadric.hh b/src/orange/surfaces/GeneralQuadric.hh new file mode 100644 index 0000000000..b122151205 --- /dev/null +++ b/src/orange/surfaces/GeneralQuadric.hh @@ -0,0 +1,213 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file GeneralQuadric.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Array.hh" +#include "base/ArrayUtils.hh" +#include "base/Span.hh" +#include "base/Types.hh" +#include "detail/QuadraticSolver.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * General quadric surface. + * + * General quadrics that cannot be simplified to other CPTE surfaces include + * hyperboloids and paraboloids; and non-axis-aligned cylinders, ellipsoids, + * and cones. + * + * \f[ + ax^2 + by^2 + cz^2 + dxy + eyz + fzx + gx + hy + iz + j = 0 + \f] + */ +class GeneralQuadric +{ + public: + //@{ + //! Type aliases + using Intersections = Array; + using SpanConstRealN = Span; + using SpanConstReal3 = Span; + //@} + + //// CLASS ATTRIBUTES //// + + //! Surface type identifier + static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type() + { + return SurfaceType::gq; + } + + //! Storage requirements + static CELER_CONSTEXPR_FUNCTION size_type size() + { + return SpanConstRealN::extent; + } + + public: + //// CONSTRUCTORS //// + + // Construct with radius + explicit inline CELER_FUNCTION GeneralQuadric(const Real3& abc, + const Real3& def, + const Real3& ghi, + real_type j); + + // Construct from raw data + explicit inline CELER_FUNCTION GeneralQuadric(SpanConstRealN); + + //// ACCESSORS //// + + //! Second-order terms + CELER_FUNCTION SpanConstReal3 second() const { return {&a_, 3}; } + + //! Cross terms (xy, yz, zx) + CELER_FUNCTION SpanConstReal3 cross() const { return {&d_, 3}; } + + //! First-order terms + CELER_FUNCTION SpanConstReal3 first() const { return {&g_, 3}; } + + //! Zeroth-order term + CELER_FUNCTION real_type zeroth() const { return j_; } + + //! Get a view to the data for type-deleted storage + CELER_FUNCTION SpanConstRealN data() const { return {&a_, 10}; } + + //// CALCULATION //// + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION SignedSense calc_sense(const Real3& pos) const; + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION Intersections calc_intersections( + const Real3& pos, const Real3& dir, SurfaceState on_surface) const; + + // Calculate outward normal at a position + inline CELER_FUNCTION Real3 calc_normal(const Real3& pos) const; + + private: + // Second-order terms (a, b, c) + real_type a_, b_, c_; + // Second-order cross terms (d, e, f) + real_type d_, e_, f_; + // First-order terms (g, h, i) + real_type g_, h_, i_; + // Constant term + real_type j_; +}; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with all coefficients. + */ +CELER_FUNCTION GeneralQuadric::GeneralQuadric(const Real3& abc, + const Real3& def, + const Real3& ghi, + real_type j) + : a_(abc[0]) + , b_(abc[1]) + , c_(abc[2]) + , d_(def[0]) + , e_(def[1]) + , f_(def[2]) + , g_(ghi[0]) + , h_(ghi[1]) + , i_(ghi[2]) + , j_(j) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Construct from raw data. + */ +CELER_FUNCTION GeneralQuadric::GeneralQuadric(SpanConstRealN data) + : a_(data[0]) + , b_(data[1]) + , c_(data[2]) + , d_(data[3]) + , e_(data[4]) + , f_(data[5]) + , g_(data[6]) + , h_(data[7]) + , i_(data[8]) + , j_(data[9]) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Determine the sense of the position relative to this surface. + */ +CELER_FUNCTION SignedSense GeneralQuadric::calc_sense(const Real3& pos) const +{ + const real_type x = pos[0]; + const real_type y = pos[1]; + const real_type z = pos[2]; + + real_type result = (a_ * x + d_ * y + f_ * z + g_) * x + + (b_ * y + e_ * z + h_) * y + (c_ * z + i_) * z + j_; + + return real_to_sense(result); +} + +//---------------------------------------------------------------------------// +/*! + * Determine the sense of the position relative to this surface. + */ +CELER_FUNCTION auto +GeneralQuadric::calc_intersections(const Real3& pos, + const Real3& dir, + SurfaceState on_surface) const + -> Intersections +{ + const real_type x = pos[0]; + const real_type y = pos[1]; + const real_type z = pos[2]; + const real_type u = dir[0]; + const real_type v = dir[1]; + const real_type w = dir[2]; + + // Quadratic values + real_type a = (a_ * u + d_ * v) * u + (b_ * v + e_ * w) * v + + (c_ * w + f_ * u) * w; + real_type half_b = real_type(0.5) + * ((2 * a_ * x + d_ * y + f_ * z + g_) * u + + (2 * b_ * y + d_ * x + e_ * z + h_) * v + + (2 * c_ * z + +e_ * y + f_ * x + i_) * w); + real_type c = ((a_ * x + d_ * y + g_) * x + (b_ * y + e_ * z + h_) * y + + (c_ * z + f_ * x + i_) * z + j_); + + return detail::QuadraticSolver::solve_general(a, half_b, c, on_surface); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate outward normal at a position. + */ +CELER_FUNCTION Real3 GeneralQuadric::calc_normal(const Real3& pos) const +{ + const real_type x = pos[0]; + const real_type y = pos[1]; + const real_type z = pos[2]; + + Real3 norm; + norm[0] = 2 * a_ * x + d_ * y + f_ * z + g_; + norm[1] = 2 * b_ * y + d_ * x + e_ * z + h_; + norm[2] = 2 * c_ * z + e_ * y + f_ * x + i_; + + normalize_direction(&norm); + return norm; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/surfaces/PlaneAligned.hh b/src/orange/surfaces/PlaneAligned.hh index 64cdb957f6..f85bd5c8bd 100644 --- a/src/orange/surfaces/PlaneAligned.hh +++ b/src/orange/surfaces/PlaneAligned.hh @@ -8,7 +8,6 @@ #pragma once #include "base/Array.hh" -#include "base/Assert.hh" #include "base/Macros.hh" #include "base/Span.hh" #include "../Types.hh" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a690ac2802..6f15b0b12e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -185,6 +185,7 @@ celeritas_setup_tests(SERIAL PREFIX orange) celeritas_add_test(orange/surfaces/detail/QuadraticSolver.test.cc) celeritas_add_test(orange/surfaces/CylCentered.test.cc) +celeritas_add_test(orange/surfaces/GeneralQuadric.test.cc) celeritas_add_test(orange/surfaces/PlaneAligned.test.cc) #-----------------------------------------------------------------------------# diff --git a/test/orange/surfaces/GeneralQuadric.test.cc b/test/orange/surfaces/GeneralQuadric.test.cc new file mode 100644 index 0000000000..a7dfab905c --- /dev/null +++ b/test/orange/surfaces/GeneralQuadric.test.cc @@ -0,0 +1,91 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file GeneralQuadric.test.cc +//---------------------------------------------------------------------------// +#include "orange/surfaces/GeneralQuadric.hh" + +#include "celeritas_test.hh" + +using celeritas::GeneralQuadric; +using celeritas::no_intersection; +using celeritas::SignedSense; +using celeritas::SurfaceState; + +using celeritas::ipow; +using celeritas::Real3; +using celeritas::real_type; + +using Intersections = GeneralQuadric::Intersections; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class GeneralQuadricTest : public celeritas::Test +{ + protected: + void SetUp() override {} +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// +/*! + * This shape is an ellipsoid: + * - with radii {3, 2, 1}, + * - rotated by 60 degrees about the x axis, + * - then 30 degrees about z, + * - and finally translated by {1, 2, 3}. + */ +TEST_F(GeneralQuadricTest, all) +{ + const Real3 second{10.3125, 22.9375, 15.75}; + const Real3 cross{-21.867141445557, -20.25, 11.69134295109}; + const Real3 first{-11.964745962156, -9.1328585544429, -65.69134295109}; + real_type zeroth = 77.652245962156; + + GeneralQuadric gq{second, cross, first, zeroth}; + EXPECT_VEC_SOFT_EQ(second, gq.second()); + EXPECT_VEC_SOFT_EQ(cross, gq.cross()); + EXPECT_VEC_SOFT_EQ(first, gq.first()); + EXPECT_SOFT_EQ(zeroth, gq.zeroth()); + + EXPECT_EQ(SignedSense::outside, gq.calc_sense({4, 5, 5})); + EXPECT_EQ(SignedSense::inside, gq.calc_sense({1, 2, 3})); + + const Real3 center{1, 2, 3}; + const Real3 on_surface{3.598076211353292, 3.5, 3}; + const Real3 inward{-0.8660254037844386, -0.5, 0}; + const Real3 outward{0.8660254037844386, 0.5, 0}; + + Intersections distances; + + // On surface, inward + distances = gq.calc_intersections(on_surface, inward, SurfaceState::on); + EXPECT_SOFT_EQ(6.0, distances[0]); + EXPECT_SOFT_EQ(no_intersection(), distances[1]); + + // "Not on surface", inward + distances = gq.calc_intersections(on_surface, inward, SurfaceState::off); + EXPECT_SOFT_EQ(1e-16, distances[0]); + EXPECT_SOFT_EQ(6.0, distances[1]); + + // On surface, outward + distances = gq.calc_intersections(on_surface, outward, SurfaceState::on); + EXPECT_SOFT_EQ(no_intersection(), distances[0]); + EXPECT_SOFT_EQ(no_intersection(), distances[1]); + + // In center + distances = gq.calc_intersections(center, outward, SurfaceState::off); + EXPECT_SOFT_EQ(3.0, distances[1]); + EXPECT_SOFT_EQ(no_intersection(), distances[0]); + + // Outside, hitting both + const Real3 pos{-2.464101615137754, 0, 3}; + distances = gq.calc_intersections(pos, outward, SurfaceState::off); + EXPECT_SOFT_EQ(1.0, distances[0]); + EXPECT_SOFT_EQ(7.0, distances[1]); +} From 0c2a01f44f5f4cf6df0995a9b29ca24c17a8e420 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Wed, 20 Oct 2021 09:48:20 -0400 Subject: [PATCH 05/22] Add stream printing to surfaces --- src/CMakeLists.txt | 1 + src/orange/surfaces/SurfaceIO.cc | 53 ++++++++++++++++++++ src/orange/surfaces/SurfaceIO.hh | 28 +++++++++++ src/orange/surfaces/SurfaceTypeTraits.hh | 62 ++++++++++++++++++++++++ 4 files changed, 144 insertions(+) create mode 100644 src/orange/surfaces/SurfaceIO.cc create mode 100644 src/orange/surfaces/SurfaceIO.hh create mode 100644 src/orange/surfaces/SurfaceTypeTraits.hh diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5d839902b1..2114fb8049 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -62,6 +62,7 @@ list(APPEND SOURCES comm/detail/LoggerMessage.cc geometry/detail/ScopedTimeAndRedirect.cc orange/Types.cc + orange/surfaces/SurfaceIO.cc io/ImportProcess.cc io/ImportPhysicsTable.cc io/ImportPhysicsVector.cc diff --git a/src/orange/surfaces/SurfaceIO.cc b/src/orange/surfaces/SurfaceIO.cc new file mode 100644 index 0000000000..542954019d --- /dev/null +++ b/src/orange/surfaces/SurfaceIO.cc @@ -0,0 +1,53 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceIO.cc +//---------------------------------------------------------------------------// +#include "SurfaceIO.hh" + +#include +#include "base/SpanIO.hh" +#include "CylCentered.hh" +#include "GeneralQuadric.hh" +#include "PlaneAligned.hh" + +namespace celeritas +{ +#define ORANGE_INSTANTIATE_SHAPE_STREAM(SHAPE) \ + template std::ostream& operator<<(std::ostream&, const SHAPE&); \ + template std::ostream& operator<<(std::ostream&, const SHAPE&); \ + template std::ostream& operator<<(std::ostream&, const SHAPE&) + +//---------------------------------------------------------------------------// +template +std::ostream& operator<<(std::ostream& os, const CylCentered& s) +{ + os << "Cyl " << to_char(T) << ": r=" << std::sqrt(s.radius_sq()); + return os; +} + +ORANGE_INSTANTIATE_SHAPE_STREAM(CylCentered); +//---------------------------------------------------------------------------// +std::ostream& operator<<(std::ostream& os, const GeneralQuadric& s) +{ + os << "GQuadric: " << s.second() << ' ' << s.cross() << ' ' << s.first() + << ' ' << s.zeroth(); + + return os; +} + +//---------------------------------------------------------------------------// +template +std::ostream& operator<<(std::ostream& os, const PlaneAligned& s) +{ + os << "Plane: " << to_char(T) << '=' << s.position(); + return os; +} + +ORANGE_INSTANTIATE_SHAPE_STREAM(PlaneAligned); +//---------------------------------------------------------------------------// + +#undef ORANGE_INSTANTIATE_SHAPE_STREAM +} // namespace celeritas diff --git a/src/orange/surfaces/SurfaceIO.hh b/src/orange/surfaces/SurfaceIO.hh new file mode 100644 index 0000000000..4be6c0cacc --- /dev/null +++ b/src/orange/surfaces/SurfaceIO.hh @@ -0,0 +1,28 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceIO.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include "SurfaceTypeTraits.hh" +#include "../Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +//!@{ +//! Print surfaces to a stream. +template +std::ostream& operator<<(std::ostream& os, const PlaneAligned& s); + +template +std::ostream& operator<<(std::ostream& os, const CylCentered& s); + +std::ostream& operator<<(std::ostream& os, const GeneralQuadric& s); +//!@} +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/surfaces/SurfaceTypeTraits.hh b/src/orange/surfaces/SurfaceTypeTraits.hh new file mode 100644 index 0000000000..e5ad7d3fb8 --- /dev/null +++ b/src/orange/surfaces/SurfaceTypeTraits.hh @@ -0,0 +1,62 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceTypeTraits.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "../Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +// FORWARD DECLARATIONS +//---------------------------------------------------------------------------// +template +class PlaneAligned; +template +class CylCentered; +class GeneralQuadric; + +//---------------------------------------------------------------------------// +/*! + * Map surface enumeration to surface type. + */ +template +struct SurfaceTypeTraits; + +#define ORANGE_SURFACE_TRAITS(ENUM_VALUE, CLS) \ + template<> \ + struct SurfaceTypeTraits \ + { \ + using type = CLS; \ + } + +// clang-format off +ORANGE_SURFACE_TRAITS(px, PlaneAligned); +ORANGE_SURFACE_TRAITS(py, PlaneAligned); +ORANGE_SURFACE_TRAITS(pz, PlaneAligned); +ORANGE_SURFACE_TRAITS(cxc, CylCentered); +ORANGE_SURFACE_TRAITS(cyc, CylCentered); +ORANGE_SURFACE_TRAITS(czc, CylCentered); +#if 0 +ORANGE_SURFACE_TRAITS(sc, CenteredSphere); +ORANGE_SURFACE_TRAITS(cx, CylAligned); +ORANGE_SURFACE_TRAITS(cy, CylAligned); +ORANGE_SURFACE_TRAITS(cz, CylAligned); +ORANGE_SURFACE_TRAITS(p, Plane); +ORANGE_SURFACE_TRAITS(s, Sphere); +ORANGE_SURFACE_TRAITS(kx, ConeAligned); +ORANGE_SURFACE_TRAITS(ky, ConeAligned); +ORANGE_SURFACE_TRAITS(kz, ConeAligned); +ORANGE_SURFACE_TRAITS(sq, SimpleQuadric); +#endif +ORANGE_SURFACE_TRAITS(gq, GeneralQuadric); +// clang-format on + +#undef ORANGE_SURFACE_TRAITS + +//---------------------------------------------------------------------------// +} // namespace celeritas From a67045fb75253f49c372b9c5b6a0be89f2cdc520 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Wed, 20 Oct 2021 21:06:39 -0400 Subject: [PATCH 06/22] Add surface inserter, data interface, and surface action --- src/CMakeLists.txt | 1 + src/base/Collection.hh | 5 +- src/base/Span.hh | 5 +- src/geometry/GeoData.hh | 2 +- src/orange/Data.hh | 210 ++++++++++++++++++++ src/orange/Types.cc | 37 ++++ src/orange/Types.hh | 12 +- src/orange/construct/SurfaceInserter.cc | 49 +++++ src/orange/construct/SurfaceInserter.hh | 84 ++++++++ src/orange/surfaces/CylCentered.hh | 15 +- src/orange/surfaces/GeneralQuadric.hh | 23 +-- src/orange/surfaces/PlaneAligned.hh | 19 +- src/orange/surfaces/SurfaceAction.hh | 28 +++ src/orange/surfaces/Surfaces.hh | 94 +++++++++ src/orange/surfaces/detail/SurfaceAction.hh | 131 ++++++++++++ test/CMakeLists.txt | 1 + test/orange/surfaces/CylCentered.test.cc | 2 +- test/orange/surfaces/SurfaceAction.test.cc | 101 ++++++++++ 18 files changed, 769 insertions(+), 50 deletions(-) create mode 100644 src/orange/Data.hh create mode 100644 src/orange/construct/SurfaceInserter.cc create mode 100644 src/orange/construct/SurfaceInserter.hh create mode 100644 src/orange/surfaces/SurfaceAction.hh create mode 100644 src/orange/surfaces/Surfaces.hh create mode 100644 src/orange/surfaces/detail/SurfaceAction.hh create mode 100644 test/orange/surfaces/SurfaceAction.test.cc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2114fb8049..c1b507e4aa 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -62,6 +62,7 @@ list(APPEND SOURCES comm/detail/LoggerMessage.cc geometry/detail/ScopedTimeAndRedirect.cc orange/Types.cc + orange/construct/SurfaceInserter.cc orange/surfaces/SurfaceIO.cc io/ImportProcess.cc io/ImportPhysicsTable.cc diff --git a/src/base/Collection.hh b/src/base/Collection.hh index 9f71aececf..8896a132f7 100644 --- a/src/base/Collection.hh +++ b/src/base/Collection.hh @@ -41,8 +41,9 @@ namespace celeritas * data collections are typically accessed by thread ID. \c ParamsData are * immutable and always "mirrored" on both host and device. Sometimes it's * sensible to partition \c ParamsData into discrete helper structs (stored by - * value), each with a group of collections: each of these should be called a - * \c DataGroup. + * value), each with a group of collections, and perhaps another struct that + * has non-templated scalars (since the default assignment operator is less + * work than manually copying scalars in a templated assignment operator. * * A collection group has the following requirements to be compatible with the \c diff --git a/src/base/Span.hh b/src/base/Span.hh index ad017489d5..3eda5a3b89 100644 --- a/src/base/Span.hh +++ b/src/base/Span.hh @@ -51,7 +51,7 @@ class Span //!@} //! Size (may be dynamic) - static constexpr std::size_t extent = Extent; + static constexpr CELER_FUNCTION std::size_t extent = Extent; public: //// CONSTRUCTION //// @@ -145,6 +145,9 @@ class Span detail::SpanImpl s_; }; +template +constexpr CELER_FUNCTION std::size_t Span::extent; + //---------------------------------------------------------------------------// // HELPER FUNCTIONS //---------------------------------------------------------------------------// diff --git a/src/geometry/GeoData.hh b/src/geometry/GeoData.hh index 53c4527287..8734a8e097 100644 --- a/src/geometry/GeoData.hh +++ b/src/geometry/GeoData.hh @@ -88,7 +88,7 @@ struct GeoStateData //// METHODS //// - //! True if assigned + //! True if sizes are consistent and states are assigned explicit CELER_FUNCTION operator bool() const { return this->size() > 0 && dir.size() == this->size() diff --git a/src/orange/Data.hh b/src/orange/Data.hh new file mode 100644 index 0000000000..ccbd814438 --- /dev/null +++ b/src/orange/Data.hh @@ -0,0 +1,210 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file Data.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Collection.hh" +#include "base/CollectionBuilder.hh" +#include "base/OpaqueId.hh" +#include "Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +// PARAMS +//---------------------------------------------------------------------------// +/*! + * Data for surface definitions. + * + * Surfaces each have a compile-time number of real data needed to define them. + * (These usually are the nonzero coefficients of the quadric equation.) A + * surface ID points to an offset into the `data` field. These surface IDs are + * *global* over all universes. + */ +template +struct SurfaceData +{ + //// TYPES //// + + template + using Items = Collection; + + //// DATA //// + + Items types; + Items> offsets; + Collection reals; + + //// METHODS //// + + //! True if sizes are valid + explicit CELER_FUNCTION operator bool() const + { + return !types.empty() && offsets.size() == types.size() + && reals.size() >= types.size(); + } + + //! Assign from another set of data + template + SurfaceData& operator=(const SurfaceData& other) + { + CELER_EXPECT(other); + types = other.types; + offsets = other.offsets; + reals = other.reals; + return *this; + } +}; + +//---------------------------------------------------------------------------// +/*! + * Data for a single volume definition. + * + * A volume is a CSG tree of surfaces. Each surface defines an "inside" space + * and "outside" space that correspond to "negative" and "positive" values of + * the quadric expression's evaluation. Left of a plane is negative, for + * example, and evaluates to "false" or "inside" or "negative". The CSG tree is + * flattened into a + */ +struct VolumeDef +{ + ItemRange faces; + ItemRange logic; + + logic_int num_intersections{0}; + logic_int flags{0}; +}; + +//---------------------------------------------------------------------------// +/*! + * Data for volume definitions. + */ +template +struct VolumeData +{ + //// TYPES //// + + template + using Items = Collection; + + //// DATA //// + + Items defs; + + // Storage + Collection faces; + Collection logic; + + //// METHODS //// + + //! True if sizes are valid + explicit CELER_FUNCTION operator bool() const { return defs.size() > 0; } + + //! Assign from another set of data + template + VolumeData& operator=(const VolumeData& other) + { + CELER_EXPECT(other); + + defs = other.defs; + faces = other.faces; + logic = other.logic; + + return *this; + } +}; + +//---------------------------------------------------------------------------// +/*! + * Data for universe definitions. + */ +template +struct UniverseData +{ + //// TYPES //// + + template + using Items = Collection; + + //// DATA //// + + //// METHODS //// + + //! True if sizes are valid + explicit CELER_FUNCTION operator bool() const { return false; } + + //! Assign from another set of data + template + UniverseData& operator=(const UniverseData& other) + { + CELER_EXPECT(other); + return *this; + } +}; + +//---------------------------------------------------------------------------// +/*! + * Scalar values particular to an ORANGE geometry instance. + */ +struct OrangeParamsScalars +{ + size_type max_level{}; + size_type max_faces{}; + size_type max_intersections{}; + + // TODO: fuzziness/length scale +}; + +//---------------------------------------------------------------------------// +/*! + * Data to persistent data used by ORANGE implementation. + */ +template +struct OrangeParamsData +{ + //// DATA //// + + SurfaceData surfaces; + VolumeData volumes; + UniverseData universes; + + OrangeParamsScalars scalars; + + //// METHODS //// + + //! True if assigned + explicit CELER_FUNCTION operator bool() const + { + return surfaces && volumes && universes; + } + + //! Assign from another set of data + template + OrangeParamsData& operator=(const OrangeParamsData& other) + { + CELER_EXPECT(other); + surfaces = other.surfaces; + volumes = other.volumes; + universes = other.universes; + scalars = other.scalars; + return *this; + } +}; + +//---------------------------------------------------------------------------// +// STATE +//---------------------------------------------------------------------------// +/*! + * Data to persistent data used by ORANGE implementation. + */ +template +struct OrangeStateData +{ +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/Types.cc b/src/orange/Types.cc index ef172b9b3a..c4a7372635 100644 --- a/src/orange/Types.cc +++ b/src/orange/Types.cc @@ -7,9 +7,46 @@ //---------------------------------------------------------------------------// #include "Types.hh" +#include "base/Assert.hh" + namespace celeritas { //---------------------------------------------------------------------------// +/*! + * Get a string corresponding to a surface type. + */ +const char* to_cstring(SurfaceType value) +{ + CELER_EXPECT(value != SurfaceType::size_); + + static const char* const strings[] = { + "px", + "py", + "pz", + "cxc", + "cyc", + "czc", +#if 0 + "sc", + "cx", + "cy", + "cz", + "p", + "s", + "kx", + "ky", + "kz", + "sq", +#endif + "gq", + }; + static_assert( + static_cast(SurfaceType::size_) * sizeof(const char*) + == sizeof(strings), + "Enum strings are incorrect"); + + return strings[static_cast(value)]; +} //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/orange/Types.hh b/src/orange/Types.hh index 95c2fa69ec..a13cdf4ebc 100644 --- a/src/orange/Types.hh +++ b/src/orange/Types.hh @@ -34,14 +34,13 @@ using UniverseId = OpaqueId; /*! * Whether a position is logically "inside" or "outside" a surface. * - * For a plane, "POS" (outside/true) is equivalent to + * For a plane, "outside" (true) is the "positive" sense and equivalent to * \f[ \vec x \cdot \vec n >= 0 * \f] - * and "inside" (NEG) is to the left of the plane's normal. Likewise, for a + * and "inside" is to the left of the plane's normal. Likewise, for a * sphere, "inside" is where the dot product of the position and outward normal - * is negative. These are *opposite* the signs from KENO, where `-` has the - * sense of negating a closed region of space (i.e. the interior of a sphere). + * is negative. */ enum class Sense : bool { @@ -132,8 +131,7 @@ enum class SurfaceState : bool // CLASSES //---------------------------------------------------------------------------// /*! - * \struct Initialization - * \brief Volume ID and surface ID after initialization + * Volume ID and surface ID after initialization. * * Possible configurations for the initialization result ('X' means 'has * a valid ID', i.e. evaluates to true): @@ -217,6 +215,8 @@ CELER_CONSTEXPR_FUNCTION Sense to_sense(SignedSense s) * Sentinel value indicating "no intersection". * * \todo There is probably a better place to put this since it's not a "type". + * \todo A value of zero might also work since zero-length steps are + * prohibited. */ CELER_CONSTEXPR_FUNCTION real_type no_intersection() { diff --git a/src/orange/construct/SurfaceInserter.cc b/src/orange/construct/SurfaceInserter.cc new file mode 100644 index 0000000000..76da09bf9b --- /dev/null +++ b/src/orange/construct/SurfaceInserter.cc @@ -0,0 +1,49 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceInserter.cc +//---------------------------------------------------------------------------// +#include "SurfaceInserter.hh" + +#include "base/CollectionBuilder.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct with a reference to empty surfaces. + */ +SurfaceInserter::SurfaceInserter(Data* surfaces) : surface_data_(surfaces) +{ + CELER_EXPECT(surface_data_ && surface_data_->types.empty() + && surface_data_->offsets.empty() + && surface_data_->reals.empty()); +} + +//---------------------------------------------------------------------------// +/*! + * Insert a generic surface. + */ +SurfaceId SurfaceInserter::operator()(GenericSurfaceRef generic_surf) +{ + CELER_EXPECT(generic_surf); + + // TODO: surface deduplication goes here + + auto types = make_builder(&surface_data_->types); + auto offsets = make_builder(&surface_data_->offsets); + auto reals = make_builder(&surface_data_->reals); + + SurfaceId::size_type new_id = types.size(); + types.push_back(generic_surf.type); + offsets.push_back(OpaqueId(reals.size())); + reals.insert_back(generic_surf.data.begin(), generic_surf.data.end()); + + CELER_ENSURE(types.size() == offsets.size()); + return SurfaceId{new_id}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/construct/SurfaceInserter.hh b/src/orange/construct/SurfaceInserter.hh new file mode 100644 index 0000000000..9daab7945a --- /dev/null +++ b/src/orange/construct/SurfaceInserter.hh @@ -0,0 +1,84 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceInserter.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "../Data.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct surfaces on the host. + * + * Currently this simply appends the surface to the Data, but the full + * robust geometry implementation will implement "soft" surface deduplication. + * + * \code + SurfaceInserter insert_surface(¶ms.surfaces); + auto id = insert_surface(PlaneX(123)); + auto id2 = insert_surface(PlaneX(123.0000001)); // TODO: equals id + \endcode + */ +class SurfaceInserter +{ + public: + //!@{ + //! Type aliases + using Data = SurfaceData; + //!@} + + //! Type-deleted reference to a surface + struct GenericSurfaceRef + { + SurfaceType type; + Span data; + + inline explicit operator bool() const; + }; + + public: + // Construct with reference to surfaces to build + explicit SurfaceInserter(Data* surfaces_); + + // Add a new surface + template + inline SurfaceId operator()(const T& surface); + + // Append a generic surface view to the vector + SurfaceId operator()(GenericSurfaceRef generic_surf); + + private: + Data* surface_data_; +}; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Whether a generic surface reference is valid. + */ +SurfaceInserter::GenericSurfaceRef::operator bool() const +{ + return !data.empty() && type != SurfaceType::size_; +} + +//---------------------------------------------------------------------------// +/*! + * Add a surface (type-deleted) with the given coefficients + */ +template +SurfaceId SurfaceInserter::operator()(const T& surface) +{ + static_assert(sizeof(typename T::Intersections) > 0, + "Template parameter must be a surface class"); + + return (*this)(GenericSurfaceRef{surface.surface_type(), surface.data()}); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/surfaces/CylCentered.hh b/src/orange/surfaces/CylCentered.hh index d1d41156fe..919b41e09b 100644 --- a/src/orange/surfaces/CylCentered.hh +++ b/src/orange/surfaces/CylCentered.hh @@ -41,7 +41,7 @@ class CylCentered //@{ //! Type aliases using Intersections = Array; - using SpanConstRealN = Span; + using Storage = Span; //@} //// CLASS ATTRIBUTES //// @@ -49,12 +49,6 @@ class CylCentered // Surface type identifier static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type(); - //! Storage requirements - static CELER_CONSTEXPR_FUNCTION size_type size() - { - return SpanConstRealN::extent; - } - public: //// CONSTRUCTORS //// @@ -62,7 +56,7 @@ class CylCentered explicit inline CELER_FUNCTION CylCentered(real_type radius); // Construct from raw data - explicit inline CELER_FUNCTION CylCentered(SpanConstRealN); + explicit inline CELER_FUNCTION CylCentered(Storage); //// ACCESSORS //// @@ -70,7 +64,7 @@ class CylCentered CELER_FUNCTION real_type radius_sq() const { return radius_sq_; } //! Get a view to the data for type-deleted storage - CELER_FUNCTION SpanConstRealN data() const { return {&radius_sq_, 1}; } + CELER_FUNCTION Storage data() const { return {&radius_sq_, 1}; } //// CALCULATION //// @@ -131,8 +125,7 @@ CELER_FUNCTION CylCentered::CylCentered(real_type radius) * Construct from raw data. */ template -CELER_FUNCTION CylCentered::CylCentered(SpanConstRealN data) - : radius_sq_(data[0]) +CELER_FUNCTION CylCentered::CylCentered(Storage data) : radius_sq_(data[0]) { } diff --git a/src/orange/surfaces/GeneralQuadric.hh b/src/orange/surfaces/GeneralQuadric.hh index b122151205..416d58e2cb 100644 --- a/src/orange/surfaces/GeneralQuadric.hh +++ b/src/orange/surfaces/GeneralQuadric.hh @@ -33,7 +33,7 @@ class GeneralQuadric //@{ //! Type aliases using Intersections = Array; - using SpanConstRealN = Span; + using Storage = Span; using SpanConstReal3 = Span; //@} @@ -45,12 +45,6 @@ class GeneralQuadric return SurfaceType::gq; } - //! Storage requirements - static CELER_CONSTEXPR_FUNCTION size_type size() - { - return SpanConstRealN::extent; - } - public: //// CONSTRUCTORS //// @@ -61,7 +55,7 @@ class GeneralQuadric real_type j); // Construct from raw data - explicit inline CELER_FUNCTION GeneralQuadric(SpanConstRealN); + explicit inline CELER_FUNCTION GeneralQuadric(Storage); //// ACCESSORS //// @@ -78,7 +72,7 @@ class GeneralQuadric CELER_FUNCTION real_type zeroth() const { return j_; } //! Get a view to the data for type-deleted storage - CELER_FUNCTION SpanConstRealN data() const { return {&a_, 10}; } + CELER_FUNCTION Storage data() const { return {&a_, 10}; } //// CALCULATION //// @@ -130,7 +124,7 @@ CELER_FUNCTION GeneralQuadric::GeneralQuadric(const Real3& abc, /*! * Construct from raw data. */ -CELER_FUNCTION GeneralQuadric::GeneralQuadric(SpanConstRealN data) +CELER_FUNCTION GeneralQuadric::GeneralQuadric(Storage data) : a_(data[0]) , b_(data[1]) , c_(data[2]) @@ -180,14 +174,13 @@ GeneralQuadric::calc_intersections(const Real3& pos, // Quadratic values real_type a = (a_ * u + d_ * v) * u + (b_ * v + e_ * w) * v + (c_ * w + f_ * u) * w; - real_type half_b = real_type(0.5) - * ((2 * a_ * x + d_ * y + f_ * z + g_) * u - + (2 * b_ * y + d_ * x + e_ * z + h_) * v - + (2 * c_ * z + +e_ * y + f_ * x + i_) * w); + real_type b = (2 * a_ * x + d_ * y + f_ * z + g_) * u + + (2 * b_ * y + d_ * x + e_ * z + h_) * v + + (2 * c_ * z + e_ * y + f_ * x + i_) * w; real_type c = ((a_ * x + d_ * y + g_) * x + (b_ * y + e_ * z + h_) * y + (c_ * z + f_ * x + i_) * z + j_); - return detail::QuadraticSolver::solve_general(a, half_b, c, on_surface); + return detail::QuadraticSolver::solve_general(a, b / 2, c, on_surface); } //---------------------------------------------------------------------------// diff --git a/src/orange/surfaces/PlaneAligned.hh b/src/orange/surfaces/PlaneAligned.hh index f85bd5c8bd..1de631defc 100644 --- a/src/orange/surfaces/PlaneAligned.hh +++ b/src/orange/surfaces/PlaneAligned.hh @@ -25,7 +25,7 @@ class PlaneAligned //@{ //! Type aliases using Intersections = Array; - using SpanConstRealN = Span; + using Storage = Span; //@} //// CLASS ATTRIBUTES //// @@ -33,12 +33,6 @@ class PlaneAligned // Surface type identifier static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type(); - //! Storage requirements - static CELER_CONSTEXPR_FUNCTION size_type size() - { - return SpanConstRealN::extent; - } - public: //// CONSTRUCTORS //// @@ -46,7 +40,7 @@ class PlaneAligned explicit inline CELER_FUNCTION PlaneAligned(real_type position); // Construct from raw data - explicit inline CELER_FUNCTION PlaneAligned(SpanConstRealN); + explicit inline CELER_FUNCTION PlaneAligned(Storage); //// ACCESSORS //// @@ -54,7 +48,7 @@ class PlaneAligned CELER_FUNCTION real_type position() const { return position_; } //! Get a view to the data for type-deleted storage - CELER_FUNCTION SpanConstRealN data() const { return {&position_, 1}; } + CELER_FUNCTION Storage data() const { return {&position_, 1}; } //// CALCULATION //// @@ -101,7 +95,7 @@ CELER_CONSTEXPR_FUNCTION SurfaceType PlaneAligned::surface_type() //---------------------------------------------------------------------------// /*! - * Construct with radius. + * Construct from axis intercept. */ template CELER_FUNCTION PlaneAligned::PlaneAligned(real_type position) @@ -114,8 +108,7 @@ CELER_FUNCTION PlaneAligned::PlaneAligned(real_type position) * Construct from raw data. */ template -CELER_FUNCTION PlaneAligned::PlaneAligned(SpanConstRealN data) - : position_(data[0]) +CELER_FUNCTION PlaneAligned::PlaneAligned(Storage data) : position_(data[0]) { } @@ -156,7 +149,7 @@ PlaneAligned::calc_intersections(const Real3& pos, * Calculate outward normal at a position. */ template -CELER_FUNCTION Real3 PlaneAligned::calc_normal(const Real3& pos) const +CELER_FUNCTION Real3 PlaneAligned::calc_normal(const Real3&) const { Real3 norm{0, 0, 0}; diff --git a/src/orange/surfaces/SurfaceAction.hh b/src/orange/surfaces/SurfaceAction.hh new file mode 100644 index 0000000000..09d2f874c8 --- /dev/null +++ b/src/orange/surfaces/SurfaceAction.hh @@ -0,0 +1,28 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceAction.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Macros.hh" +#include "Surfaces.hh" +#include "detail/SurfaceAction.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Helper function for creating a SurfaceAction instance. + */ +template +inline CELER_FUNCTION detail::SurfaceAction + make_surface_action(const Surfaces& surfaces, F&& action) +{ + return detail::SurfaceAction{surfaces, std::forward(action)}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/surfaces/Surfaces.hh b/src/orange/surfaces/Surfaces.hh new file mode 100644 index 0000000000..2b5a24a6a0 --- /dev/null +++ b/src/orange/surfaces/Surfaces.hh @@ -0,0 +1,94 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file Surfaces.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "../Data.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Operations on one or more surfaces. + * + * These are all surfaces in the entire geometry. + */ +class Surfaces +{ + public: + //@{ + //! Type aliases + using Data = SurfaceData; + //@} + + public: + // Construct with defaults + explicit inline CELER_FUNCTION Surfaces(const Data&); + + // Number of surfaces + inline CELER_FUNCTION SurfaceId::size_type num_surfaces() const; + + // Get the type of an indexed surface + inline CELER_FUNCTION SurfaceType surface_type(SurfaceId) const; + + // Convert a stored surface to a class + template + inline CELER_FUNCTION T make_surface(SurfaceId) const; + + private: + const Data& data_; +}; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with defaults. + */ +CELER_FUNCTION Surfaces::Surfaces(const Data& data) : data_(data) +{ + CELER_EXPECT(data_); +} + +//---------------------------------------------------------------------------// +/*! + * Number of surfaces. + */ +CELER_FUNCTION SurfaceId::size_type Surfaces::num_surfaces() const +{ + return data_.types.size(); +} + +//---------------------------------------------------------------------------// +/*! + * Get the type of an indexed surface. + */ +CELER_FUNCTION SurfaceType Surfaces::surface_type(SurfaceId sid) const +{ + CELER_EXPECT(sid < this->num_surfaces()); + return data_.types[sid]; +} + +//---------------------------------------------------------------------------// +/*! + * Convert a stored surface to a class. + */ +template +CELER_FUNCTION T Surfaces::make_surface(SurfaceId sid) const +{ + static_assert(T::Storage::extent > 0, + "Template parameter must be a surface class"); + CELER_EXPECT(sid < this->num_surfaces()); + OpaqueId start = data_.offsets[sid]; + OpaqueId end{start.get() + + static_cast(T::Storage::extent)}; + CELER_ASSERT(end.unchecked_get() <= data_.reals.size()); + return T{data_.reals[ItemRange{start, end}]}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/surfaces/detail/SurfaceAction.hh b/src/orange/surfaces/detail/SurfaceAction.hh new file mode 100644 index 0000000000..5a785e04c7 --- /dev/null +++ b/src/orange/surfaces/detail/SurfaceAction.hh @@ -0,0 +1,131 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceAction.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +#include "base/Macros.hh" +#include "orange/Types.hh" +#include "../PlaneAligned.hh" +#include "../CylCentered.hh" +#include "../GeneralQuadric.hh" +#include "../SurfaceTypeTraits.hh" +#include "../Surfaces.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Helper class for applying an action functor to a generic surface. + * + * The templated operator() of the given functor F must be a Surface class. The + * `result_type` type alias here uses GeneralQuadric to represent the "most + * generic" type the functor accepts. + */ +template +class SurfaceAction +{ + public: + //@{ + //! Type aliases + using result_type + = decltype(std::declval()(std::declval())); + //@} + + public: + // Construct from surfaces and action + inline CELER_FUNCTION SurfaceAction(const Surfaces& surfaces, F&& action); + + // Apply to the surface specified by a surface ID + inline CELER_FUNCTION result_type operator()(SurfaceId id); + + //! Access the resulting action + CELER_FUNCTION const F& action() const { return action_; } + + private: + //// DATA //// + Surfaces surfaces_; + F action_; + + //// METHODS //// + template + inline CELER_FUNCTION result_type apply_impl(SurfaceId id); +}; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with reference to surfaces and action to apply. + */ +template +CELER_FUNCTION +SurfaceAction::SurfaceAction(const Surfaces& surfaces, F&& action) + : surfaces_(surfaces), action_(std::forward(action)) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Apply to the surface specified by the given surface ID. + */ +template +CELER_FUNCTION auto SurfaceAction::operator()(SurfaceId id) -> result_type +{ + CELER_EXPECT(id < surfaces_.num_surfaces()); +#define ORANGE_SURF_APPLY_IMPL(TYPE) \ + case (SurfaceType::TYPE): \ + return this->apply_impl(id) + + switch (surfaces_.surface_type(id)) + { + ORANGE_SURF_APPLY_IMPL(px); + ORANGE_SURF_APPLY_IMPL(py); + ORANGE_SURF_APPLY_IMPL(pz); + ORANGE_SURF_APPLY_IMPL(cxc); + ORANGE_SURF_APPLY_IMPL(cyc); + ORANGE_SURF_APPLY_IMPL(czc); +#if 0 + ORANGE_SURF_APPLY_IMPL(sc); + ORANGE_SURF_APPLY_IMPL(cx); + ORANGE_SURF_APPLY_IMPL(cy); + ORANGE_SURF_APPLY_IMPL(cz); + ORANGE_SURF_APPLY_IMPL(p); + ORANGE_SURF_APPLY_IMPL(s); + ORANGE_SURF_APPLY_IMPL(kx); + ORANGE_SURF_APPLY_IMPL(ky); + ORANGE_SURF_APPLY_IMPL(kz); + ORANGE_SURF_APPLY_IMPL(sq); +#endif + ORANGE_SURF_APPLY_IMPL(gq); + case SurfaceType::size_: + CELER_ASSERT_UNREACHABLE(); + } +#undef ORANGE_SURF_APPLY_IMPL + CELER_ASSERT_UNREACHABLE(); +} + +//---------------------------------------------------------------------------// +// PRIVATE INLINE FUNCTIONS +//---------------------------------------------------------------------------// +/*! + * Apply to the surface specified by a surface ID. + */ +template +template +CELER_FUNCTION auto SurfaceAction::apply_impl(SurfaceId id) -> result_type +{ + using Surface_t = typename SurfaceTypeTraits::type; + return action_(surfaces_.make_surface(id)); +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6f15b0b12e..974b01157b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -187,6 +187,7 @@ celeritas_add_test(orange/surfaces/detail/QuadraticSolver.test.cc) celeritas_add_test(orange/surfaces/CylCentered.test.cc) celeritas_add_test(orange/surfaces/GeneralQuadric.test.cc) celeritas_add_test(orange/surfaces/PlaneAligned.test.cc) +celeritas_add_test(orange/surfaces/SurfaceAction.test.cc) #-----------------------------------------------------------------------------# # I/O (ROOT) diff --git a/test/orange/surfaces/CylCentered.test.cc b/test/orange/surfaces/CylCentered.test.cc index 542d40cda8..c5bae73df0 100644 --- a/test/orange/surfaces/CylCentered.test.cc +++ b/test/orange/surfaces/CylCentered.test.cc @@ -32,7 +32,7 @@ using VecReal = std::vector; //---------------------------------------------------------------------------// TEST(TestCCylX, construction) { - EXPECT_EQ(1, CCylX::size()); + EXPECT_EQ(1, CCylX::Storage::extent); EXPECT_EQ(2, CCylX::Intersections{}.size()); CCylX c(4.0); diff --git a/test/orange/surfaces/SurfaceAction.test.cc b/test/orange/surfaces/SurfaceAction.test.cc new file mode 100644 index 0000000000..c9a7af5d5c --- /dev/null +++ b/test/orange/surfaces/SurfaceAction.test.cc @@ -0,0 +1,101 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceAction.test.cc +//---------------------------------------------------------------------------// +#include "orange/surfaces/SurfaceAction.hh" + +#include +#include +#include + +#include "base/CollectionMirror.hh" +#include "base/Range.hh" +#include "orange/Data.hh" +#include "orange/construct/SurfaceInserter.hh" +#include "orange/surfaces/Surfaces.hh" +#include "orange/surfaces/SurfaceIO.hh" +#include "celeritas_test.hh" +// #include "SurfaceAction.test.hh" + +using namespace celeritas; +// using namespace celeritas_test; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class SurfaceActionTest : public celeritas::Test +{ + protected: + using SurfaceDataMirror = CollectionMirror; + + void SetUp() override + { + SurfaceData surface_data; + SurfaceInserter insert(&surface_data); + insert(PlaneX(1)); + insert(PlaneY(2)); + insert(PlaneZ(3)); + insert(CCylX(5)); + insert(CCylY(6)); + insert(CCylZ(7)); + insert(GeneralQuadric({0, 1, 2}, {3, 4, 5}, {6, 7, 8}, 9)); + + surf_params_ = SurfaceDataMirror{std::move(surface_data)}; + } + + SurfaceDataMirror surf_params_; +}; + +//---------------------------------------------------------------------------// +// HELPERS +//---------------------------------------------------------------------------// + +struct ToString +{ + template + std::string operator()(S&& surf) const + { + std::ostringstream os; + os << surf; + return os.str(); + } + + // Make sure this test class is being move-constructed + ToString() = default; + ToString(const ToString&) = delete; + ToString(ToString&&) = default; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(SurfaceActionTest, string) +{ + // Create functor + Surfaces surfaces(surf_params_.host()); + auto surf_to_string = make_surface_action(surfaces, ToString{}); + + // Loop over all surfaces and apply + std::vector strings; + for (auto id : range(SurfaceId{surfaces.num_surfaces()})) + { + strings.push_back(surf_to_string(id)); + } + + // clang-format off + const std::string expected_strings[] = { + "Plane: x=1", + "Plane: y=2", + "Plane: z=3", + "Cyl x: r=5", + "Cyl y: r=6", + "Cyl z: r=7", + "GQuadric: {0,1,2} {3,4,5} {6,7,8} 9"}; + // clang-format on + EXPECT_VEC_EQ(expected_strings, strings); +} From 76d0fb97ddb646a6846c0fbd10733809fa69dc8f Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sun, 24 Oct 2021 08:33:42 -0400 Subject: [PATCH 07/22] Add min_element with indices --- src/base/Algorithms.hh | 16 ++++++++++++++++ test/base/Algorithms.test.cc | 16 ++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/src/base/Algorithms.hh b/src/base/Algorithms.hh index ae0374887d..5fbd7970a5 100644 --- a/src/base/Algorithms.hh +++ b/src/base/Algorithms.hh @@ -80,6 +80,22 @@ CELER_CONSTEXPR_FUNCTION const T& max(const T& a, const T& b) noexcept return (b > a) ? b : a; } +//---------------------------------------------------------------------------// +/*! + * Return the index of the lowest value. + */ +template +inline CELER_FUNCTION size_type min_element(const T* ptr, size_type length) +{ + size_type result = 0; + for (size_type i = 1; i < length; ++i) + { + if (ptr[i] < ptr[result]) + result = i; + } + return result; +} + //---------------------------------------------------------------------------// // Replace/extend //---------------------------------------------------------------------------// diff --git a/test/base/Algorithms.test.cc b/test/base/Algorithms.test.cc index e9ecfe1a1f..746d4a78e5 100644 --- a/test/base/Algorithms.test.cc +++ b/test/base/Algorithms.test.cc @@ -61,3 +61,19 @@ TEST(AlgorithmsTest, lower_bound) } } } + +TEST(AlgorithmsTest, min_element) +{ + std::vector v; + // Empty vector will return 0, which is off-the-end + EXPECT_EQ(0, celeritas::min_element(v.data(), v.size())); + + v = {100}; + EXPECT_EQ(0, celeritas::min_element(v.data(), v.size())); + + v = {10, 2, 100, 3, -1}; + EXPECT_EQ(4, celeritas::min_element(v.data(), v.size())); + + v[2] = -100; + EXPECT_EQ(2, celeritas::min_element(v.data(), v.size())); +} From 23401c05df5f509dd4bab5766873a6d63dcaac41 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sun, 24 Oct 2021 08:39:18 -0400 Subject: [PATCH 08/22] Third-guess myself and change back to iterators --- src/base/Algorithms.hh | 18 +++++++++++------- test/base/Algorithms.test.cc | 12 ++++++++---- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/base/Algorithms.hh b/src/base/Algorithms.hh index 5fbd7970a5..a60fac304e 100644 --- a/src/base/Algorithms.hh +++ b/src/base/Algorithms.hh @@ -82,16 +82,20 @@ CELER_CONSTEXPR_FUNCTION const T& max(const T& a, const T& b) noexcept //---------------------------------------------------------------------------// /*! - * Return the index of the lowest value. + * Return the lower of two values. */ -template -inline CELER_FUNCTION size_type min_element(const T* ptr, size_type length) +template +inline CELER_FUNCTION ForwardIt min_element(ForwardIt iter, ForwardIt last) { - size_type result = 0; - for (size_type i = 1; i < length; ++i) + // Avoid incrementing past the end + if (iter == last) + return last; + + ForwardIt result = iter++; + for (; iter != last; ++iter) { - if (ptr[i] < ptr[result]) - result = i; + if (*iter < *result) + result = iter; } return result; } diff --git a/test/base/Algorithms.test.cc b/test/base/Algorithms.test.cc index 746d4a78e5..14e0058afc 100644 --- a/test/base/Algorithms.test.cc +++ b/test/base/Algorithms.test.cc @@ -65,15 +65,19 @@ TEST(AlgorithmsTest, lower_bound) TEST(AlgorithmsTest, min_element) { std::vector v; + auto min_element_idx = [&v]() { + return celeritas::min_element(v.begin(), v.end()) - v.begin(); + }; + // Empty vector will return 0, which is off-the-end - EXPECT_EQ(0, celeritas::min_element(v.data(), v.size())); + EXPECT_EQ(0, min_element_idx()); v = {100}; - EXPECT_EQ(0, celeritas::min_element(v.data(), v.size())); + EXPECT_EQ(0, min_element_idx()); v = {10, 2, 100, 3, -1}; - EXPECT_EQ(4, celeritas::min_element(v.data(), v.size())); + EXPECT_EQ(4, min_element_idx()); v[2] = -100; - EXPECT_EQ(2, celeritas::min_element(v.data(), v.size())); + EXPECT_EQ(2, min_element_idx()); } From 49c10cb4dbb89be8d23ca1f82f7a176cc6d26530 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sun, 24 Oct 2021 21:22:41 -0400 Subject: [PATCH 09/22] Add CUDA surface action test --- src/base/Span.hh | 4 +- src/orange/Types.hh | 9 ++ test/CMakeLists.txt | 3 +- test/orange/surfaces/SurfaceAction.test.cc | 144 +++++++++++++++++- test/orange/surfaces/SurfaceAction.test.cu | 48 ++++++ test/orange/surfaces/SurfaceAction.test.hh | 167 +++++++++++++++++++++ test/physics/base/Particle.test.cc | 2 - test/physics/base/Particle.test.hh | 7 + 8 files changed, 377 insertions(+), 7 deletions(-) create mode 100644 test/orange/surfaces/SurfaceAction.test.cu create mode 100644 test/orange/surfaces/SurfaceAction.test.hh diff --git a/src/base/Span.hh b/src/base/Span.hh index 3eda5a3b89..993e5b77dd 100644 --- a/src/base/Span.hh +++ b/src/base/Span.hh @@ -51,7 +51,7 @@ class Span //!@} //! Size (may be dynamic) - static constexpr CELER_FUNCTION std::size_t extent = Extent; + static constexpr std::size_t extent = Extent; public: //// CONSTRUCTION //// @@ -146,7 +146,7 @@ class Span }; template -constexpr CELER_FUNCTION std::size_t Span::extent; +constexpr std::size_t Span::extent; //---------------------------------------------------------------------------// // HELPER FUNCTIONS diff --git a/src/orange/Types.hh b/src/orange/Types.hh index a13cdf4ebc..9186937a3b 100644 --- a/src/orange/Types.hh +++ b/src/orange/Types.hh @@ -210,6 +210,15 @@ CELER_CONSTEXPR_FUNCTION Sense to_sense(SignedSense s) return Sense(static_cast(s) >= 0); } +//---------------------------------------------------------------------------// +/*! + * Convert a signed sense to a surface state. + */ +CELER_CONSTEXPR_FUNCTION SurfaceState to_surface_state(SignedSense s) +{ + return s == SignedSense::on ? SurfaceState::on : SurfaceState::off; +} + //---------------------------------------------------------------------------// /*! * Sentinel value indicating "no intersection". diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 974b01157b..3648878bf0 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -187,7 +187,8 @@ celeritas_add_test(orange/surfaces/detail/QuadraticSolver.test.cc) celeritas_add_test(orange/surfaces/CylCentered.test.cc) celeritas_add_test(orange/surfaces/GeneralQuadric.test.cc) celeritas_add_test(orange/surfaces/PlaneAligned.test.cc) -celeritas_add_test(orange/surfaces/SurfaceAction.test.cc) + +celeritas_cudaoptional_test(orange/surfaces/SurfaceAction) #-----------------------------------------------------------------------------# # I/O (ROOT) diff --git a/test/orange/surfaces/SurfaceAction.test.cc b/test/orange/surfaces/SurfaceAction.test.cc index c9a7af5d5c..d0a14bf44d 100644 --- a/test/orange/surfaces/SurfaceAction.test.cc +++ b/test/orange/surfaces/SurfaceAction.test.cc @@ -7,6 +7,9 @@ //---------------------------------------------------------------------------// #include "orange/surfaces/SurfaceAction.hh" +#include +#include +#include #include #include #include @@ -17,11 +20,47 @@ #include "orange/construct/SurfaceInserter.hh" #include "orange/surfaces/Surfaces.hh" #include "orange/surfaces/SurfaceIO.hh" +#include "random/distributions/IsotropicDistribution.hh" +#include "random/distributions/UniformBoxDistribution.hh" #include "celeritas_test.hh" -// #include "SurfaceAction.test.hh" +#include "SurfaceAction.test.hh" using namespace celeritas; -// using namespace celeritas_test; +using namespace celeritas_test; + +namespace celeritas +{ +std::ostream& operator<<(std::ostream& os, Sense s) +{ + return os << to_char(s); +} +} // namespace celeritas + +namespace +{ +#if 0 +std::vector string_to_senses(std::string s) +{ + std::vector result(s.size()); + std::transform(s.begin(), s.end(), result.begin(), [](char c) { + CELER_EXPECT(c == '+' || c == '-'); + return c == '+' ? Sense::outside : Sense::inside; + }); + return result; +} +#endif + +std::string senses_to_string(Span s) +{ + std::string result(s.size(), ' '); + std::transform(s.begin(), s.end(), result.begin(), [](Sense s) { + return to_char(s); + }); + return result; +} + +constexpr real_type noint = celeritas::no_intersection(); +} // namespace //---------------------------------------------------------------------------// // TEST HARNESS @@ -47,7 +86,26 @@ class SurfaceActionTest : public celeritas::Test surf_params_ = SurfaceDataMirror{std::move(surface_data)}; } + void fill_uniform_box(celeritas::Span pos) + { + UniformBoxDistribution<> sample_box{{-3, -2, -1}, {6, 8, 10}}; + for (Real3& d : pos) + { + d = sample_box(rng_); + } + } + + void fill_isotropic(celeritas::Span dir) + { + IsotropicDistribution<> sample_isotropic; + for (Real3& d : dir) + { + d = sample_isotropic(rng_); + } + } + SurfaceDataMirror surf_params_; + std::mt19937 rng_; }; //---------------------------------------------------------------------------// @@ -99,3 +157,85 @@ TEST_F(SurfaceActionTest, string) // clang-format on EXPECT_VEC_EQ(expected_strings, strings); } + +TEST_F(SurfaceActionTest, host_distances) +{ + // Create states and sample uniform box, isotropic direction + OrangeMiniStateData states; + resize(&states, surf_params_.host(), 1024); + OrangeMiniStateData state_ref; + state_ref = states; + this->fill_uniform_box(state_ref.pos[AllItems{}]); + this->fill_isotropic(state_ref.dir[AllItems{}]); + + CalcSenseDistanceLauncher<> calc_thread{surf_params_.host(), state_ref}; + for (auto tid : range(ThreadId{states.size()})) + { + calc_thread(tid); + } + + auto test_threads = celeritas::range(ThreadId{10}); + // PRINT_EXPECTED(senses_to_string(state_ref.sense[test_threads])); + // PRINT_EXPECTED(state_ref.distance[test_threads]); + + const char expected_senses[] + = {'-', '-', '+', '+', '-', '-', '+', '-', '+', '-'}; + const double expected_distance[] = {noint, + noint, + noint, + noint, + 8.623486582635, + 8.115429697208, + noint, + noint, + 22.04764572126, + noint}; + EXPECT_VEC_EQ(expected_senses, + senses_to_string(state_ref.sense[test_threads])); + EXPECT_VEC_SOFT_EQ(expected_distance, state_ref.distance[test_threads]); +} + +TEST_F(SurfaceActionTest, TEST_IF_CELERITAS_CUDA(device_distances)) +{ + OrangeMiniStateData device_states; + { + // Initialize on host + OrangeMiniStateData host_states; + resize(&host_states, surf_params_.host(), 1024); + this->fill_uniform_box(host_states.pos[AllItems{}]); + this->fill_isotropic(host_states.dir[AllItems{}]); + + // Copy starting position/direction to device + device_states = host_states; + } + + // Launch kernel + SATestInput input; + input.params = surf_params_.device(); + input.states = device_states; + sa_test(input); + + { + // Copy result back to host + OrangeMiniStateData host_states; + host_states = device_states; + auto test_threads = celeritas::range(ThreadId{10}); + + const char expected_senses[] + = {'-', '-', '+', '+', '-', '-', '+', '-', '+', '-'}; + const double expected_distance[] = {noint, + noint, + noint, + noint, + 8.623486582635, + 8.115429697208, + noint, + noint, + 22.04764572126, + noint}; + EXPECT_VEC_EQ(expected_senses, + senses_to_string(host_states.sense[test_threads])); + EXPECT_VEC_SOFT_EQ(expected_distance, + host_states.distance[test_threads]); + } +} diff --git a/test/orange/surfaces/SurfaceAction.test.cu b/test/orange/surfaces/SurfaceAction.test.cu new file mode 100644 index 0000000000..3fc69007ae --- /dev/null +++ b/test/orange/surfaces/SurfaceAction.test.cu @@ -0,0 +1,48 @@ +//---------------------------------*-CUDA-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceAction.test.cu +//---------------------------------------------------------------------------// +#include "SurfaceAction.test.hh" + +#include "base/KernelParamCalculator.cuda.hh" + +namespace celeritas_test +{ +namespace +{ +//---------------------------------------------------------------------------// +// KERNELS +//---------------------------------------------------------------------------// + +__global__ void sa_test_kernel(SATestInput input) +{ + auto tid = celeritas::KernelParamCalculator::thread_id(); + if (tid.get() >= input.states.size()) + return; + + // Calculate distances in parallel + CalcSenseDistanceLauncher<> calc_thread{input.params, input.states}; + calc_thread(tid); +} +} // namespace + +//---------------------------------------------------------------------------// +// TESTING INTERFACE +//---------------------------------------------------------------------------// +//! Run on device and return results +void sa_test(SATestInput input) +{ + static const celeritas::KernelParamCalculator calc_launch_params( + sa_test_kernel, "sa_test"); + auto params = calc_launch_params(input.states.size()); + sa_test_kernel<<>>(input); + + CELER_CUDA_CHECK_ERROR(); + CELER_CUDA_CALL(cudaDeviceSynchronize()); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas_test diff --git a/test/orange/surfaces/SurfaceAction.test.hh b/test/orange/surfaces/SurfaceAction.test.hh new file mode 100644 index 0000000000..777c1f6e26 --- /dev/null +++ b/test/orange/surfaces/SurfaceAction.test.hh @@ -0,0 +1,167 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file SurfaceAction.test.hh +//---------------------------------------------------------------------------// +#include +#include "orange/Data.hh" +#include "orange/surfaces/Surfaces.hh" +#include "orange/surfaces/SurfaceAction.hh" + +namespace celeritas_test +{ +using celeritas::MemSpace; +using celeritas::Ownership; +using celeritas::Real3; +using celeritas::real_type; +using celeritas::ThreadId; + +using celeritas::Sense; +using celeritas::SurfaceState; + +//---------------------------------------------------------------------------// +/*! + * Small subset of state for testing. + */ +template +struct OrangeMiniStateData +{ + template + using StateItems = celeritas::StateCollection; + + StateItems pos; + StateItems dir; + StateItems sense; + StateItems distance; + + //! True if sizes are consistent and nonzero + explicit CELER_FUNCTION operator bool() const + { + // clang-format off + return pos.size() > 0 + && dir.size() == pos.size() + && sense.size() == pos.size() + && distance.size() == pos.size(); + // clang-format on + } + + //! Assign from another set of data + template + OrangeMiniStateData& operator=(OrangeMiniStateData& other) + { + CELER_EXPECT(other); + pos = other.pos; + dir = other.dir; + sense = other.sense; + distance = other.distance; + CELER_ENSURE(*this); + return *this; + } + + //! State size + CELER_FUNCTION ThreadId::size_type size() const { return pos.size(); } +}; + +//---------------------------------------------------------------------------// +/*! + * Resize states in host code. + */ +template +inline void +resize(OrangeMiniStateData* data, + const celeritas::SurfaceData&, + celeritas::size_type size) +{ + CELER_EXPECT(data); + CELER_EXPECT(size > 0); + make_builder(&data->pos).resize(size); + make_builder(&data->dir).resize(size); + make_builder(&data->sense).resize(size); + make_builder(&data->distance).resize(size); + + CELER_ENSURE(*data); +} + +//---------------------------------------------------------------------------// +// HELPER FUNCTIONS +//---------------------------------------------------------------------------// +//! Calculate the sense of a surface at a given position. +struct CalcSenseDistance +{ + const Real3& pos; + const Real3& dir; + Sense* sense; + real_type* distance; + + template + CELER_FUNCTION void operator()(S&& surf) + { + // Calculate sense + auto signed_sense = surf.calc_sense(this->pos); + *this->sense = to_sense(signed_sense); + + // Calculate nearest distance + auto intersect = surf.calc_intersections( + this->pos, this->dir, to_surface_state(signed_sense)); + for (real_type distance : intersect) + { + CELER_ASSERT(distance > 0); + } + *this->distance + = *celeritas::min_element(intersect.begin(), intersect.end()); + } +}; + +//! Calculate distance to a single surface +template +struct CalcSenseDistanceLauncher +{ + celeritas::SurfaceData params; + OrangeMiniStateData states; + + CELER_FUNCTION void operator()(ThreadId tid) const + { + celeritas::Surfaces surfaces(this->params); + + auto calc_sense_dist = celeritas::make_surface_action( + surfaces, + CalcSenseDistance{this->states.pos[tid], + this->states.dir[tid], + &this->states.sense[tid], + &this->states.distance[tid]}); + + celeritas::SurfaceId sid{tid.get() % surfaces.num_surfaces()}; + calc_sense_dist(sid); + } +}; + +//---------------------------------------------------------------------------// +// TESTING INTERFACE +//---------------------------------------------------------------------------// +//! Input data +struct SATestInput +{ + using ParamsRef + = celeritas::SurfaceData; + using StateRef + = OrangeMiniStateData; + + ParamsRef params; + StateRef states; +}; + +//---------------------------------------------------------------------------// +//! Run on device and return results +void sa_test(SATestInput); + +#if !CELERITAS_USE_CUDA +inline void sa_test(SATestInput) +{ + CELER_NOT_CONFIGURED("CUDA"); +} +#endif + +//---------------------------------------------------------------------------// +} // namespace celeritas_test diff --git a/test/physics/base/Particle.test.cc b/test/physics/base/Particle.test.cc index b1ac757cb9..a5ffc5f91d 100644 --- a/test/physics/base/Particle.test.cc +++ b/test/physics/base/Particle.test.cc @@ -226,9 +226,7 @@ TEST_F(ParticleDeviceTest, TEST_IF_CELERITAS_CUDA(calc_props)) // Run GPU test PTVTestOutput result; -#if CELERITAS_USE_CUDA result = ptv_test(input); -#endif // Check results const double expected_props[] = {0.5, diff --git a/test/physics/base/Particle.test.hh b/test/physics/base/Particle.test.hh index 7ba521fe97..2e08c72333 100644 --- a/test/physics/base/Particle.test.hh +++ b/test/physics/base/Particle.test.hh @@ -42,5 +42,12 @@ struct PTVTestOutput //! Run on device and return results PTVTestOutput ptv_test(PTVTestInput); +#if !CELERITAS_USE_CUDA +inline PTVTestOutput ptv_test(PTVTestInput) +{ + CELER_NOT_CONFIGURED("CUDA"); +} +#endif + //---------------------------------------------------------------------------// } // namespace celeritas_test From 1ad2c644c9bd5705abd9f22478f806499232a04e Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 25 Oct 2021 20:06:41 -0400 Subject: [PATCH 10/22] Define CELER_UNREACHABLE on device Before: ``` /home/s3j/.local/src/celeritas/src/orange/surfaces/detail/SurfaceAction.hh(112): warning: calling a __host__ function("__builtin_unreachable") from a __host__ __device__ function("celeritas::detail::SurfaceAction< ::celeritas_test::CalcSenseDistance> ::operator ()") is not allowed ``` Using `__trap` didn't reduce the number of instructions, but simply omitting the warning works. --- src/base/Macros.hh | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/base/Macros.hh b/src/base/Macros.hh index 627c446e48..5cb7755453 100644 --- a/src/base/Macros.hh +++ b/src/base/Macros.hh @@ -92,12 +92,19 @@ * * See https://clang.llvm.org/docs/LanguageExtensions.html#builtin-unreachable * or https://msdn.microsoft.com/en-us/library/1b3fsfxw.aspx= + * or + * https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#__builtin_unreachable + * + * (The 'unreachable' and 'assume' compiler optimizations for CUDA are only + * available in API version 11.1 or higher, which is encoded as major*1000 + + * minor*10). * * \note This macro should not generally be used; instead, the macro \c * CELER_ASSERT_UNREACHABLE() defined in base/Assert.hh should be used instead * (to provide a more detailed error message in case the point *is* reached). */ -#if defined(__clang__) || defined(__GNUC__) +#if (!defined(__CUDA_ARCH__) && (defined(__clang__) || defined(__GNUC__))) \ + || (defined(__CUDA_ARCH__) && CUDART_VERSION >= 11010) # define CELER_UNREACHABLE __builtin_unreachable() #elif defined(_MSC_VER) # define CELER_UNREACHABLE __assume(false) From dd710dcd7f869d31ec51b9c112f3c60693a62dae Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 26 Oct 2021 11:15:04 -0400 Subject: [PATCH 11/22] fixup! Add surface inserter, data interface, and surface action --- src/orange/surfaces/Surfaces.hh | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/orange/surfaces/Surfaces.hh b/src/orange/surfaces/Surfaces.hh index 2b5a24a6a0..272d8517ec 100644 --- a/src/orange/surfaces/Surfaces.hh +++ b/src/orange/surfaces/Surfaces.hh @@ -13,7 +13,7 @@ namespace celeritas { //---------------------------------------------------------------------------// /*! - * Operations on one or more surfaces. + * Access stored surface data. * * These are all surfaces in the entire geometry. */ @@ -22,12 +22,13 @@ class Surfaces public: //@{ //! Type aliases - using Data = SurfaceData; + using SurfaceDataRef + = SurfaceData; //@} public: - // Construct with defaults - explicit inline CELER_FUNCTION Surfaces(const Data&); + // Construct with reference to persistent data + explicit inline CELER_FUNCTION Surfaces(const SurfaceDataRef&); // Number of surfaces inline CELER_FUNCTION SurfaceId::size_type num_surfaces() const; @@ -40,16 +41,16 @@ class Surfaces inline CELER_FUNCTION T make_surface(SurfaceId) const; private: - const Data& data_; + const SurfaceDataRef& data_; }; //---------------------------------------------------------------------------// // INLINE FUNCTION DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with defaults. + * Construct with reference to persistent data. */ -CELER_FUNCTION Surfaces::Surfaces(const Data& data) : data_(data) +CELER_FUNCTION Surfaces::Surfaces(const SurfaceDataRef& data) : data_(data) { CELER_EXPECT(data_); } From f82b255329120d3d0270d4798a5de248c29c73c5 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 26 Oct 2021 11:31:36 -0400 Subject: [PATCH 12/22] fixup! Add aligned plane --- test/orange/surfaces/PlaneAligned.test.cc | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/test/orange/surfaces/PlaneAligned.test.cc b/test/orange/surfaces/PlaneAligned.test.cc index c9b6c5f082..13b6ec8cb6 100644 --- a/test/orange/surfaces/PlaneAligned.test.cc +++ b/test/orange/surfaces/PlaneAligned.test.cc @@ -69,8 +69,6 @@ TEST_F(PlaneAlignedTest, x_plane) calc_intersection(p, {0.9999, 0.0, 0.0}, px, SurfaceState::on)); EXPECT_EQ(no_intersection(), calc_intersection(p, {1.0001, 0.0, 0.0}, mx, SurfaceState::on)); - EXPECT_SOFT_EQ(0.0, calc_intersection(p, {1.0, 0.0, 0.0}, px)); - EXPECT_SOFT_EQ(0.0, calc_intersection(p, {1.0, 0.0, 0.0}, mx)); EXPECT_SOFT_EQ(0.01, calc_intersection(p, {0.99, 0.0, 0.0}, px)); EXPECT_SOFT_EQ(0.01, calc_intersection(p, {1.01, 0.0, 0.0}, mx)); EXPECT_SOFT_EQ(no_intersection(), @@ -99,10 +97,10 @@ TEST_F(PlaneAlignedTest, y_plane) Real3 px{1.0, 0.0, 0.0}; EXPECT_EQ(no_intersection(), - calc_intersection(p, {0, 1 - 1e-8, 0.0}, py, SurfaceState::on)); + calc_intersection(p, {0, -1 + 1e-8, 0.0}, py, SurfaceState::on)); EXPECT_SOFT_EQ(0.01, calc_intersection(p, {0.0, -1.01, 0.0}, py)); EXPECT_SOFT_EQ(no_intersection(), - calc_intersection(p, {0.0, 0.99, 0.0}, my)); + calc_intersection(p, {0.0, -1.1, 0.0}, my)); EXPECT_EQ(no_intersection(), calc_intersection(p, {-1.01, 1.0, 0.0}, px)); } @@ -125,8 +123,6 @@ TEST_F(PlaneAlignedTest, plane_z) EXPECT_EQ(no_intersection(), calc_intersection(p, {0.0, 0.0, -1e-8}, pz, SurfaceState::on)); - EXPECT_SOFT_EQ(0.0, calc_intersection(p, {0.0, 0.0, 0.0}, pz)); - EXPECT_SOFT_EQ(0.0, calc_intersection(p, {0.0, 0.0, 0.0}, mz)); EXPECT_SOFT_EQ(0.01, calc_intersection(p, {0.0, 0.0, -0.01}, pz)); EXPECT_SOFT_EQ(0.01, calc_intersection(p, {0.0, 0.0, 0.01}, mz)); EXPECT_EQ(no_intersection(), calc_intersection(p, {-1.01, 0.0, 0.0}, px)); From d562fdccfa75038f8b78eb262d0680f2340d58ec Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 26 Oct 2021 18:54:00 -0400 Subject: [PATCH 13/22] Document `size_` members --- src/base/Range.hh | 3 ++- src/orange/Types.hh | 12 ++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/base/Range.hh b/src/base/Range.hh index 9a91d9b849..008e577fd1 100644 --- a/src/base/Range.hh +++ b/src/base/Range.hh @@ -47,7 +47,8 @@ using RangeIter = detail::range_iter; * * Here, T can be any of: * - an integer, - * - an enum that has a "size_" member, or + * - an enum that has contiguous zero-indexed values and a "size_" enumeration + * value indicating how many, or * - an OpaqueId. * * It is OK to dereference the end iterator! The result should just be the diff --git a/src/orange/Types.hh b/src/orange/Types.hh index 9186937a3b..e3e812280a 100644 --- a/src/orange/Types.hh +++ b/src/orange/Types.hh @@ -54,10 +54,10 @@ enum class Sense : bool */ enum class Axis { - x, - y, - z, - size_ + x, //!< X axis/I index coordinate + y, //!< Y axis/J index coordinate + z, //!< Z axis/K index coordinate + size_ //!< Sentinel value for looping over axes }; //---------------------------------------------------------------------------// @@ -87,8 +87,8 @@ enum class SurfaceType : unsigned char kz, //!< Cone parallel to Z axis sq, //!< Simple quadric #endif - gq, //!< General quadric - size_ + gq, //!< General quadric + size_ //!< Sentinel value for number of surface types }; //---------------------------------------------------------------------------// From 129b3aac55ad97c0493d5f3cd5d32780b3eed131 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 26 Oct 2021 18:59:45 -0400 Subject: [PATCH 14/22] Defer unused extra data/type definitions --- src/orange/Data.hh | 147 -------------------------------------------- src/orange/Types.hh | 3 - 2 files changed, 150 deletions(-) diff --git a/src/orange/Data.hh b/src/orange/Data.hh index ccbd814438..036a02e78f 100644 --- a/src/orange/Data.hh +++ b/src/orange/Data.hh @@ -8,7 +8,6 @@ #pragma once #include "base/Collection.hh" -#include "base/CollectionBuilder.hh" #include "base/OpaqueId.hh" #include "Types.hh" @@ -60,151 +59,5 @@ struct SurfaceData } }; -//---------------------------------------------------------------------------// -/*! - * Data for a single volume definition. - * - * A volume is a CSG tree of surfaces. Each surface defines an "inside" space - * and "outside" space that correspond to "negative" and "positive" values of - * the quadric expression's evaluation. Left of a plane is negative, for - * example, and evaluates to "false" or "inside" or "negative". The CSG tree is - * flattened into a - */ -struct VolumeDef -{ - ItemRange faces; - ItemRange logic; - - logic_int num_intersections{0}; - logic_int flags{0}; -}; - -//---------------------------------------------------------------------------// -/*! - * Data for volume definitions. - */ -template -struct VolumeData -{ - //// TYPES //// - - template - using Items = Collection; - - //// DATA //// - - Items defs; - - // Storage - Collection faces; - Collection logic; - - //// METHODS //// - - //! True if sizes are valid - explicit CELER_FUNCTION operator bool() const { return defs.size() > 0; } - - //! Assign from another set of data - template - VolumeData& operator=(const VolumeData& other) - { - CELER_EXPECT(other); - - defs = other.defs; - faces = other.faces; - logic = other.logic; - - return *this; - } -}; - -//---------------------------------------------------------------------------// -/*! - * Data for universe definitions. - */ -template -struct UniverseData -{ - //// TYPES //// - - template - using Items = Collection; - - //// DATA //// - - //// METHODS //// - - //! True if sizes are valid - explicit CELER_FUNCTION operator bool() const { return false; } - - //! Assign from another set of data - template - UniverseData& operator=(const UniverseData& other) - { - CELER_EXPECT(other); - return *this; - } -}; - -//---------------------------------------------------------------------------// -/*! - * Scalar values particular to an ORANGE geometry instance. - */ -struct OrangeParamsScalars -{ - size_type max_level{}; - size_type max_faces{}; - size_type max_intersections{}; - - // TODO: fuzziness/length scale -}; - -//---------------------------------------------------------------------------// -/*! - * Data to persistent data used by ORANGE implementation. - */ -template -struct OrangeParamsData -{ - //// DATA //// - - SurfaceData surfaces; - VolumeData volumes; - UniverseData universes; - - OrangeParamsScalars scalars; - - //// METHODS //// - - //! True if assigned - explicit CELER_FUNCTION operator bool() const - { - return surfaces && volumes && universes; - } - - //! Assign from another set of data - template - OrangeParamsData& operator=(const OrangeParamsData& other) - { - CELER_EXPECT(other); - surfaces = other.surfaces; - volumes = other.volumes; - universes = other.universes; - scalars = other.scalars; - return *this; - } -}; - -//---------------------------------------------------------------------------// -// STATE -//---------------------------------------------------------------------------// -/*! - * Data to persistent data used by ORANGE implementation. - */ -template -struct OrangeStateData -{ -}; - //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/orange/Types.hh b/src/orange/Types.hh index e3e812280a..92fd0129d3 100644 --- a/src/orange/Types.hh +++ b/src/orange/Types.hh @@ -25,9 +25,6 @@ using logic_int = unsigned short int; //! Identifier for a surface in a universe using SurfaceId = OpaqueId; -//! Identifier for a relocatable set of volumes -using UniverseId = OpaqueId; - //---------------------------------------------------------------------------// // ENUMERATIONS //---------------------------------------------------------------------------// From 159f7b6284c05584f3dc14a3849d839b96608285 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 26 Oct 2021 19:14:41 -0400 Subject: [PATCH 15/22] fixup! Defer unused extra data/type definitions --- test/orange/surfaces/SurfaceAction.test.hh | 1 + 1 file changed, 1 insertion(+) diff --git a/test/orange/surfaces/SurfaceAction.test.hh b/test/orange/surfaces/SurfaceAction.test.hh index 777c1f6e26..a6ac1db96d 100644 --- a/test/orange/surfaces/SurfaceAction.test.hh +++ b/test/orange/surfaces/SurfaceAction.test.hh @@ -6,6 +6,7 @@ //! \file SurfaceAction.test.hh //---------------------------------------------------------------------------// #include +#include "base/CollectionBuilder.hh" #include "orange/Data.hh" #include "orange/surfaces/Surfaces.hh" #include "orange/surfaces/SurfaceAction.hh" From 5d38b71ff640d911464912b5d1c6f2f187f87958 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 26 Oct 2021 19:57:18 -0400 Subject: [PATCH 16/22] Use CUDA 10 on emmet and add surface assertion --- scripts/build/emmet.sh | 2 +- src/orange/surfaces/Surfaces.hh | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/scripts/build/emmet.sh b/scripts/build/emmet.sh index 6ef997f77a..6fb26b1778 100755 --- a/scripts/build/emmet.sh +++ b/scripts/build/emmet.sh @@ -13,7 +13,7 @@ mkdir ${BUILD_DIR} 2>/dev/null \ cd ${BUILD_DIR} module purge -module load cuda +module load cuda/10 CELERITAS_ENV=$SPACK_ROOT/var/spack/environments/celeritas/.spack-env/view export PATH=$CELERITAS_ENV/bin:${PATH} export CMAKE_PREFIX_PATH=$CELERITAS_ENV:${CMAKE_PREFIX_PATH} diff --git a/src/orange/surfaces/Surfaces.hh b/src/orange/surfaces/Surfaces.hh index 272d8517ec..115a1ec5b3 100644 --- a/src/orange/surfaces/Surfaces.hh +++ b/src/orange/surfaces/Surfaces.hh @@ -84,6 +84,8 @@ CELER_FUNCTION T Surfaces::make_surface(SurfaceId sid) const static_assert(T::Storage::extent > 0, "Template parameter must be a surface class"); CELER_EXPECT(sid < this->num_surfaces()); + CELER_EXPECT(this->surface_type(sid) == T::surface_type()); + OpaqueId start = data_.offsets[sid]; OpaqueId end{start.get() + static_cast(T::Storage::extent)}; From 25d54d6fc6e63c1b81e823943a6442e163bff894 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Wed, 27 Oct 2021 22:49:43 -0400 Subject: [PATCH 17/22] Add sqrt(2) and sqrt(3) constexpr constants --- src/base/Constants.hh | 8 +++++--- test/base/Constants.test.cc | 9 +++++++++ 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/src/base/Constants.hh b/src/base/Constants.hh index 7c0aaa5f8f..d69c2c88d9 100644 --- a/src/base/Constants.hh +++ b/src/base/Constants.hh @@ -47,9 +47,11 @@ namespace constants */ //!@{ -//! Mathemetical constant -constexpr real_type pi = 3.14159265358979323846; // truncated -constexpr real_type euler = 2.71828182845904523536; +//! Mathemetical constants (truncated) +constexpr real_type pi = 3.14159265358979323846; +constexpr real_type euler = 2.71828182845904523536; +constexpr real_type sqrt_two = 1.41421356237309504880; +constexpr real_type sqrt_three = 1.73205080756887729353; //!@} //!@{ diff --git a/test/base/Constants.test.cc b/test/base/Constants.test.cc index b9e56c6dc0..e905d6e57e 100644 --- a/test/base/Constants.test.cc +++ b/test/base/Constants.test.cc @@ -8,6 +8,7 @@ #include "base/Constants.hh" #include "base/Units.hh" +#include #include "celeritas_test.hh" using namespace celeritas::units; @@ -18,6 +19,14 @@ using celeritas::real_type; // TESTS //---------------------------------------------------------------------------// +TEST(UnitsTest, mathematical) +{ + EXPECT_DOUBLE_EQ(euler, std::exp(1.0)); + EXPECT_DOUBLE_EQ(pi, std::acos(-1.0)); + EXPECT_DOUBLE_EQ(sqrt_two, std::sqrt(2.0)); + EXPECT_DOUBLE_EQ(sqrt_three, std::sqrt(3.0)); +} + TEST(UnitsTest, equivalence) { EXPECT_DOUBLE_EQ(ampere * ampere * second * second * second * second From 3cd051b30552da7ca8d76e50ca99e247d72dfcd6 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 25 Oct 2021 19:18:51 -0400 Subject: [PATCH 18/22] Add extent to Array This is a nonstandard extension that lets us more easily interrogate array sizes. --- src/base/Array.hh | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/base/Array.hh b/src/base/Array.hh index a28a07f47d..de7ee08d9f 100644 --- a/src/base/Array.hh +++ b/src/base/Array.hh @@ -38,6 +38,9 @@ struct Array using const_iterator = const_pointer; //!@} + //! Static size (not in std::array) + static constexpr size_type extent = N; + //// DATA //// T data_[N]; //!< Storage @@ -82,6 +85,9 @@ struct Array //!@} }; +template +constexpr size_type Array::extent; + //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// From 47dacb51ed9fad8a5d3c4e4f325ea119fb43048b Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 28 Oct 2021 06:29:11 -0400 Subject: [PATCH 19/22] Add sphere surface --- src/orange/Types.cc | 2 + src/orange/Types.hh | 4 +- src/orange/surfaces/Sphere.hh | 147 ++++++++++++++++++++ src/orange/surfaces/SurfaceIO.cc | 8 ++ src/orange/surfaces/SurfaceIO.hh | 8 +- src/orange/surfaces/SurfaceTypeTraits.hh | 3 + src/orange/surfaces/detail/SurfaceAction.hh | 3 + test/CMakeLists.txt | 1 + test/orange/surfaces/GeneralQuadric.test.cc | 4 + test/orange/surfaces/Sphere.test.cc | 95 +++++++++++++ test/orange/surfaces/SurfaceAction.test.cc | 41 +++--- 11 files changed, 292 insertions(+), 24 deletions(-) create mode 100644 src/orange/surfaces/Sphere.hh create mode 100644 test/orange/surfaces/Sphere.test.cc diff --git a/src/orange/Types.cc b/src/orange/Types.cc index c4a7372635..ff1af8e2af 100644 --- a/src/orange/Types.cc +++ b/src/orange/Types.cc @@ -32,7 +32,9 @@ const char* to_cstring(SurfaceType value) "cy", "cz", "p", +#endif "s", +#if 0 "kx", "ky", "kz", diff --git a/src/orange/Types.hh b/src/orange/Types.hh index 92fd0129d3..7521eb7614 100644 --- a/src/orange/Types.hh +++ b/src/orange/Types.hh @@ -78,7 +78,9 @@ enum class SurfaceType : unsigned char cy, //!< Cylinder parallel to Y axis cz, //!< Cylinder parallel to Z axis p, //!< General plane - s, //!< Sphere +#endif + s, //!< Sphere +#if 0 kx, //!< Cone parallel to X axis ky, //!< Cone parallel to Y axis kz, //!< Cone parallel to Z axis diff --git a/src/orange/surfaces/Sphere.hh b/src/orange/surfaces/Sphere.hh new file mode 100644 index 0000000000..a4e9b20194 --- /dev/null +++ b/src/orange/surfaces/Sphere.hh @@ -0,0 +1,147 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file Sphere.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Algorithms.hh" +#include "base/Array.hh" +#include "base/ArrayUtils.hh" +#include "base/Span.hh" +#include "base/Types.hh" +#include "detail/QuadraticSolver.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Sphere centered at an arbitrary point. + */ +class Sphere +{ + public: + //@{ + //! Type aliases + using Intersections = Array; + using Storage = Span; + using SpanConstReal3 = Span; + //@} + + //// CLASS ATTRIBUTES //// + + //! Surface type identifier + static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type() + { + return SurfaceType::s; + } + + public: + //// CONSTRUCTORS //// + + // Construct with origin and radius + explicit inline CELER_FUNCTION + Sphere(const Real3& origin, real_type radius); + + // Construct from raw data + explicit inline CELER_FUNCTION Sphere(Storage); + + //// ACCESSORS //// + + //! Second-order terms + CELER_FUNCTION const Real3& origin() const { return origin_; } + + //! Square of the radius + CELER_FUNCTION real_type radius_sq() const { return radius_sq_; } + + //! Get a view to the data for type-deleted storage + CELER_FUNCTION Storage data() const { return {origin_.data(), 4}; } + + //// CALCULATION //// + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION SignedSense calc_sense(const Real3& pos) const; + + // Determine the sense of the position relative to this surface + inline CELER_FUNCTION Intersections calc_intersections( + const Real3& pos, const Real3& dir, SurfaceState on_surface) const; + + // Calculate outward normal at a position + inline CELER_FUNCTION Real3 calc_normal(const Real3& pos) const; + + private: + // Spatial position + Real3 origin_; + // Square of the radius + real_type radius_sq_; +}; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with all coefficients. + */ +CELER_FUNCTION Sphere::Sphere(const Real3& origin, real_type radius) + : origin_(origin), radius_sq_(ipow<2>(radius)) +{ + CELER_EXPECT(radius > 0); +} + +//---------------------------------------------------------------------------// +/*! + * Construct from raw data. + */ +CELER_FUNCTION Sphere::Sphere(Storage data) + : origin_{data[0], data[1], data[2]}, radius_sq_{data[3]} +{ +} + +//---------------------------------------------------------------------------// +/*! + * Determine the sense of the position relative to this surface. + */ +CELER_FUNCTION SignedSense Sphere::calc_sense(const Real3& pos) const +{ + Real3 tpos{pos[0] - origin_[0], pos[1] - origin_[1], pos[2] - origin_[2]}; + + return real_to_sense(dot_product(tpos, tpos) - radius_sq_); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate all possible straight-line intersections with this surface. + */ +CELER_FUNCTION auto Sphere::calc_intersections(const Real3& pos, + const Real3& dir, + SurfaceState on_surface) const + -> Intersections +{ + Real3 tpos{pos[0] - origin_[0], pos[1] - origin_[1], pos[2] - origin_[2]}; + + detail::QuadraticSolver solve_quadric(real_type(1), dot_product(tpos, dir)); + if (on_surface == SurfaceState::off) + { + return solve_quadric(dot_product(tpos, tpos) - radius_sq_); + } + else + { + return solve_quadric(); + } +} + +//---------------------------------------------------------------------------// +/*! + * Calculate outward normal at a position. + */ +CELER_FUNCTION Real3 Sphere::calc_normal(const Real3& pos) const +{ + Real3 tpos{pos[0] - origin_[0], pos[1] - origin_[1], pos[2] - origin_[2]}; + normalize_direction(&tpos); + return tpos; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/surfaces/SurfaceIO.cc b/src/orange/surfaces/SurfaceIO.cc index 542954019d..5b03e3ce02 100644 --- a/src/orange/surfaces/SurfaceIO.cc +++ b/src/orange/surfaces/SurfaceIO.cc @@ -12,6 +12,7 @@ #include "CylCentered.hh" #include "GeneralQuadric.hh" #include "PlaneAligned.hh" +#include "Sphere.hh" namespace celeritas { @@ -48,6 +49,13 @@ std::ostream& operator<<(std::ostream& os, const PlaneAligned& s) ORANGE_INSTANTIATE_SHAPE_STREAM(PlaneAligned); //---------------------------------------------------------------------------// +std::ostream& operator<<(std::ostream& os, const Sphere& s) +{ + os << "Sphere: r=" << std::sqrt(s.radius_sq()) << " at " + << make_span(s.origin()); + return os; +} +//---------------------------------------------------------------------------// #undef ORANGE_INSTANTIATE_SHAPE_STREAM } // namespace celeritas diff --git a/src/orange/surfaces/SurfaceIO.hh b/src/orange/surfaces/SurfaceIO.hh index 4be6c0cacc..d38a064be6 100644 --- a/src/orange/surfaces/SurfaceIO.hh +++ b/src/orange/surfaces/SurfaceIO.hh @@ -16,13 +16,15 @@ namespace celeritas //---------------------------------------------------------------------------// //!@{ //! Print surfaces to a stream. -template -std::ostream& operator<<(std::ostream& os, const PlaneAligned& s); - template std::ostream& operator<<(std::ostream& os, const CylCentered& s); std::ostream& operator<<(std::ostream& os, const GeneralQuadric& s); + +template +std::ostream& operator<<(std::ostream& os, const PlaneAligned& s); + +std::ostream& operator<<(std::ostream& os, const Sphere& s); //!@} //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/orange/surfaces/SurfaceTypeTraits.hh b/src/orange/surfaces/SurfaceTypeTraits.hh index e5ad7d3fb8..22965f9bdd 100644 --- a/src/orange/surfaces/SurfaceTypeTraits.hh +++ b/src/orange/surfaces/SurfaceTypeTraits.hh @@ -19,6 +19,7 @@ class PlaneAligned; template class CylCentered; class GeneralQuadric; +class Sphere; //---------------------------------------------------------------------------// /*! @@ -47,7 +48,9 @@ ORANGE_SURFACE_TRAITS(cx, CylAligned); ORANGE_SURFACE_TRAITS(cy, CylAligned); ORANGE_SURFACE_TRAITS(cz, CylAligned); ORANGE_SURFACE_TRAITS(p, Plane); +#endif ORANGE_SURFACE_TRAITS(s, Sphere); +#if 0 ORANGE_SURFACE_TRAITS(kx, ConeAligned); ORANGE_SURFACE_TRAITS(ky, ConeAligned); ORANGE_SURFACE_TRAITS(kz, ConeAligned); diff --git a/src/orange/surfaces/detail/SurfaceAction.hh b/src/orange/surfaces/detail/SurfaceAction.hh index 5a785e04c7..d1d8dd98e5 100644 --- a/src/orange/surfaces/detail/SurfaceAction.hh +++ b/src/orange/surfaces/detail/SurfaceAction.hh @@ -14,6 +14,7 @@ #include "../PlaneAligned.hh" #include "../CylCentered.hh" #include "../GeneralQuadric.hh" +#include "../Sphere.hh" #include "../SurfaceTypeTraits.hh" #include "../Surfaces.hh" @@ -98,7 +99,9 @@ CELER_FUNCTION auto SurfaceAction::operator()(SurfaceId id) -> result_type ORANGE_SURF_APPLY_IMPL(cy); ORANGE_SURF_APPLY_IMPL(cz); ORANGE_SURF_APPLY_IMPL(p); +#endif ORANGE_SURF_APPLY_IMPL(s); +#if 0 ORANGE_SURF_APPLY_IMPL(kx); ORANGE_SURF_APPLY_IMPL(ky); ORANGE_SURF_APPLY_IMPL(kz); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3648878bf0..f7faf91f72 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -187,6 +187,7 @@ celeritas_add_test(orange/surfaces/detail/QuadraticSolver.test.cc) celeritas_add_test(orange/surfaces/CylCentered.test.cc) celeritas_add_test(orange/surfaces/GeneralQuadric.test.cc) celeritas_add_test(orange/surfaces/PlaneAligned.test.cc) +celeritas_add_test(orange/surfaces/Sphere.test.cc) celeritas_cudaoptional_test(orange/surfaces/SurfaceAction) diff --git a/test/orange/surfaces/GeneralQuadric.test.cc b/test/orange/surfaces/GeneralQuadric.test.cc index a7dfab905c..69b6866132 100644 --- a/test/orange/surfaces/GeneralQuadric.test.cc +++ b/test/orange/surfaces/GeneralQuadric.test.cc @@ -42,6 +42,10 @@ class GeneralQuadricTest : public celeritas::Test */ TEST_F(GeneralQuadricTest, all) { + EXPECT_EQ(celeritas::SurfaceType::gq, GeneralQuadric::surface_type()); + EXPECT_EQ(10, GeneralQuadric::Storage::extent); + EXPECT_EQ(2, GeneralQuadric::Intersections::extent); + const Real3 second{10.3125, 22.9375, 15.75}; const Real3 cross{-21.867141445557, -20.25, 11.69134295109}; const Real3 first{-11.964745962156, -9.1328585544429, -65.69134295109}; diff --git a/test/orange/surfaces/Sphere.test.cc b/test/orange/surfaces/Sphere.test.cc new file mode 100644 index 0000000000..ea3689b6ef --- /dev/null +++ b/test/orange/surfaces/Sphere.test.cc @@ -0,0 +1,95 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file Sphere.test.cc +//---------------------------------------------------------------------------// +#include "orange/surfaces/Sphere.hh" + +#include "base/Constants.hh" +#include "celeritas_test.hh" + +using celeritas::no_intersection; +using celeritas::SignedSense; +using celeritas::Sphere; +using celeritas::SurfaceState; + +using celeritas::ipow; +using celeritas::Real3; +using celeritas::real_type; + +using Intersections = Sphere::Intersections; + +constexpr real_type sqrt_third = 1 / celeritas::constants::sqrt_three; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class SphereTest : public celeritas::Test +{ + protected: + void SetUp() override {} +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// +/*! + * This shape is an ellipsoid: + * - with radii {3, 2, 1}, + * - rotated by 60 degrees about the x axis, + * - then 30 degrees about z, + * - and finally translated by {1, 2, 3}. + */ +TEST_F(SphereTest, all) +{ + EXPECT_EQ(celeritas::SurfaceType::s, Sphere::surface_type()); + EXPECT_EQ(4, Sphere::Storage::extent); + EXPECT_EQ(2, Sphere::Intersections::extent); + + const Real3 origin{-1.1, 2.2, -3.3}; + real_type radius = 4.4; + + Sphere s{origin, radius}; + EXPECT_VEC_SOFT_EQ(origin, s.origin()); + EXPECT_SOFT_EQ(radius * radius, s.radius_sq()); + + EXPECT_EQ(SignedSense::outside, s.calc_sense({4, 5, 5})); + EXPECT_EQ(SignedSense::inside, s.calc_sense({1, 2, -3})); + + const Real3 on_surface{origin[0] + radius * sqrt_third, + origin[1] + radius * sqrt_third, + origin[2] + radius * sqrt_third}; + const Real3 inward{-sqrt_third, -sqrt_third, -sqrt_third}; + const Real3 outward{sqrt_third, sqrt_third, sqrt_third}; + + Intersections distances; + + // On surface, inward + distances = s.calc_intersections(on_surface, inward, SurfaceState::on); + EXPECT_SOFT_EQ(2 * radius, distances[0]); + EXPECT_SOFT_EQ(no_intersection(), distances[1]); + + // "Not on surface", inward + distances = s.calc_intersections(on_surface, inward, SurfaceState::off); + EXPECT_SOFT_EQ(1e-16, distances[0]); + EXPECT_SOFT_EQ(2 * radius, distances[1]); + + // On surface, outward + distances = s.calc_intersections(on_surface, outward, SurfaceState::on); + EXPECT_SOFT_EQ(no_intersection(), distances[0]); + EXPECT_SOFT_EQ(no_intersection(), distances[1]); + + // At center + distances = s.calc_intersections(origin, outward, SurfaceState::off); + EXPECT_SOFT_EQ(radius, distances[1]); + EXPECT_SOFT_EQ(no_intersection(), distances[0]); + + // Outside, hitting both + distances = s.calc_intersections( + Real3{-6.5, 2.2, -3.3}, Real3{1, 0, 0}, SurfaceState::off); + EXPECT_SOFT_EQ(1.0, distances[0]); + EXPECT_SOFT_EQ(1 + 2 * radius, distances[1]); +} diff --git a/test/orange/surfaces/SurfaceAction.test.cc b/test/orange/surfaces/SurfaceAction.test.cc index d0a14bf44d..f2880b6a60 100644 --- a/test/orange/surfaces/SurfaceAction.test.cc +++ b/test/orange/surfaces/SurfaceAction.test.cc @@ -58,8 +58,6 @@ std::string senses_to_string(Span s) }); return result; } - -constexpr real_type noint = celeritas::no_intersection(); } // namespace //---------------------------------------------------------------------------// @@ -81,6 +79,7 @@ class SurfaceActionTest : public celeritas::Test insert(CCylX(5)); insert(CCylY(6)); insert(CCylZ(7)); + insert(Sphere({1, 2, 3}, 1.5)); insert(GeneralQuadric({0, 1, 2}, {3, 4, 5}, {6, 7, 8}, 9)); surf_params_ = SurfaceDataMirror{std::move(surface_data)}; @@ -153,6 +152,7 @@ TEST_F(SurfaceActionTest, string) "Cyl x: r=5", "Cyl y: r=6", "Cyl z: r=7", + "Sphere: r=1.5 at {1,2,3}", "GQuadric: {0,1,2} {3,4,5} {6,7,8} 9"}; // clang-format on EXPECT_VEC_EQ(expected_strings, strings); @@ -179,17 +179,18 @@ TEST_F(SurfaceActionTest, host_distances) // PRINT_EXPECTED(state_ref.distance[test_threads]); const char expected_senses[] - = {'-', '-', '+', '+', '-', '-', '+', '-', '+', '-'}; - const double expected_distance[] = {noint, - noint, - noint, - noint, + = {'-', '-', '+', '+', '-', '-', '+', '+', '-', '-'}; + const double expected_distance[] = {inf, + inf, + inf, + inf, 8.623486582635, 8.115429697208, - noint, - noint, - 22.04764572126, - noint}; + inf, + 5.436749550654, + 0.9761596300109, + 5.848454015622}; + EXPECT_VEC_EQ(expected_senses, senses_to_string(state_ref.sense[test_threads])); EXPECT_VEC_SOFT_EQ(expected_distance, state_ref.distance[test_threads]); @@ -222,17 +223,17 @@ TEST_F(SurfaceActionTest, TEST_IF_CELERITAS_CUDA(device_distances)) auto test_threads = celeritas::range(ThreadId{10}); const char expected_senses[] - = {'-', '-', '+', '+', '-', '-', '+', '-', '+', '-'}; - const double expected_distance[] = {noint, - noint, - noint, - noint, + = {'-', '-', '+', '+', '-', '-', '+', '+', '-', '-'}; + const double expected_distance[] = {inf, + inf, + inf, + inf, 8.623486582635, 8.115429697208, - noint, - noint, - 22.04764572126, - noint}; + inf, + 5.436749550654, + 0.9761596300109, + 5.848454015622}; EXPECT_VEC_EQ(expected_senses, senses_to_string(host_states.sense[test_threads])); EXPECT_VEC_SOFT_EQ(expected_distance, From 4ff31f15ae6db4eefbf04c57a6088aa3b0312db1 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 28 Oct 2021 07:00:17 -0400 Subject: [PATCH 20/22] Apply review feedback --- src/base/Algorithms.hh | 2 +- src/orange/surfaces/CylCentered.hh | 4 ++-- src/orange/surfaces/GeneralQuadric.hh | 4 ++-- src/orange/surfaces/PlaneAligned.hh | 4 ++-- src/orange/surfaces/Sphere.hh | 2 +- test/orange/surfaces/CylCentered.test.cc | 12 ------------ 6 files changed, 8 insertions(+), 20 deletions(-) diff --git a/src/base/Algorithms.hh b/src/base/Algorithms.hh index a60fac304e..bb579a51a0 100644 --- a/src/base/Algorithms.hh +++ b/src/base/Algorithms.hh @@ -82,7 +82,7 @@ CELER_CONSTEXPR_FUNCTION const T& max(const T& a, const T& b) noexcept //---------------------------------------------------------------------------// /*! - * Return the lower of two values. + * Return an iterator to the lowest value in the range. */ template inline CELER_FUNCTION ForwardIt min_element(ForwardIt iter, ForwardIt last) diff --git a/src/orange/surfaces/CylCentered.hh b/src/orange/surfaces/CylCentered.hh index 919b41e09b..df0d415ebc 100644 --- a/src/orange/surfaces/CylCentered.hh +++ b/src/orange/surfaces/CylCentered.hh @@ -71,7 +71,7 @@ class CylCentered // Determine the sense of the position relative to this surface inline CELER_FUNCTION SignedSense calc_sense(const Real3& pos) const; - // Determine the sense of the position relative to this surface + // Calculate all possible straight-line intersections with this surface inline CELER_FUNCTION Intersections calc_intersections( const Real3& pos, const Real3& dir, SurfaceState on_surface) const; @@ -144,7 +144,7 @@ CELER_FUNCTION SignedSense CylCentered::calc_sense(const Real3& pos) const //---------------------------------------------------------------------------// /*! - * Determine the sense of the position relative to this surface. + * Calculate all possible straight-line intersections with this surface. */ template CELER_FUNCTION auto diff --git a/src/orange/surfaces/GeneralQuadric.hh b/src/orange/surfaces/GeneralQuadric.hh index 416d58e2cb..31cc48440e 100644 --- a/src/orange/surfaces/GeneralQuadric.hh +++ b/src/orange/surfaces/GeneralQuadric.hh @@ -79,7 +79,7 @@ class GeneralQuadric // Determine the sense of the position relative to this surface inline CELER_FUNCTION SignedSense calc_sense(const Real3& pos) const; - // Determine the sense of the position relative to this surface + // Calculate all possible straight-line intersections with this surface inline CELER_FUNCTION Intersections calc_intersections( const Real3& pos, const Real3& dir, SurfaceState on_surface) const; @@ -156,7 +156,7 @@ CELER_FUNCTION SignedSense GeneralQuadric::calc_sense(const Real3& pos) const //---------------------------------------------------------------------------// /*! - * Determine the sense of the position relative to this surface. + * Calculate all possible straight-line intersections with this surface. */ CELER_FUNCTION auto GeneralQuadric::calc_intersections(const Real3& pos, diff --git a/src/orange/surfaces/PlaneAligned.hh b/src/orange/surfaces/PlaneAligned.hh index 1de631defc..c26b3bd02d 100644 --- a/src/orange/surfaces/PlaneAligned.hh +++ b/src/orange/surfaces/PlaneAligned.hh @@ -55,7 +55,7 @@ class PlaneAligned // Determine the sense of the position relative to this surface inline CELER_FUNCTION SignedSense calc_sense(const Real3& pos) const; - // Determine the sense of the position relative to this surface + // Calculate all possible straight-line intersections with this surface inline CELER_FUNCTION Intersections calc_intersections( const Real3& pos, const Real3& dir, SurfaceState on_surface) const; @@ -124,7 +124,7 @@ CELER_FUNCTION SignedSense PlaneAligned::calc_sense(const Real3& pos) const //---------------------------------------------------------------------------// /*! - * Determine the sense of the position relative to this surface. + * Calculate all possible straight-line intersections with this surface. */ template CELER_FUNCTION auto diff --git a/src/orange/surfaces/Sphere.hh b/src/orange/surfaces/Sphere.hh index a4e9b20194..8b776f5f3f 100644 --- a/src/orange/surfaces/Sphere.hh +++ b/src/orange/surfaces/Sphere.hh @@ -64,7 +64,7 @@ class Sphere // Determine the sense of the position relative to this surface inline CELER_FUNCTION SignedSense calc_sense(const Real3& pos) const; - // Determine the sense of the position relative to this surface + // Calculate all possible straight-line intersections with this surface inline CELER_FUNCTION Intersections calc_intersections( const Real3& pos, const Real3& dir, SurfaceState on_surface) const; diff --git a/test/orange/surfaces/CylCentered.test.cc b/test/orange/surfaces/CylCentered.test.cc index c5bae73df0..843032983d 100644 --- a/test/orange/surfaces/CylCentered.test.cc +++ b/test/orange/surfaces/CylCentered.test.cc @@ -272,18 +272,6 @@ TEST(TestCCylZ, calc_intersections_on_surface) EXPECT_EQ(no_intersection(), distances[0]); EXPECT_EQ(no_intersection(), distances[1]); - // Tangent - distances = cyl.calc_intersections( - Real3{-1 + eps, 0, 0}, Real3{0, 1, 0}, SurfaceState::on); - EXPECT_EQ(no_intersection(), distances[0]); - EXPECT_EQ(no_intersection(), distances[1]); - - // Heading in, slightly inside - distances = cyl.calc_intersections( - Real3{-1 + eps, 0, 0}, Real3{1, 0, 0}, SurfaceState::on); - EXPECT_SOFT_NEAR(2.0 - eps, distances[0], eps); - EXPECT_EQ(no_intersection(), distances[1]); - // Heading away, slightly outside distances = cyl.calc_intersections( Real3{1.0 - eps, 0, 0}, Real3{1, 0, 0}, SurfaceState::on); From a598fcb6110c20b40bbe6b218052306eb100da38 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 28 Oct 2021 08:26:16 -0400 Subject: [PATCH 21/22] Remove unused types --- src/orange/Types.hh | 31 +------------------------------ 1 file changed, 1 insertion(+), 30 deletions(-) diff --git a/src/orange/Types.hh b/src/orange/Types.hh index 7521eb7614..1d8fb8695a 100644 --- a/src/orange/Types.hh +++ b/src/orange/Types.hh @@ -126,35 +126,6 @@ enum class SurfaceState : bool on = true }; -//---------------------------------------------------------------------------// -// CLASSES -//---------------------------------------------------------------------------// -/*! - * Volume ID and surface ID after initialization. - * - * Possible configurations for the initialization result ('X' means 'has - * a valid ID', i.e. evaluates to true): - * - * Vol | Surface | Description - * :----: | :-----: | :------------------------------- - * | | Failed to find new volume - * X | | Initialized - * X | X | Crossed surface into new volume - * | X | Initialized on a surface (reject) - */ -struct Initialization -{ - VolumeId volume; - SurfaceId surface; - Sense sense = Sense::inside; // Sense if on a surface - - //! Whether initialization succeeded - explicit CELER_FUNCTION operator bool() const - { - return static_cast(volume); - } -}; - //---------------------------------------------------------------------------// // HELPER FUNCTIONS (HOST/DEVICE) //---------------------------------------------------------------------------// @@ -224,7 +195,7 @@ CELER_CONSTEXPR_FUNCTION SurfaceState to_surface_state(SignedSense s) * * \todo There is probably a better place to put this since it's not a "type". * \todo A value of zero might also work since zero-length steps are - * prohibited. + * prohibited. But we'll need custom `min` and `min_element` in that case. */ CELER_CONSTEXPR_FUNCTION real_type no_intersection() { From d311572f0d4a4a6f069a87fc6590c278b4670144 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 28 Oct 2021 11:23:16 -0400 Subject: [PATCH 22/22] fixup! Add general quadric and test placeholder --- src/orange/surfaces/GeneralQuadric.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/orange/surfaces/GeneralQuadric.hh b/src/orange/surfaces/GeneralQuadric.hh index 31cc48440e..330b7e0ea1 100644 --- a/src/orange/surfaces/GeneralQuadric.hh +++ b/src/orange/surfaces/GeneralQuadric.hh @@ -19,7 +19,7 @@ namespace celeritas /*! * General quadric surface. * - * General quadrics that cannot be simplified to other CPTE surfaces include + * General quadrics that cannot be simplified to other ORANGE surfaces include * hyperboloids and paraboloids; and non-axis-aligned cylinders, ellipsoids, * and cones. *