corecel/grid/NonuniformGrid.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 corecel/grid/NonuniformGrid.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Macros.hh"
0011 #include "corecel/Types.hh"
0012 #include "corecel/data/Collection.hh"
0013 #include "corecel/math/Algorithms.hh"
0014
0015 namespace celeritas
0016 {
0017 //---------------------------------------------------------------------------//
0018 /*!
0019  * Interact with a nonuniform grid of increasing values.
0020  *
0021  * This should have the same interface (aside from constructor) as
0022  * UniformGrid.
0023  */
0024 template<class T>
0025 class NonuniformGrid
0026 {
0027   public:
0028     //!@{
0029     //! \name Type aliases
0030     using value_type = T;
0031     using Storage
0032         = Collection<value_type, Ownership::const_reference, MemSpace::native>;
0033     using ItemRangeT = ItemRange<value_type>;
0034     using SpanConstT = typename Storage::SpanConstT;
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     // Construct a span referring to the grid points
0070     inline CELER_FUNCTION SpanConstT values() const;
0071
0072   private:
0073     Storage const& storage_;
0074     ItemRangeT offset_;
0075 };
0076
0077 //---------------------------------------------------------------------------//
0078 // INLINE DEFINITIONS
0079 //---------------------------------------------------------------------------//
0080 /*!
0081  * Construct with a range indexing into backend storage.
0082  */
0083 template<class T>
0084 CELER_FUNCTION NonuniformGrid<T>::NonuniformGrid(ItemRangeT const& values,
0085                                                  Storage const& storage)
0086     : storage_{storage}, offset_{values}
0087 {
0088     CELER_EXPECT(offset_.size() >= 2);
0089     CELER_EXPECT(*offset_.end() <= storage.size());
0090     CELER_EXPECT(this->front() <= this->back());  // Approximation for "sorted"
0091 }
0092
0093 //---------------------------------------------------------------------------//
0094 /*!
0095  * Get the value at the given grid point.
0096  */
0097 template<class T>
0098 CELER_FUNCTION auto
0099 NonuniformGrid<T>::operator[](size_type i) const -> value_type
0100 {
0101     CELER_EXPECT(i < offset_.size());
0102     return storage_[offset_[i]];
0103 }
0104
0105 //---------------------------------------------------------------------------//
0106 /*!
0107  * Find the value bin such that storage[result] <= value < data[result + 1].
0108  *
0109  * The given value *must* be in range, because out-of-bounds values usually
0110  * require different treatment (e.g. clipping to the boundary values rather
0111  * than interpolating). It's easier to test the exceptional cases (final grid
0112  * point) outside of the grid view.
0113  */
0114 template<class T>
0115 CELER_FUNCTION size_type NonuniformGrid<T>::find(value_type value) const
0116 {
0117     CELER_EXPECT(value >= this->front() && value < this->back());
0118
0119     using ItemIdT = ItemId<T>;
0120     auto iter = celeritas::lower_bound(
0121         offset_.begin(),
0122         offset_.end(),
0123         value,
0124         [&v = storage_](ItemIdT i, T value) { return v[i] < value; });
0125     CELER_ASSERT(iter != offset_.end());
0126
0127     if (value != storage_[*iter])
0128     {
0129         // Exactly on end grid point, or not on a grid point at all: move to
0130         // previous bin
0131         --iter;
0132     }
0133
0134     return iter - offset_.begin();
0135 }
0136
0137 //---------------------------------------------------------------------------//
0138 /*!
0139  * Construct a span referring to the grid points.
0140  */
0141 template<class T>
0142 CELER_FUNCTION auto NonuniformGrid<T>::values() const -> SpanConstT
0143 {
0144     return storage_[offset_];
0145 }
0146
0147 //---------------------------------------------------------------------------//
0148 }  // namespace celeritas