corecel/grid/NonuniformGrid.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 corecel/grid/NonuniformGrid.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Assert.hh"
0011 #include "corecel/Macros.hh"
0012 #include "corecel/Types.hh"
0013 #include "corecel/data/Collection.hh"
0014 #include "corecel/math/Algorithms.hh"
0015
0016 namespace celeritas
0017 {
0018 //---------------------------------------------------------------------------//
0019 /*!
0020  * Interact with a nonuniform grid of increasing values.
0021  *
0022  * This should have the same interface (aside from constructor) as
0023  * UniformGrid.
0024  */
0025 template<class T>
0026 class NonuniformGrid
0027 {
0028   public:
0029     //!@{
0030     //! \name Type aliases
0031     using value_type = T;
0032     using Storage
0033         = Collection<value_type, Ownership::const_reference, MemSpace::native>;
0034     using ItemRangeT = ItemRange<value_type>;
0035     //!@}
0036
0037   public:
0038     // Construct with storage
0039     inline CELER_FUNCTION
0040     NonuniformGrid(ItemRangeT const& values, Storage const& storage);
0041
0042     //! Number of grid points
0043     CELER_FORCEINLINE_FUNCTION size_type size() const
0044     {
0045         return offset_.size();
0046     }
0047
0048     //! Minimum/first value
0049     CELER_FORCEINLINE_FUNCTION value_type front() const
0050     {
0051         return storage_[*offset_.begin()];
0052     }
0053
0054     //! Maximum/last value
0055     CELER_FORCEINLINE_FUNCTION value_type back() const
0056     {
0057         return storage_[*(offset_.end() - 1)];
0058     }
0059
0060     // Calculate the value at the given grid point
0061     inline CELER_FUNCTION value_type operator[](size_type i) const;
0062
0063     // Find the index of the given value (*must* be in bounds)
0064     inline CELER_FUNCTION size_type find(value_type value) const;
0065
0066     //! Low-level access to offsets for downstream utilities
0067     CELER_FORCEINLINE_FUNCTION ItemRangeT offset() const { return offset_; }
0068
0069   private:
0070     Storage const& storage_;
0071     ItemRangeT offset_;
0072 };
0073
0074 //---------------------------------------------------------------------------//
0075 // INLINE DEFINITIONS
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Construct with a range indexing into backend storage.
0079  */
0080 template<class T>
0081 CELER_FUNCTION NonuniformGrid<T>::NonuniformGrid(ItemRangeT const& values,
0082                                                  Storage const& storage)
0083     : storage_{storage}, offset_{values}
0084 {
0085     CELER_EXPECT(offset_.size() >= 2);
0086     CELER_EXPECT(*offset_.end() <= storage.size());
0087     CELER_EXPECT(this->front() <= this->back());  // Approximation for "sorted"
0088 }
0089
0090 //---------------------------------------------------------------------------//
0091 /*!
0092  * Get the value at the given grid point.
0093  */
0094 template<class T>
0095 CELER_FUNCTION auto
0096 NonuniformGrid<T>::operator[](size_type i) const -> value_type
0097 {
0098     CELER_EXPECT(i < offset_.size());
0099     return storage_[offset_[i]];
0100 }
0101
0102 //---------------------------------------------------------------------------//
0103 /*!
0104  * Find the value bin such that storage[result] <= value < data[result + 1].
0105  *
0106  * The given value *must* be in range, because out-of-bounds values usually
0107  * require different treatment (e.g. clipping to the boundary values rather
0108  * than interpolating). It's easier to test the exceptional cases (final grid
0109  * point) outside of the grid view.
0110  */
0111 template<class T>
0112 CELER_FUNCTION size_type NonuniformGrid<T>::find(value_type value) const
0113 {
0114     CELER_EXPECT(value >= this->front() && value < this->back());
0115
0116     using ItemIdT = ItemId<T>;
0117     auto iter = celeritas::lower_bound(
0118         offset_.begin(),
0119         offset_.end(),
0120         value,
0121         [&v = storage_](ItemIdT i, T value) { return v[i] < value; });
0122     CELER_ASSERT(iter != offset_.end());
0123
0124     if (value != storage_[*iter])
0125     {
0126         // Exactly on end grid point, or not on a grid point at all: move to
0127         // previous bin
0128         --iter;
0129     }
0130
0131     return iter - offset_.begin();
0132 }
0133
0134 //---------------------------------------------------------------------------//
0135 }  // namespace celeritas