optical/surface/SurfacePhysicsUtils.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 celeritas/optical/surface/SurfacePhysicsUtils.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "corecel/math/ArrayUtils.hh"
0011 #include "celeritas/optical/Types.hh"
0012
0013 #if !CELER_DEVICE_COMPILE
0014 #    include "corecel/io/Logger.hh"
0015 #    include "corecel/io/Repr.hh"
0016 #endif
0017
0018 namespace celeritas
0019 {
0020 namespace optical
0021 {
0022 //---------------------------------------------------------------------------//
0023 /*!
0024  * Whether a track is entering the surface defined by the given normal.
0025  *
0026  * The surface normal convention used in Celeritas optical physics is that
0027  * the normal direction points opposite the incident track direction. This
0028  * function makes checks for this condition explicit in the code.
0029  */
0030 inline CELER_FUNCTION bool is_entering_surface(Real3 const& dir,
0031                                                Real3 const& normal)
0032 {
0033     return dot_product(dir, normal) < 0;
0034 }
0035
0036 //---------------------------------------------------------------------------//
0037 /*!
0038  * Get the local surface ID given the current local material and direction.
0039  */
0040 inline CELER_FUNCTION LocalSurfaceId local_surf_id(LocalPositionId pos,
0041                                                    LocalDirection dir)
0042 {
0043     CELER_EXPECT(
0044         pos && (pos > LocalPositionId{0} || dir == LocalDirection::forward));
0045     return id_cast<LocalSurfaceId>(
0046         static_cast<int>(pos.unchecked_get()) + static_cast<int>(dir) - 1);
0047 }
0048
0049 //---------------------------------------------------------------------------//
0050 /*!
0051  * Get the next track surface position in the given direction.
0052  *
0053  * Type-safe operation to ensure direction is only added in track-local frames.
0054  * Uses unsigned underflow when moving reverse (dir = -1) while on a
0055  * pre-surface (pos = 0) to wrap to an invalid position value.
0056  */
0057 inline CELER_FUNCTION LocalPositionId next_local_pos_id(LocalPositionId pos,
0058                                                         LocalDirection dir)
0059 {
0060     CELER_EXPECT(pos);
0061     return LocalPositionId(pos.unchecked_get() + to_signed_offset(dir));
0062 }
0063
0064 //---------------------------------------------------------------------------//
0065 /*!
0066  * Calculate subsurface direction from a track's geometry direction.
0067  */
0068 inline CELER_FUNCTION LocalDirection calc_subsurface_direction(
0069     Real3 const& geo_dir, Real3 const& normal)
0070 {
0071     return static_cast<LocalDirection>(is_entering_surface(geo_dir, normal));
0072 }
0073
0074 //---------------------------------------------------------------------------//
0075 /*!
0076  * Calculate geometric reflection of an incident vector about a normal.
0077  */
0078 [[nodiscard]] inline CELER_FUNCTION Real3 geometric_reflected_from(
0079     Real3 dir, Real3 const& normal)
0080 {
0081     axpy(-2 * dot_product(dir, normal), normal, &dir);
0082     return dir;
0083 }
0084
0085 //---------------------------------------------------------------------------//
0086 /*!
0087  * Sample a valid facet normal by wrapping a roughness sampler.
0088  *
0089  * Some facet normal calculators might not produce surface normals valid for
0090  * optical physics surface crossings (see \c is_entering_surface ). This
0091  * functor will construct and repeatedly sample the distribution until the
0092  * track is exiting the sampled facet normal.
0093  *
0094  * \par Example
0095  *
0096  * This uses a gaussian roughness distribution to accept only facets that have
0097  * a direction where the track is exiting the sampled normal.
0098  *
0099  * \code
0100    EnteringSurfaceNormalSampler<GaussianRoughnessSampler> sample_facet{
0101        track.geometry().dir(),
0102        s_phys.global_normal(),
0103        data.sigma_alpha[sub_model_id]};
0104    auto rng = track.rng();
0105    s_phys.facet_normal(sample_facet(rng));
0106  * \endcode
0107  */
0108 template<class Calculator>
0109 class EnteringSurfaceNormalSampler
0110 {
0111   public:
0112     template<class... Args>
0113     CELER_FUNCTION EnteringSurfaceNormalSampler(
0114         Real3 const& dir, Real3 const& normal, Args&&... args)
0115         : dir_{dir}, sample_normal_{normal, celeritas::forward<Args>(args)...}
0116     {
0117         CELER_EXPECT(is_entering_surface(normal, dir));
0118     }
0119
0120     // Repeatedly sample facet normal until satisfies entering surface
0121     template<class Engine>
0122     CELER_FUNCTION Real3 operator()(Engine& rng)
0123     {
0124         Real3 local_normal;
0125         int loop_guard{256};
0126         do
0127         {
0128             --loop_guard;
0129 #if !CELER_DEVICE_COMPILE
0130             CELER_VALIDATE(
0131                 loop_guard > 0,
0132                 << "failed to sample a microfacet direction into which "
0133                 << repr(dir_) << " is entering (last normal sampled: "
0134                 << repr(local_normal) << ")");
0135 #else
0136             if (CELER_UNLIKELY(loop_guard == 0))
0137                 return {0, 0, 0};
0138 #endif
0139             local_normal = sample_normal_(rng);
0140         } while (!is_entering_surface(local_normal, dir_));
0141         return local_normal;
0142     }
0143
0144   private:
0145     Real3 const& dir_;
0146     Calculator sample_normal_;
0147 };
0148
0149 //---------------------------------------------------------------------------//
0150 }  // namespace optical
0151 }  // namespace celeritas