univ/detail/SurfaceFunctors.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file orange/univ/detail/SurfaceFunctors.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/cont/Array.hh"
0012 #include "corecel/math/Algorithms.hh"
0013 #include "corecel/math/ArrayUtils.hh"
0014 #include "corecel/math/NumericLimits.hh"
0015
0016 #include "Types.hh"
0017
0018 namespace celeritas
0019 {
0020 namespace detail
0021 {
0022 //---------------------------------------------------------------------------//
0023 //! Calculate the sense of a surface at a given position.
0024 struct CalcSense
0025 {
0026     Real3 const& pos;
0027
0028     template<class S>
0029     CELER_FUNCTION SignedSense operator()(S const& surf)
0030     {
0031         return surf.calc_sense(this->pos);
0032     }
0033 };
0034
0035 //---------------------------------------------------------------------------//
0036 //! Get the number of intersections of a surface
0037 struct NumIntersections
0038 {
0039     template<class ST>
0040     CELER_CONSTEXPR_FUNCTION size_type operator()(ST) const noexcept
0041     {
0042         using S = typename ST::type;
0043         return typename S::Intersections{}.size();
0044     }
0045 };
0046
0047 //---------------------------------------------------------------------------//
0048 //! Calculate the outward normal at a position.
0049 struct CalcNormal
0050 {
0051     Real3 const& pos;
0052
0053     template<class S>
0054     CELER_FUNCTION Real3 operator()(S const& surf)
0055     {
0056         return surf.calc_normal(this->pos);
0057     }
0058 };
0059
0060 //---------------------------------------------------------------------------//
0061 /*!
0062  * Calculate the smallest distance from a point to the surface.
0063  *
0064  * For certain surface types (spheres, cylinders, planes), defined such that
0065  * the normal is *outward* (positive when "outside", negative when "inside"),
0066  * the nearest distance to the surface can be calculated quite trivially.
0067  */
0068 struct CalcSafetyDistance
0069 {
0070     Real3 const& pos;
0071
0072     //! Operate on a surface
0073     template<class S>
0074     CELER_FUNCTION real_type operator()(S const& surf)
0075     {
0076         if (!S::simple_safety())
0077         {
0078             // Not a surface that satisfies our simplifying constraints: return
0079             // a conservative answer.
0080             return 0;
0081         }
0082
0083         // Calculate outward normal
0084         Real3 dir = surf.calc_normal(this->pos);
0085         if (CELER_UNLIKELY(std::isnan(dir[0])))
0086         {
0087             // Magnitude may have been zero, e.g. taking the safety at the
0088             // center of a sphere/cyl, so assume it's a long way off.
0089             return numeric_limits<real_type>::infinity();
0090         }
0091         CELER_ASSERT(is_soft_unit_vector(dir));
0092
0093         auto sense = surf.calc_sense(this->pos);
0094         // If sense is "positive" (on or outside), flip direction to inward so
0095         // that the vector points toward the surface
0096         if (sense == SignedSense::outside)
0097         {
0098             for (real_type& d : dir)
0099             {
0100                 d *= -1;
0101             }
0102         }
0103         else if (sense == SignedSense::on)
0104         {
0105             return 0;
0106         }
0107
0108         // Return the closest intersection
0109         auto intersect
0110             = surf.calc_intersections(this->pos, dir, SurfaceState::off);
0111         return *celeritas::min_element(intersect.begin(), intersect.end());
0112     }
0113 };
0114
0115 //---------------------------------------------------------------------------//
0116 /*!
0117  * Fill an array with valid distances-to-intersection.
0118  *
0119  * \tparam F Predicate for returning whether the distance is allowable
0120  *
0121  * This assumes that each call is to the next face index, starting with face
0122  * zero.
0123  */
0124 template<class F>
0125 class CalcIntersections
0126 {
0127   public:
0128     //! Construct from the particle point, direction, face ID, and temp storage
0129     CELER_FUNCTION CalcIntersections(F is_valid_isect,
0130                                      Real3 const& pos,
0131                                      Real3 const& dir,
0132                                      FaceId on_face,
0133                                      bool is_simple,
0134                                      TempNextFace const& next_face)
0135         : is_valid_isect_(celeritas::forward<F>(is_valid_isect))
0136         , pos_(pos)
0137         , dir_(dir)
0138         , on_face_idx_(on_face.unchecked_get())
0139         , fill_isect_(!is_simple)
0140         , face_(next_face.face)
0141         , distance_(next_face.distance)
0142         , isect_(next_face.isect)
0143     {
0144         CELER_EXPECT(face_ && distance_);
0145     }
0146
0147     //! Operate on a surface
0148     template<class S>
0149     CELER_FUNCTION void operator()(S const& surf)
0150     {
0151         auto on_surface = (on_face_idx_ == face_idx_) ? SurfaceState::on
0152                                                       : SurfaceState::off;
0153         if constexpr (typename S::Intersections{}.size() == 1)
0154         {
0155             if (on_surface == SurfaceState::on)
0156             {
0157                 // On surface so cannot reintersect
0158                 ++face_idx_;
0159                 return;
0160             }
0161         }
0162
0163         // Calculate distance to surface along this direction
0164         auto all_dist = surf.calc_intersections(pos_, dir_, on_surface);
0165
0166         // Copy possible intersections and this surface to the output
0167         for (real_type dist : all_dist)
0168         {
0169             CELER_ASSERT(dist > 0);
0170             if (is_valid_isect_(dist))
0171             {
0172                 // Save intersection in the list
0173                 face_[isect_idx_] = FaceId{face_idx_};
0174                 distance_[isect_idx_] = dist;
0175                 if (fill_isect_)
0176                 {
0177                     isect_[isect_idx_] = isect_idx_;
0178                 }
0179                 ++isect_idx_;
0180             }
0181         }
0182         // Increment to next face
0183         ++face_idx_;
0184     }
0185
0186     CELER_FUNCTION size_type face_idx() const { return face_idx_; }
0187     CELER_FUNCTION size_type isect_idx() const { return isect_idx_; }
0188
0189   private:
0190     //// DATA ////
0191
0192     F is_valid_isect_;
0193     Real3 const& pos_;
0194     Real3 const& dir_;
0195     size_type const on_face_idx_;
0196     bool const fill_isect_;
0197     FaceId* const face_;
0198     real_type* const distance_;
0199     size_type* const isect_;
0200     size_type face_idx_{0};
0201     size_type isect_idx_{0};
0202 };
0203
0204 template<class F, class... Args>
0205 CELER_FUNCTION CalcIntersections(F&&, Args&&... args) -> CalcIntersections<F>;
0206
0207 //---------------------------------------------------------------------------//
0208 }  // namespace detail
0209 }  // namespace celeritas