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/CMakeLists.txt b/src/CMakeLists.txt index 9330cd4a11..c1b507e4aa 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -61,6 +61,9 @@ list(APPEND SOURCES comm/ScopedMpiInit.cc 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 io/ImportPhysicsVector.cc diff --git a/src/base/Algorithms.hh b/src/base/Algorithms.hh index ae0374887d..bb579a51a0 100644 --- a/src/base/Algorithms.hh +++ b/src/base/Algorithms.hh @@ -80,6 +80,26 @@ CELER_CONSTEXPR_FUNCTION const T& max(const T& a, const T& b) noexcept return (b > a) ? b : a; } +//---------------------------------------------------------------------------// +/*! + * Return an iterator to the lowest value in the range. + */ +template +inline CELER_FUNCTION ForwardIt min_element(ForwardIt iter, ForwardIt last) +{ + // Avoid incrementing past the end + if (iter == last) + return last; + + ForwardIt result = iter++; + for (; iter != last; ++iter) + { + if (*iter < *result) + result = iter; + } + return result; +} + //---------------------------------------------------------------------------// // Replace/extend //---------------------------------------------------------------------------// 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 //---------------------------------------------------------------------------// 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/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/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) 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/base/Span.hh b/src/base/Span.hh index ad017489d5..993e5b77dd 100644 --- a/src/base/Span.hh +++ b/src/base/Span.hh @@ -145,6 +145,9 @@ class Span detail::SpanImpl s_; }; +template +constexpr 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/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/Data.hh b/src/orange/Data.hh new file mode 100644 index 0000000000..036a02e78f --- /dev/null +++ b/src/orange/Data.hh @@ -0,0 +1,63 @@ +//----------------------------------*-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/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; + } +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/Types.cc b/src/orange/Types.cc new file mode 100644 index 0000000000..ff1af8e2af --- /dev/null +++ b/src/orange/Types.cc @@ -0,0 +1,54 @@ +//----------------------------------*-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" + +#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", +#endif + "s", +#if 0 + "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 new file mode 100644 index 0000000000..1d8fb8695a --- /dev/null +++ b/src/orange/Types.hh @@ -0,0 +1,245 @@ +//----------------------------------*-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; + +//---------------------------------------------------------------------------// +// ENUMERATIONS +//---------------------------------------------------------------------------// +/*! + * Whether a position is logically "inside" or "outside" a surface. + * + * For a plane, "outside" (true) is the "positive" sense and equivalent to + * \f[ + \vec x \cdot \vec n >= 0 + * \f] + * 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. + */ +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, //!< X axis/I index coordinate + y, //!< Y axis/J index coordinate + z, //!< Z axis/K index coordinate + size_ //!< Sentinel value for looping over axes +}; + +//---------------------------------------------------------------------------// +/*! + * 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 +#endif + s, //!< Sphere +#if 0 + 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_ //!< Sentinel value for number of surface types +}; + +//---------------------------------------------------------------------------// +/*! + * Evaluated quadric expression allowing for distinct 'on surface' state. + * + * 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 +}; + +//---------------------------------------------------------------------------// +/*! + * 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 +}; + +//---------------------------------------------------------------------------// +// 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); +} + +//---------------------------------------------------------------------------// +/*! + * 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". + * + * \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. But we'll need custom `min` and `min_element` in that case. + */ +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/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 new file mode 100644 index 0000000000..df0d415ebc --- /dev/null +++ b/src/orange/surfaces/CylCentered.hh @@ -0,0 +1,222 @@ +//----------------------------------*-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/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 Storage = Span; + //@} + + //// CLASS ATTRIBUTES //// + + // Surface type identifier + static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type(); + + public: + //// CONSTRUCTORS //// + + // Construct with radius + explicit inline CELER_FUNCTION CylCentered(real_type radius); + + // Construct from raw data + explicit inline CELER_FUNCTION CylCentered(Storage); + + //// 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 Storage 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; + + // 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; + + // 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(Storage 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_); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate all possible straight-line intersections with 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/GeneralQuadric.hh b/src/orange/surfaces/GeneralQuadric.hh new file mode 100644 index 0000000000..330b7e0ea1 --- /dev/null +++ b/src/orange/surfaces/GeneralQuadric.hh @@ -0,0 +1,206 @@ +//----------------------------------*-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 ORANGE 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 Storage = Span; + using SpanConstReal3 = Span; + //@} + + //// CLASS ATTRIBUTES //// + + //! Surface type identifier + static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type() + { + return SurfaceType::gq; + } + + 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(Storage); + + //// 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 Storage 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; + + // 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; + + // 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(Storage 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); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate all possible straight-line intersections with 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 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, b / 2, 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 new file mode 100644 index 0000000000..c26b3bd02d --- /dev/null +++ b/src/orange/surfaces/PlaneAligned.hh @@ -0,0 +1,161 @@ +//----------------------------------*-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/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 Storage = Span; + //@} + + //// CLASS ATTRIBUTES //// + + // Surface type identifier + static CELER_CONSTEXPR_FUNCTION SurfaceType surface_type(); + + public: + //// CONSTRUCTORS //// + + // Construct with radius + explicit inline CELER_FUNCTION PlaneAligned(real_type position); + + // Construct from raw data + explicit inline CELER_FUNCTION PlaneAligned(Storage); + + //// 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 Storage 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; + + // 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; + + // 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 from axis intercept. + */ +template +CELER_FUNCTION PlaneAligned::PlaneAligned(real_type position) + : position_(position) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Construct from raw data. + */ +template +CELER_FUNCTION PlaneAligned::PlaneAligned(Storage 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_); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate all possible straight-line intersections with 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&) const +{ + Real3 norm{0, 0, 0}; + + norm[t_index()] = 1.; + return norm; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/orange/surfaces/Sphere.hh b/src/orange/surfaces/Sphere.hh new file mode 100644 index 0000000000..8b776f5f3f --- /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; + + // 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; + + // 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/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/SurfaceIO.cc b/src/orange/surfaces/SurfaceIO.cc new file mode 100644 index 0000000000..5b03e3ce02 --- /dev/null +++ b/src/orange/surfaces/SurfaceIO.cc @@ -0,0 +1,61 @@ +//----------------------------------*-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" +#include "Sphere.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); +//---------------------------------------------------------------------------// +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 new file mode 100644 index 0000000000..d38a064be6 --- /dev/null +++ b/src/orange/surfaces/SurfaceIO.hh @@ -0,0 +1,30 @@ +//----------------------------------*-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 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 new file mode 100644 index 0000000000..22965f9bdd --- /dev/null +++ b/src/orange/surfaces/SurfaceTypeTraits.hh @@ -0,0 +1,65 @@ +//----------------------------------*-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; +class Sphere; + +//---------------------------------------------------------------------------// +/*! + * 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); +#endif +ORANGE_SURFACE_TRAITS(s, Sphere); +#if 0 +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 diff --git a/src/orange/surfaces/Surfaces.hh b/src/orange/surfaces/Surfaces.hh new file mode 100644 index 0000000000..115a1ec5b3 --- /dev/null +++ b/src/orange/surfaces/Surfaces.hh @@ -0,0 +1,97 @@ +//----------------------------------*-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 +{ +//---------------------------------------------------------------------------// +/*! + * Access stored surface data. + * + * These are all surfaces in the entire geometry. + */ +class Surfaces +{ + public: + //@{ + //! Type aliases + using SurfaceDataRef + = SurfaceData; + //@} + + public: + // 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; + + // 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 SurfaceDataRef& data_; +}; + +//---------------------------------------------------------------------------// +// INLINE FUNCTION DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with reference to persistent data. + */ +CELER_FUNCTION Surfaces::Surfaces(const SurfaceDataRef& 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()); + CELER_EXPECT(this->surface_type(sid) == T::surface_type()); + + 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/QuadraticSolver.hh b/src/orange/surfaces/detail/QuadraticSolver.hh new file mode 100644 index 0000000000..fdbd9086c5 --- /dev/null +++ b/src/orange/surfaces/detail/QuadraticSolver.hh @@ -0,0 +1,194 @@ +//----------------------------------*-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, SurfaceState 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, + SurfaceState on_surface) + -> Intersections +{ + if (std::fabs(a) >= min_a()) + { + // Not along the surface + QuadraticSolver solve(a, half_b); + return on_surface == SurfaceState::on ? solve() : solve(c); + } + else + { + // Travelling parallel to the quadric's surface + if (on_surface == SurfaceState::off) + { + 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/src/orange/surfaces/detail/SurfaceAction.hh b/src/orange/surfaces/detail/SurfaceAction.hh new file mode 100644 index 0000000000..d1d8dd98e5 --- /dev/null +++ b/src/orange/surfaces/detail/SurfaceAction.hh @@ -0,0 +1,134 @@ +//----------------------------------*-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 "../Sphere.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); +#endif + ORANGE_SURF_APPLY_IMPL(s); +#if 0 + 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 ecabd09e8a..f7faf91f72 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -178,6 +178,19 @@ 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) +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) + #-----------------------------------------------------------------------------# # I/O (ROOT) diff --git a/test/base/Algorithms.test.cc b/test/base/Algorithms.test.cc index e9ecfe1a1f..14e0058afc 100644 --- a/test/base/Algorithms.test.cc +++ b/test/base/Algorithms.test.cc @@ -61,3 +61,23 @@ 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, min_element_idx()); + + v = {100}; + EXPECT_EQ(0, min_element_idx()); + + v = {10, 2, 100, 3, -1}; + EXPECT_EQ(4, min_element_idx()); + + v[2] = -100; + EXPECT_EQ(2, min_element_idx()); +} 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 diff --git a/test/orange/surfaces/CylCentered.test.cc b/test/orange/surfaces/CylCentered.test.cc new file mode 100644 index 0000000000..843032983d --- /dev/null +++ b/test/orange/surfaces/CylCentered.test.cc @@ -0,0 +1,443 @@ +//----------------------------------*-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::Storage::extent); + 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]); + + // 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/GeneralQuadric.test.cc b/test/orange/surfaces/GeneralQuadric.test.cc new file mode 100644 index 0000000000..69b6866132 --- /dev/null +++ b/test/orange/surfaces/GeneralQuadric.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 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) +{ + 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}; + 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]); +} diff --git a/test/orange/surfaces/PlaneAligned.test.cc b/test/orange/surfaces/PlaneAligned.test.cc new file mode 100644 index 0000000000..13b6ec8cb6 --- /dev/null +++ b/test/orange/surfaces/PlaneAligned.test.cc @@ -0,0 +1,129 @@ +//----------------------------------*-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.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, -1.1, 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.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)); +} 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 new file mode 100644 index 0000000000..f2880b6a60 --- /dev/null +++ b/test/orange/surfaces/SurfaceAction.test.cc @@ -0,0 +1,242 @@ +//----------------------------------*-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 +#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 "random/distributions/IsotropicDistribution.hh" +#include "random/distributions/UniformBoxDistribution.hh" +#include "celeritas_test.hh" +#include "SurfaceAction.test.hh" + +using namespace celeritas; +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; +} +} // namespace + +//---------------------------------------------------------------------------// +// 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(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)}; + } + + 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_; +}; + +//---------------------------------------------------------------------------// +// 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", + "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); +} + +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[] = {inf, + inf, + inf, + inf, + 8.623486582635, + 8.115429697208, + 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]); +} + +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[] = {inf, + inf, + inf, + inf, + 8.623486582635, + 8.115429697208, + 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, + 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..a6ac1db96d --- /dev/null +++ b/test/orange/surfaces/SurfaceAction.test.hh @@ -0,0 +1,168 @@ +//----------------------------------*-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 "base/CollectionBuilder.hh" +#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/orange/surfaces/detail/QuadraticSolver.test.cc b/test/orange/surfaces/detail/QuadraticSolver.test.cc new file mode 100644 index 0000000000..e56c69fdf6 --- /dev/null +++ b/test/orange/surfaces/detail/QuadraticSolver.test.cc @@ -0,0 +1,164 @@ +//----------------------------------*-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::SurfaceState; +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, SurfaceState::off); + + 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, SurfaceState::off); + + EXPECT_SOFT_EQ(no_intersection(), x[0]); + EXPECT_SOFT_NEAR(1.0e3, x[1], 1e-7); +} 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