univ/detail/SurfaceFunctors.hh

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