Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-16 09:03:19

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/OrangeTrackView.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008 
0009 #include "corecel/Macros.hh"
0010 #include "corecel/Types.hh"
0011 #include "corecel/cont/Array.hh"
0012 #include "corecel/sys/ThreadId.hh"
0013 
0014 #include "LevelStateAccessor.hh"
0015 #include "OrangeData.hh"
0016 #include "OrangeTypes.hh"
0017 #include "transform/TransformVisitor.hh"
0018 #include "univ/SimpleUnitTracker.hh"
0019 #include "univ/TrackerVisitor.hh"
0020 #include "univ/UniverseTypeTraits.hh"
0021 #include "univ/detail/Types.hh"
0022 
0023 #include "detail/UniverseIndexer.hh"
0024 
0025 #if !CELER_DEVICE_COMPILE
0026 #    include "corecel/io/Logger.hh"
0027 #    include "corecel/io/Repr.hh"
0028 #endif
0029 
0030 namespace celeritas
0031 {
0032 //---------------------------------------------------------------------------//
0033 /*!
0034  * Navigate through an ORANGE geometry on a single thread.
0035  *
0036  * Since the navigation relies on computationally expensive calls and must
0037  * ensure a consistent state between physical and logical boundaries, there is
0038  * an ordering followed by Celeritas' internal calls to each track's state.
0039  * Access (\c pos, \c dir, \c volume/surface/is_outside/is_on_boundary) is
0040  * valid at any time.
0041  *
0042  * The required ordering is:
0043  *
0044  * - Initialization, via the assignment operator
0045  * - Locating the boundary crossing along the current direction, \c
0046  *   find_next_step
0047  * - Locating the closest point on the boundary in any direction, \c
0048  *   find_safety
0049  * - Movement within a volume, not crossing a boundary: \c move_internal or \c
0050  *   move_to_boundary
0051  * - If on a boundary, logically moving to the adjacent volume ("relocation")
0052  *   with \c cross_boundary
0053  *
0054  * At any time, \c set_dir may be called, but then \c find_next_step must again
0055  * be called before any subsequent \c move or \c cross action above .
0056  *
0057  * The main point is that \c find_next_step depends on the current
0058  * straight-line direction, \c move_to_boundary and \c move_internal (with
0059  * a step length) depends on that distance, and
0060  * \c cross_boundary depends on being on the boundary with a knowledge of the
0061  * post-boundary state.
0062  *
0063  * \todo \c move_internal with a position \em should depend on the safety
0064  * distance, but that check is not yet implemented.
0065  */
0066 class OrangeTrackView
0067 {
0068   public:
0069     //!@{
0070     //! \name Type aliases
0071     using ParamsRef = NativeCRef<OrangeParamsData>;
0072     using StateRef = NativeRef<OrangeStateData>;
0073     using Initializer_t = GeoTrackInitializer;
0074     //!@}
0075 
0076     //! Helper struct for initializing from an existing geometry state
0077     struct DetailedInitializer
0078     {
0079         OrangeTrackView const& other;  //!< Existing geometry
0080         Real3 const& dir;  //!< New direction
0081     };
0082 
0083   public:
0084     // Construct from params and state
0085     inline CELER_FUNCTION OrangeTrackView(ParamsRef const& params,
0086                                           StateRef const& states,
0087                                           TrackSlotId tid);
0088 
0089     // Initialize the state
0090     inline CELER_FUNCTION OrangeTrackView& operator=(Initializer_t const& init);
0091     // Initialize the state from a parent state and new direction
0092     inline CELER_FUNCTION OrangeTrackView&
0093     operator=(DetailedInitializer const& init);
0094 
0095     //// ACCESSORS ////
0096 
0097     // The current position
0098     inline CELER_FUNCTION Real3 const& pos() const;
0099     // The current direction
0100     inline CELER_FUNCTION Real3 const& dir() const;
0101 
0102     // The current volume ID (null if outside)
0103     inline CELER_FUNCTION VolumeId volume_id() const;
0104     // Get the physical volume ID in the current cell
0105     inline CELER_FUNCTION VolumeInstanceId volume_instance_id() const;
0106     // The current level
0107     inline CELER_FUNCTION LevelId const& level() const;
0108     // Get the volume instance ID for all levels
0109     inline CELER_FUNCTION void volume_instance_id(Span<VolumeInstanceId>) const;
0110 
0111     // The current surface ID
0112     inline CELER_FUNCTION SurfaceId surface_id() const;
0113     // After 'find_next_step', the next straight-line surface
0114     inline CELER_FUNCTION SurfaceId next_surface_id() const;
0115     // Whether the track is outside the valid geometry region
0116     inline CELER_FUNCTION bool is_outside() const;
0117     // Whether the track is exactly on a surface
0118     inline CELER_FUNCTION bool is_on_boundary() const;
0119     //! Whether the last operation resulted in an error
0120     CELER_FORCEINLINE_FUNCTION bool failed() const { return failed_; }
0121 
0122     //// OPERATIONS ////
0123 
0124     // Find the distance to the next boundary
0125     inline CELER_FUNCTION Propagation find_next_step();
0126 
0127     // Find the distance to the next boundary, up to and including a step
0128     inline CELER_FUNCTION Propagation find_next_step(real_type max_step);
0129 
0130     // Find the distance to the nearest boundary in any direction
0131     inline CELER_FUNCTION real_type find_safety();
0132 
0133     // Find the distance to the nearest nearby boundary in any direction
0134     inline CELER_FUNCTION real_type find_safety(real_type max_step);
0135 
0136     // Move to the boundary in preparation for crossing it
0137     inline CELER_FUNCTION void move_to_boundary();
0138 
0139     // Move within the volume
0140     inline CELER_FUNCTION void move_internal(real_type step);
0141 
0142     // Move within the volume to a specific point
0143     inline CELER_FUNCTION void move_internal(Real3 const& pos);
0144 
0145     // Cross from one side of the current surface to the other
0146     inline CELER_FUNCTION void cross_boundary();
0147 
0148     // Change direction
0149     inline CELER_FUNCTION void set_dir(Real3 const& newdir);
0150 
0151   private:
0152     //// TYPES ////
0153 
0154     using LSA = LevelStateAccessor;
0155 
0156     //// DATA ////
0157 
0158     ParamsRef const& params_;
0159     StateRef const& states_;
0160     TrackSlotId track_slot_;
0161     bool failed_{false};
0162 
0163     //// PRIVATE STATE MUTATORS ////
0164 
0165     // The current level
0166     inline CELER_FUNCTION void level(LevelId);
0167 
0168     // The boundary on the current surface level
0169     inline CELER_FUNCTION void boundary(BoundaryResult);
0170 
0171     // The next step distance, as stored on the state
0172     inline CELER_FUNCTION void next_step(real_type dist);
0173 
0174     // The next surface to be encounted
0175     inline CELER_FUNCTION void next_surf(detail::OnLocalSurface const&);
0176 
0177     // The level of the next surface to be encounted
0178     inline CELER_FUNCTION void next_surface_level(LevelId);
0179 
0180     //// PRIVATE STATE ACCESSORS ////
0181 
0182     // The current surface level
0183     inline CELER_FUNCTION LevelId const& surface_level() const;
0184 
0185     // The local surface on the current surface level
0186     inline CELER_FUNCTION LocalSurfaceId const& surf() const;
0187 
0188     // The sense on the current surface level
0189     inline CELER_FUNCTION Sense const& sense() const;
0190 
0191     // The boundary on the current surface level
0192     inline CELER_FUNCTION BoundaryResult const& boundary() const;
0193 
0194     // The next step distance, as stored on the state
0195     inline CELER_FUNCTION real_type const& next_step() const;
0196 
0197     // The next surface to be encounted
0198     inline CELER_FUNCTION detail::OnLocalSurface next_surf() const;
0199 
0200     // The level of the next surface to be encounted
0201     inline CELER_FUNCTION LevelId const& next_surface_level() const;
0202 
0203     //// HELPER FUNCTIONS ////
0204 
0205     // Iterate over lower levels to find the next step
0206     inline CELER_FUNCTION Propagation
0207     find_next_step_impl(detail::Intersection isect);
0208 
0209     // Create local sense reference
0210     inline CELER_FUNCTION Span<SenseValue> make_temp_sense() const;
0211 
0212     // Create local distance
0213     inline CELER_FUNCTION detail::TempNextFace make_temp_next() const;
0214 
0215     inline CELER_FUNCTION detail::LocalState
0216     make_local_state(LevelId level) const;
0217 
0218     // Whether the next distance-to-boundary has been found
0219     inline CELER_FUNCTION bool has_next_step() const;
0220 
0221     // Whether the next surface has been found
0222     inline CELER_FUNCTION bool has_next_surface() const;
0223 
0224     // Invalidate the next distance-to-boundary and surface
0225     inline CELER_FUNCTION void clear_next();
0226 
0227     // Assign the surface on the current level
0228     inline CELER_FUNCTION void
0229     surface(LevelId level, detail::OnLocalSurface surf);
0230 
0231     // Clear the surface on the current level
0232     inline CELER_FUNCTION void clear_surface();
0233 
0234     // Make a LevelStateAccessor for the current thread and level
0235     inline CELER_FUNCTION LSA make_lsa() const;
0236 
0237     // Make a LevelStateAccessor for the current thread and a given level
0238     inline CELER_FUNCTION LSA make_lsa(LevelId level) const;
0239 
0240     // Get the daughter ID for the volume in the universe (or null)
0241     inline CELER_FUNCTION DaughterId get_daughter(LSA const& lsa);
0242 
0243     // Get the transform ID for the given daughter.
0244     inline CELER_FUNCTION TransformId get_transform(DaughterId daughter_id);
0245 
0246     // Get the transform ID to move from this level to the one below.
0247     inline CELER_FUNCTION TransformId get_transform(LevelId lev);
0248 };
0249 
0250 //---------------------------------------------------------------------------//
0251 // MEMBER FUNCTIONS
0252 //---------------------------------------------------------------------------//
0253 /*!
0254  * Construct from persistent and state data.
0255  */
0256 CELER_FUNCTION
0257 OrangeTrackView::OrangeTrackView(ParamsRef const& params,
0258                                  StateRef const& states,
0259                                  TrackSlotId tid)
0260     : params_(params), states_(states), track_slot_(tid)
0261 {
0262     CELER_EXPECT(params_);
0263     CELER_EXPECT(states_);
0264     CELER_EXPECT(track_slot_ < states.size());
0265 }
0266 
0267 //---------------------------------------------------------------------------//
0268 /*!
0269  * Construct the state.
0270  *
0271  * Expensive. This function should only be called to initialize an event from a
0272  * starting location and direction. Secondaries will initialize their states
0273  * from a copy of the parent.
0274  */
0275 CELER_FUNCTION OrangeTrackView&
0276 OrangeTrackView::operator=(Initializer_t const& init)
0277 {
0278     CELER_EXPECT(is_soft_unit_vector(init.dir));
0279 
0280     failed_ = false;
0281 
0282     // Create local state
0283     detail::LocalState local;
0284     local.pos = init.pos;
0285     local.dir = init.dir;
0286     local.volume = {};
0287     local.surface = {};
0288     local.temp_sense = this->make_temp_sense();
0289 
0290     // Helpers for applying parent-to-daughter transformations
0291     TransformVisitor apply_transform{params_};
0292     auto transform_down_local = [&local](auto&& t) {
0293         local.pos = t.transform_down(local.pos);
0294         local.dir = t.rotate_down(local.dir);
0295     };
0296 
0297     // Recurse into daughter universes starting with the outermost universe
0298     UniverseId univ_id = top_universe_id();
0299     DaughterId daughter_id;
0300     LevelId level{0};
0301     do
0302     {
0303         TrackerVisitor visit_tracker{params_};
0304         auto tinit = visit_tracker(
0305             [&local](auto&& t) { return t.initialize(local); }, univ_id);
0306 
0307         if (!tinit.volume || tinit.surface)
0308         {
0309 #if !CELER_DEVICE_COMPILE
0310             auto msg = CELER_LOG_LOCAL(error);
0311             msg << "Failed to initialize geometry state: ";
0312             if (!tinit.volume)
0313             {
0314                 msg << "could not find associated volume";
0315             }
0316             else
0317             {
0318                 msg << "started on a surface ("
0319                     << tinit.surface.id().unchecked_get() << ")";
0320             }
0321             msg << " in universe " << univ_id.unchecked_get()
0322                 << " at local position " << repr(local.pos);
0323 #endif
0324             // Mark as failed and place in local "exterior" to end the search
0325             // but preserve the current level information
0326             failed_ = true;
0327             tinit.volume = orange_exterior_volume;
0328         }
0329 
0330         auto lsa = this->make_lsa(level);
0331         lsa.vol() = tinit.volume;
0332         lsa.pos() = local.pos;
0333         lsa.dir() = local.dir;
0334         lsa.universe() = univ_id;
0335 
0336         daughter_id = visit_tracker(
0337             [&tinit](auto&& t) { return t.daughter(tinit.volume); }, univ_id);
0338 
0339         if (daughter_id)
0340         {
0341             auto const& daughter = params_.daughters[daughter_id];
0342             // Apply "transform down" based on stored transform
0343             apply_transform(transform_down_local, daughter.trans_id);
0344             // Update universe and increase level depth
0345             univ_id = daughter.universe_id;
0346             ++level;
0347         }
0348 
0349     } while (daughter_id);
0350 
0351     // Save found level
0352     this->level(level);
0353 
0354     // Reset surface/boundary information
0355     this->boundary(BoundaryResult::exiting);
0356     this->clear_surface();
0357     this->clear_next();
0358 
0359     CELER_ENSURE(!this->has_next_step());
0360     return *this;
0361 }
0362 
0363 //---------------------------------------------------------------------------//
0364 /*!
0365  * Construct the state from a direction and a copy of the parent state.
0366  */
0367 CELER_FUNCTION
0368 OrangeTrackView& OrangeTrackView::operator=(DetailedInitializer const& init)
0369 {
0370     CELER_EXPECT(is_soft_unit_vector(init.dir));
0371 
0372     failed_ = false;
0373 
0374     if (this != &init.other)
0375     {
0376         // Copy init track's position and logical state
0377         this->level(states_.level[init.other.track_slot_]);
0378         this->surface(init.other.surface_level(),
0379                       {init.other.surf(), init.other.sense()});
0380         this->boundary(init.other.boundary());
0381 
0382         for (auto lev : range(LevelId{this->level() + 1}))
0383         {
0384             // Copy all data accessed via LSA
0385             auto lsa = this->make_lsa(lev);
0386             lsa = init.other.make_lsa(lev);
0387         }
0388     }
0389 
0390     // Clear the next step information since we're changing direction or
0391     // initializing a new state
0392     this->clear_next();
0393 
0394     // Transform direction from global to local
0395     Real3 localdir = init.dir;
0396     auto apply_transform = TransformVisitor{params_};
0397     auto rotate_down
0398         = [&localdir](auto&& t) { localdir = t.rotate_down(localdir); };
0399     for (auto lev : range(LevelId{this->level()}))
0400     {
0401         auto lsa = this->make_lsa(lev);
0402         lsa.dir() = localdir;
0403         apply_transform(rotate_down,
0404                         this->get_transform(this->get_daughter(lsa)));
0405     }
0406 
0407     // Save final level
0408     this->make_lsa().dir() = localdir;
0409 
0410     CELER_ENSURE(!this->has_next_step());
0411     return *this;
0412 }
0413 
0414 //---------------------------------------------------------------------------//
0415 /*!
0416  * The current position.
0417  */
0418 CELER_FUNCTION Real3 const& OrangeTrackView::pos() const
0419 {
0420     return this->make_lsa(LevelId{0}).pos();
0421 }
0422 
0423 //---------------------------------------------------------------------------//
0424 /*!
0425  * The current direction.
0426  */
0427 CELER_FUNCTION Real3 const& OrangeTrackView::dir() const
0428 {
0429     return this->make_lsa(LevelId{0}).dir();
0430 }
0431 
0432 //---------------------------------------------------------------------------//
0433 /*!
0434  * The current volume ID.
0435  *
0436  * \note It is allowable to call this function when "outside", because the
0437  * outside in ORANGE is just a special volume. Other geometries may not have
0438  * that behavior.
0439  */
0440 CELER_FUNCTION VolumeId OrangeTrackView::volume_id() const
0441 {
0442     auto lsa = this->make_lsa();
0443     detail::UniverseIndexer ui(params_.universe_indexer_data);
0444     return ui.global_volume(lsa.universe(), lsa.vol());
0445 }
0446 
0447 //---------------------------------------------------------------------------//
0448 /*!
0449  * The current volume instance.
0450  *
0451  * \todo not implemented; VolumeId is already halfway between a "reusable
0452  * volume" and a "volume instance" anyway...
0453  */
0454 CELER_FUNCTION VolumeInstanceId OrangeTrackView::volume_instance_id() const
0455 {
0456     return {};
0457 }
0458 
0459 //---------------------------------------------------------------------------//
0460 /*!
0461  * The current level.
0462  */
0463 CELER_FORCEINLINE_FUNCTION LevelId const& OrangeTrackView::level() const
0464 {
0465     return states_.level[track_slot_];
0466 }
0467 
0468 //---------------------------------------------------------------------------//
0469 /*!
0470  * Get the volume instance ID at every level.
0471  *
0472  * The input span size must be equal to the value of "level" plus one. The
0473  * top-most level ("world" or level zero) starts at index zero and moves
0474  * downward. Note that Geant4 uses the \em reverse nomenclature.
0475  */
0476 CELER_FUNCTION void
0477 OrangeTrackView::volume_instance_id(Span<VolumeInstanceId> levels) const
0478 {
0479     CELER_EXPECT(levels.size() == this->level().get() + 1);
0480     for (auto lev : range(levels.size()))
0481     {
0482         levels[lev] = {};
0483     }
0484 }
0485 
0486 //---------------------------------------------------------------------------//
0487 /*!
0488  * The current surface ID.
0489  */
0490 CELER_FUNCTION SurfaceId OrangeTrackView::surface_id() const
0491 {
0492     if (this->is_on_boundary())
0493     {
0494         auto lsa = this->make_lsa(this->surface_level());
0495         detail::UniverseIndexer ui{params_.universe_indexer_data};
0496         return ui.global_surface(lsa.universe(), this->surf());
0497     }
0498     else
0499     {
0500         return SurfaceId{};
0501     }
0502 }
0503 
0504 //---------------------------------------------------------------------------//
0505 /*!
0506  * After 'find_next_step', the next straight-line surface.
0507  */
0508 CELER_FUNCTION SurfaceId OrangeTrackView::next_surface_id() const
0509 {
0510     CELER_EXPECT(this->has_next_surface());
0511     auto lsa = this->make_lsa(this->next_surface_level());
0512     detail::UniverseIndexer ui{params_.universe_indexer_data};
0513     return ui.global_surface(lsa.universe(), this->next_surf().id());
0514 }
0515 
0516 //---------------------------------------------------------------------------//
0517 /*!
0518  * Whether the track is outside the valid geometry region.
0519  */
0520 CELER_FUNCTION bool OrangeTrackView::is_outside() const
0521 {
0522     // Zeroth volume in outermost universe is always the exterior by
0523     // construction in ORANGE
0524     auto lsa = this->make_lsa(LevelId{0});
0525     return lsa.vol() == orange_exterior_volume;
0526 }
0527 
0528 //---------------------------------------------------------------------------//
0529 /*!
0530  * Whether the track is exactly on a surface.
0531  */
0532 CELER_FORCEINLINE_FUNCTION bool OrangeTrackView::is_on_boundary() const
0533 {
0534     return static_cast<bool>(this->surface_level());
0535 }
0536 
0537 //---------------------------------------------------------------------------//
0538 /*!
0539  * Find the distance to the next geometric boundary.
0540  */
0541 CELER_FUNCTION Propagation OrangeTrackView::find_next_step()
0542 {
0543     if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
0544     {
0545         // On a boundary, headed back in: next step is zero
0546         return {0, true};
0547     }
0548 
0549     // Find intersection at the top level: always the first simple unit
0550     auto global_isect = [this] {
0551         SimpleUnitTracker t{params_, SimpleUnitId{0}};
0552         return t.intersect(this->make_local_state(LevelId{0}));
0553     }();
0554     // Find intersection for all lower levels
0555     return this->find_next_step_impl(global_isect);
0556 }
0557 
0558 //---------------------------------------------------------------------------//
0559 /*!
0560  * Find a nearby distance to the next geometric boundary up to a distance.
0561  *
0562  * This may reduce the number of surfaces needed to check, sort, or write to
0563  * temporary memory, thereby speeding up transport.
0564  */
0565 CELER_FUNCTION Propagation OrangeTrackView::find_next_step(real_type max_step)
0566 {
0567     CELER_EXPECT(max_step > 0);
0568 
0569     if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
0570     {
0571         // On a boundary, headed back in: next step is zero
0572         return {0, true};
0573     }
0574 
0575     // Find intersection at the top level: always the first simple unit
0576     auto global_isect = [this, &max_step] {
0577         SimpleUnitTracker t{params_, SimpleUnitId{0}};
0578         return t.intersect(this->make_local_state(LevelId{0}), max_step);
0579     }();
0580     // Find intersection for all lower levels
0581     auto result = this->find_next_step_impl(global_isect);
0582     CELER_ENSURE(result.distance <= max_step);
0583     return result;
0584 }
0585 
0586 //---------------------------------------------------------------------------//
0587 /*!
0588  * Move to the next straight-line boundary but do not change volume.
0589  */
0590 CELER_FUNCTION void OrangeTrackView::move_to_boundary()
0591 {
0592     CELER_EXPECT(this->boundary() != BoundaryResult::reentrant);
0593     CELER_EXPECT(this->has_next_step());
0594     CELER_EXPECT(this->has_next_surface());
0595 
0596     // Physically move next step
0597     real_type const dist = this->next_step();
0598     for (auto i : range(this->level() + 1))
0599     {
0600         auto lsa = this->make_lsa(LevelId{i});
0601         axpy(dist, lsa.dir(), &lsa.pos());
0602     }
0603 
0604     this->surface(this->next_surface_level(), this->next_surf());
0605     this->clear_next();
0606     CELER_ENSURE(this->is_on_boundary());
0607 }
0608 
0609 //---------------------------------------------------------------------------//
0610 /*!
0611  * Move within the current volume.
0612  *
0613  * The straight-line distance *must* be less than the distance to the
0614  * boundary.
0615  */
0616 CELER_FUNCTION void OrangeTrackView::move_internal(real_type dist)
0617 {
0618     CELER_EXPECT(this->has_next_step());
0619     CELER_EXPECT(dist > 0 && dist <= this->next_step());
0620     CELER_EXPECT(dist != this->next_step() || !this->has_next_surface());
0621 
0622     // Move and update the next step
0623     for (auto i : range(this->level() + 1))
0624     {
0625         auto lsa = this->make_lsa(LevelId{i});
0626         axpy(dist, lsa.dir(), &lsa.pos());
0627     }
0628     this->next_step(this->next_step() - dist);
0629     this->clear_surface();
0630 }
0631 
0632 //---------------------------------------------------------------------------//
0633 /*!
0634  * Move within the current volume to a nearby point.
0635  *
0636  * \todo Currently it's up to the caller to make sure that the position is
0637  * "nearby". We should actually test this with an "is inside" call.
0638  */
0639 CELER_FUNCTION void OrangeTrackView::move_internal(Real3 const& pos)
0640 {
0641     // Transform all nonlocal levels
0642     auto local_pos = pos;
0643     auto apply_transform = TransformVisitor{params_};
0644     auto translate_down
0645         = [&local_pos](auto&& t) { local_pos = t.transform_down(local_pos); };
0646     for (auto lev : range(LevelId{this->level()}))
0647     {
0648         auto lsa = this->make_lsa(lev);
0649         lsa.pos() = local_pos;
0650 
0651         // Apply "transform down" based on stored transform
0652         apply_transform(translate_down,
0653                         this->get_transform(this->get_daughter(lsa)));
0654     }
0655 
0656     // Save final level
0657     auto lsa = this->make_lsa();
0658     lsa.pos() = local_pos;
0659 
0660     // Clear surface state and next-step info
0661     this->clear_surface();
0662     this->clear_next();
0663 }
0664 
0665 //---------------------------------------------------------------------------//
0666 /*!
0667  * Cross from one side of the current surface to the other.
0668  *
0669  * The position *must* be on the boundary following a move-to-boundary. This
0670  * should only be called once per boundary crossing.
0671  */
0672 CELER_FUNCTION void OrangeTrackView::cross_boundary()
0673 {
0674     CELER_EXPECT(this->is_on_boundary());
0675     CELER_EXPECT(!this->has_next_step());
0676 
0677     if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
0678     {
0679         // Direction changed while on boundary leading to no change in
0680         // volume/surface. This is logically equivalent to a reflection.
0681         this->boundary(BoundaryResult::exiting);
0682         return;
0683     }
0684 
0685     // Cross surface by flipping the sense
0686     states_.sense[track_slot_] = flip_sense(this->sense());
0687     this->boundary(BoundaryResult::exiting);
0688 
0689     // Create local state from post-crossing level and updated sense
0690     LevelId level{this->surface_level()};
0691     UniverseId universe;
0692     LocalVolumeId volume;
0693     detail::LocalState local;
0694     {
0695         auto lsa = this->make_lsa(level);
0696         universe = lsa.universe();
0697         local.pos = lsa.pos();
0698         local.dir = lsa.dir();
0699         local.volume = lsa.vol();
0700         local.surface = {this->surf(), this->sense()};
0701         local.temp_sense = this->make_temp_sense();
0702     }
0703 
0704     TrackerVisitor visit_tracker{params_};
0705 
0706     // Update the post-crossing volume by crossing the boundary of the "surface
0707     // crossing" level universe
0708     volume = visit_tracker(
0709         [&local](auto&& t) { return t.cross_boundary(local).volume; },
0710         universe);
0711     if (CELER_UNLIKELY(!volume))
0712     {
0713         // Boundary crossing failure
0714 #if !CELER_DEVICE_COMPILE
0715         CELER_LOG_LOCAL(error)
0716             << "track failed to cross local surface "
0717             << this->surf().unchecked_get() << " in universe "
0718             << universe.unchecked_get() << " at local position "
0719             << repr(local.pos) << " along local direction " << repr(local.dir);
0720 #endif
0721         // Mark as failed and place in local "exterior" to end the search
0722         // but preserve the current level information
0723         failed_ = true;
0724         volume = orange_exterior_volume;
0725     }
0726     make_lsa(level).vol() = volume;
0727 
0728     // Clear local surface before diving into daughters
0729     // TODO: this is where we'd do inter-universe mapping
0730     local.volume = {};
0731     local.surface = {};
0732 
0733     // Starting with the current level (i.e., next_surface_level), iterate
0734     // down into the deepest level: *initializing* not *crossing*
0735     auto daughter_id = visit_tracker(
0736         [volume](auto&& t) { return t.daughter(volume); }, universe);
0737     while (daughter_id)
0738     {
0739         ++level;
0740         {
0741             // Update universe, local position/direction
0742             auto const& daughter = params_.daughters[daughter_id];
0743             TransformVisitor apply_transform{params_};
0744             auto transform_down_local = [&local](auto&& t) {
0745                 local.pos = t.transform_down(local.pos);
0746                 local.dir = t.rotate_down(local.dir);
0747             };
0748             apply_transform(transform_down_local, daughter.trans_id);
0749             universe = daughter.universe_id;
0750         }
0751 
0752         // Initialize in daughter and get IDs of volume and potential daughter
0753         volume = visit_tracker(
0754             [&local](auto&& t) { return t.initialize(local).volume; },
0755             universe);
0756 
0757         if (!volume)
0758         {
0759 #if !CELER_DEVICE_COMPILE
0760             auto msg = CELER_LOG_LOCAL(error);
0761             msg << "track failed to cross boundary: could not find associated "
0762                    "volume in universe "
0763                 << universe.unchecked_get() << " at local position "
0764                 << repr(local.pos);
0765 #endif
0766             // Mark as failed and place in local "exterior" to end the search
0767             // but preserve the current level information
0768             failed_ = true;
0769             volume = orange_exterior_volume;
0770         }
0771         daughter_id = visit_tracker(
0772             [volume](auto&& t) { return t.daughter(volume); }, universe);
0773 
0774         auto lsa = make_lsa(level);
0775         lsa.vol() = volume;
0776         lsa.pos() = local.pos;
0777         lsa.dir() = local.dir;
0778         lsa.universe() = universe;
0779     }
0780 
0781     // Save final level
0782     this->level(level);
0783 
0784     CELER_ENSURE(this->is_on_boundary());
0785 }
0786 
0787 //---------------------------------------------------------------------------//
0788 /*!
0789  * Change the track's direction.
0790  *
0791  * This happens after a scattering event or movement inside a magnetic field.
0792  * It resets the calculated distance-to-boundary. It is allowed to happen on
0793  * the boundary, but changing direction so that it goes from pointing outward
0794  * to inward (or vice versa) will mean that \c cross_boundary will be a
0795  * null-op.
0796  *
0797  * TODO: This needs to be updated to handle reflections through levels
0798  */
0799 CELER_FUNCTION void OrangeTrackView::set_dir(Real3 const& newdir)
0800 {
0801     CELER_EXPECT(is_soft_unit_vector(newdir));
0802 
0803     if (this->is_on_boundary())
0804     {
0805         // Changing direction on a boundary, which may result in not leaving
0806         // current volume upon the cross_surface call
0807         auto lsa = this->make_lsa(this->surface_level());
0808 
0809         TrackerVisitor visit_tracker{params_};
0810         auto normal = visit_tracker(
0811             [pos = lsa.pos(), local_surface = this->surf()](auto&& t) {
0812                 return t.normal(pos, local_surface);
0813             },
0814             lsa.universe());
0815 
0816         // Normal is in *local* coordinates but newdir is in *global*: rotate
0817         // up to check
0818         auto apply_transform = TransformVisitor{params_};
0819         auto rotate_up = [&normal](auto&& t) { normal = t.rotate_up(normal); };
0820         for (auto level : range<int>(this->level().unchecked_get()).step(-1))
0821         {
0822             apply_transform(rotate_up, this->get_transform(LevelId(level)));
0823         }
0824 
0825         // Evaluate whether the direction dotted with the surface normal
0826         // changes (i.e. heading from inside to outside or vice versa).
0827         if ((dot_product(normal, newdir) >= 0)
0828             != (dot_product(normal, this->dir()) >= 0))
0829         {
0830             // The boundary crossing direction has changed! Reverse our
0831             // plans to change the logical state and move to a new volume.
0832             this->boundary(flip_boundary(this->boundary()));
0833         }
0834     }
0835 
0836     // Complete direction setting by transforming direction all the way down
0837     Real3 localdir = newdir;
0838     auto apply_transform = TransformVisitor{params_};
0839     auto rotate_down
0840         = [&localdir](auto&& t) { localdir = t.rotate_down(localdir); };
0841     for (auto lev : range(this->level()))
0842     {
0843         auto lsa = this->make_lsa(lev);
0844         lsa.dir() = localdir;
0845         apply_transform(rotate_down,
0846                         this->get_transform(this->get_daughter(lsa)));
0847     }
0848     // Save final level
0849     this->make_lsa().dir() = localdir;
0850 
0851     this->clear_next();
0852 }
0853 
0854 //---------------------------------------------------------------------------//
0855 // STATE ACCESSORS
0856 //---------------------------------------------------------------------------//
0857 /*!
0858  * The current level.
0859  */
0860 CELER_FORCEINLINE_FUNCTION void OrangeTrackView::level(LevelId lev)
0861 {
0862     states_.level[track_slot_] = lev;
0863 }
0864 
0865 /*!
0866  * The boundary on the current surface level.
0867  */
0868 CELER_FORCEINLINE_FUNCTION void OrangeTrackView::boundary(BoundaryResult br)
0869 {
0870     states_.boundary[track_slot_] = br;
0871 }
0872 
0873 /*!
0874  * The next step distance.
0875  */
0876 CELER_FORCEINLINE_FUNCTION void OrangeTrackView::next_step(real_type dist)
0877 {
0878     states_.next_step[track_slot_] = dist;
0879 }
0880 
0881 /*!
0882  * The next surface to be encountered.
0883  */
0884 CELER_FORCEINLINE_FUNCTION void
0885 OrangeTrackView::next_surf(detail::OnLocalSurface const& s)
0886 {
0887     states_.next_surf[track_slot_] = s.id();
0888     states_.next_sense[track_slot_] = s.unchecked_sense();
0889 }
0890 
0891 /*!
0892  * The level of the next surface to be encounted.
0893  */
0894 CELER_FORCEINLINE_FUNCTION void OrangeTrackView::next_surface_level(LevelId lev)
0895 {
0896     states_.next_level[track_slot_] = lev;
0897 }
0898 
0899 //---------------------------------------------------------------------------//
0900 // CONST STATE ACCESSORS
0901 /*!
0902  * The current surface level.
0903  */
0904 CELER_FORCEINLINE_FUNCTION LevelId const& OrangeTrackView::surface_level() const
0905 {
0906     return states_.surface_level[track_slot_];
0907 }
0908 
0909 /*!
0910  * The local surface on the current surface level.
0911  */
0912 CELER_FORCEINLINE_FUNCTION LocalSurfaceId const& OrangeTrackView::surf() const
0913 {
0914     return states_.surf[track_slot_];
0915 }
0916 
0917 /*!
0918  * The sense on the current surface level.
0919  */
0920 CELER_FORCEINLINE_FUNCTION Sense const& OrangeTrackView::sense() const
0921 {
0922     return states_.sense[track_slot_];
0923 }
0924 
0925 /*!
0926  * The boundary on the current surface level.
0927  */
0928 CELER_FORCEINLINE_FUNCTION BoundaryResult const&
0929 OrangeTrackView::boundary() const
0930 {
0931     return states_.boundary[track_slot_];
0932 }
0933 
0934 /*!
0935  * The next step distance.
0936  */
0937 CELER_FORCEINLINE_FUNCTION real_type const& OrangeTrackView::next_step() const
0938 {
0939     return states_.next_step[track_slot_];
0940 }
0941 
0942 /*!
0943  * The next surface to be encountered.
0944  */
0945 CELER_FORCEINLINE_FUNCTION detail::OnLocalSurface
0946 OrangeTrackView::next_surf() const
0947 {
0948     return {states_.next_surf[track_slot_], states_.next_sense[track_slot_]};
0949 }
0950 
0951 /*!
0952  * The level of the next surface to be encounted.
0953  */
0954 CELER_FORCEINLINE_FUNCTION LevelId const&
0955 OrangeTrackView::next_surface_level() const
0956 {
0957     return states_.next_level[track_slot_];
0958 }
0959 
0960 //---------------------------------------------------------------------------//
0961 // HELPER FUNCTIONS
0962 //---------------------------------------------------------------------------//
0963 /*!
0964  * Iterate over levels 1 to N to find the next step.
0965  *
0966  * Caller is responsible for finding the candidate next step on level 0, and
0967  * passing the resultant Intersection object as an argument.
0968  */
0969 CELER_FUNCTION Propagation
0970 OrangeTrackView::find_next_step_impl(detail::Intersection isect)
0971 {
0972     TrackerVisitor visit_tracker{params_};
0973 
0974     // The level with minimum distance to intersection
0975     LevelId min_level{0};
0976 
0977     // Find the nearest intersection from level 0 to current level
0978     // inclusive, prefering the shallowest level (i.e., lowest univ_id)
0979     for (auto levelid : range(LevelId{1}, this->level() + 1))
0980     {
0981         auto univ_id = this->make_lsa(levelid).universe();
0982         auto local_isect = visit_tracker(
0983             [local_state = this->make_local_state(levelid), &isect](auto&& t) {
0984                 return t.intersect(local_state, isect.distance);
0985             },
0986             univ_id);
0987 
0988         if (local_isect.distance < isect.distance)
0989         {
0990             isect = local_isect;
0991             min_level = levelid;
0992         }
0993     }
0994 
0995     this->next_step(isect.distance);
0996     this->next_surf(isect.surface);
0997     if (isect)
0998     {
0999         // Save level corresponding to the intersection
1000         this->next_surface_level(min_level);
1001     }
1002 
1003     Propagation result;
1004     result.distance = isect.distance;
1005     result.boundary = static_cast<bool>(isect);
1006     return result;
1007 }
1008 
1009 //---------------------------------------------------------------------------//
1010 /*!
1011  * Find the distance to the nearest boundary in any direction.
1012  *
1013  * The safety distance at a given point is the minimum safety distance over all
1014  * levels, since surface deduplication can potentionally elide bounding
1015  * surfaces at more deeply embedded levels.
1016  */
1017 CELER_FUNCTION real_type OrangeTrackView::find_safety()
1018 {
1019     CELER_EXPECT(!this->is_on_boundary());
1020 
1021     TrackerVisitor visit_tracker{params_};
1022 
1023     real_type min_safety_dist = numeric_limits<real_type>::infinity();
1024 
1025     for (auto lev : range(LevelId{this->level() + 1}))
1026     {
1027         auto lsa = this->make_lsa(lev);
1028         auto sd = visit_tracker(
1029             [&lsa](auto&& t) { return t.safety(lsa.pos(), lsa.vol()); },
1030             lsa.universe());
1031         min_safety_dist = celeritas::min(min_safety_dist, sd);
1032     }
1033     return min_safety_dist;
1034 }
1035 
1036 //---------------------------------------------------------------------------//
1037 /*!
1038  * Find the distance to the nearest nearby boundary.
1039  *
1040  * Since we currently support only "simple" safety distances, we can't
1041  * eliminate anything by checking only nearby surfaces.
1042  */
1043 CELER_FUNCTION real_type OrangeTrackView::find_safety(real_type)
1044 {
1045     return this->find_safety();
1046 }
1047 
1048 //---------------------------------------------------------------------------//
1049 /*!
1050  * Get a reference to the current volume, or to world volume if outside.
1051  */
1052 CELER_FUNCTION Span<SenseValue> OrangeTrackView::make_temp_sense() const
1053 {
1054     auto const max_faces = params_.scalars.max_faces;
1055     auto offset = track_slot_.get() * max_faces;
1056     return states_.temp_sense[AllItems<SenseValue, MemSpace::native>{}].subspan(
1057         offset, max_faces);
1058 }
1059 
1060 //---------------------------------------------------------------------------//
1061 /*!
1062  * Set up intersection scratch space.
1063  */
1064 CELER_FUNCTION detail::TempNextFace OrangeTrackView::make_temp_next() const
1065 {
1066     auto const max_isect = params_.scalars.max_intersections;
1067     auto offset = track_slot_.get() * max_isect;
1068 
1069     detail::TempNextFace result;
1070     result.face = states_.temp_face[AllItems<FaceId>{}].data() + offset;
1071     result.distance = states_.temp_distance[AllItems<real_type>{}].data()
1072                       + offset;
1073     result.isect = states_.temp_isect[AllItems<size_type>{}].data() + offset;
1074     result.size = max_isect;
1075 
1076     return result;
1077 }
1078 
1079 //---------------------------------------------------------------------------//
1080 /*!
1081  * Create a local state.
1082  */
1083 CELER_FUNCTION detail::LocalState
1084 OrangeTrackView::make_local_state(LevelId level) const
1085 {
1086     detail::LocalState local;
1087 
1088     auto lsa = this->make_lsa(level);
1089 
1090     local.pos = lsa.pos();
1091     local.dir = lsa.dir();
1092     local.volume = lsa.vol();
1093     if (level == this->surface_level())
1094     {
1095         local.surface = {this->surf(), this->sense()};
1096     }
1097     else
1098     {
1099         local.surface = {};
1100     }
1101     local.temp_sense = this->make_temp_sense();
1102     local.temp_next = this->make_temp_next();
1103     return local;
1104 }
1105 
1106 //---------------------------------------------------------------------------//
1107 /*!
1108  * Whether any next step has been calculated.
1109  */
1110 CELER_FORCEINLINE_FUNCTION bool OrangeTrackView::has_next_step() const
1111 {
1112     return this->next_step() != 0;
1113 }
1114 
1115 //---------------------------------------------------------------------------//
1116 /*!
1117  * Whether the next intersecting surface has been found.
1118  */
1119 CELER_FORCEINLINE_FUNCTION bool OrangeTrackView::has_next_surface() const
1120 {
1121     return static_cast<bool>(states_.next_surf[track_slot_]);
1122 }
1123 
1124 //---------------------------------------------------------------------------//
1125 /*!
1126  * Reset the next distance-to-boundary and surface.
1127  */
1128 CELER_FUNCTION void OrangeTrackView::clear_next()
1129 {
1130     this->next_step(0);
1131     states_.next_surf[track_slot_] = {};
1132     CELER_ENSURE(!this->has_next_step() && !this->has_next_surface());
1133 }
1134 
1135 //---------------------------------------------------------------------------//
1136 /*!
1137  * Assign the surface on the current level.
1138  */
1139 CELER_FUNCTION void
1140 OrangeTrackView::surface(LevelId level, detail::OnLocalSurface surf)
1141 {
1142     states_.surface_level[track_slot_] = level;
1143     states_.surf[track_slot_] = surf.id();
1144     states_.sense[track_slot_] = surf.unchecked_sense();
1145 }
1146 
1147 //---------------------------------------------------------------------------//
1148 /*!
1149  * Clear the surface on the current level.
1150  */
1151 CELER_FUNCTION void OrangeTrackView::clear_surface()
1152 {
1153     states_.surface_level[track_slot_] = {};
1154     CELER_ENSURE(!this->is_on_boundary());
1155 }
1156 
1157 //---------------------------------------------------------------------------//
1158 /*!
1159  * Make a LevelStateAccessor for the current thread and level.
1160  */
1161 CELER_FORCEINLINE_FUNCTION auto OrangeTrackView::make_lsa() const -> LSA
1162 {
1163     return this->make_lsa(this->level());
1164 }
1165 
1166 //---------------------------------------------------------------------------//
1167 /*!
1168  * Make a LevelStateAccessor for the current thread and a given level.
1169  */
1170 CELER_FORCEINLINE_FUNCTION auto
1171 OrangeTrackView::make_lsa(LevelId level) const -> LSA
1172 {
1173     return LSA(&states_, track_slot_, level);
1174 }
1175 
1176 //---------------------------------------------------------------------------//
1177 /*!
1178  * Get the daughter ID for the given volume in the given universe.
1179  *
1180  * \return DaughterId or {} if the current volume is a leaf.
1181  */
1182 CELER_FUNCTION DaughterId OrangeTrackView::get_daughter(LSA const& lsa)
1183 {
1184     TrackerVisitor visit_tracker{params_};
1185     return visit_tracker([&lsa](auto&& t) { return t.daughter(lsa.vol()); },
1186                          lsa.universe());
1187 }
1188 
1189 //---------------------------------------------------------------------------//
1190 /*!
1191  * Get the transform ID for the given daughter.
1192  */
1193 CELER_FUNCTION TransformId OrangeTrackView::get_transform(DaughterId daughter_id)
1194 {
1195     CELER_EXPECT(daughter_id);
1196     return params_.daughters[daughter_id].trans_id;
1197 }
1198 
1199 //---------------------------------------------------------------------------//
1200 /*!
1201  * Get the transform ID for the given daughter.
1202  */
1203 CELER_FUNCTION TransformId OrangeTrackView::get_transform(LevelId lev)
1204 {
1205     CELER_EXPECT(lev < this->level());
1206     LSA lsa(&states_, track_slot_, lev);
1207     return this->get_transform(this->get_daughter(lsa));
1208 }
1209 
1210 //---------------------------------------------------------------------------//
1211 }  // namespace celeritas