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. The grid is closed on the lower bound and open on the upper,
0023  * and the same treatment is used for coincident grid points: the higher of
0024  * the points will be chosen.
0025  */
0026 template<class T>
0027 class NonuniformGrid
0028 {
0029   public:
0030     //!@{
0031     //! \name Type aliases
0032     using value_type = T;
0033     using Storage
0034         = Collection<value_type, Ownership::const_reference, MemSpace::native>;
0035     using ItemRangeT = ItemRange<value_type>;
0036     using SpanConstT = typename Storage::SpanConstT;
0037     //!@}
0038
0039   public:
0040     // Construct with storage
0041     inline CELER_FUNCTION
0042     NonuniformGrid(ItemRangeT const& values, Storage const& storage);
0043
0044     //! Number of grid points
0045     CELER_FORCEINLINE_FUNCTION size_type size() const
0046     {
0047         return offset_.size();
0048     }
0049
0050     //! Minimum/first value
0051     CELER_FORCEINLINE_FUNCTION value_type front() const
0052     {
0053         return storage_[*offset_.begin()];
0054     }
0055
0056     //! Maximum/last value
0057     CELER_FORCEINLINE_FUNCTION value_type back() const
0058     {
0059         return storage_[*(offset_.end() - 1)];
0060     }
0061
0062     // Calculate the value at the given grid point
0063     inline CELER_FUNCTION value_type operator[](size_type i) const;
0064
0065     // Find the index of the given value (*must* be in bounds)
0066     inline CELER_FUNCTION size_type find(value_type value) const;
0067
0068     //! Low-level access to offsets for downstream utilities
0069     CELER_FORCEINLINE_FUNCTION ItemRangeT offset() const { return offset_; }
0070
0071     // Construct a span referring to the grid points
0072     inline CELER_FUNCTION SpanConstT values() const;
0073
0074   private:
0075     Storage const& storage_;
0076     ItemRangeT offset_;
0077 };
0078
0079 //---------------------------------------------------------------------------//
0080 // INLINE DEFINITIONS
0081 //---------------------------------------------------------------------------//
0082 /*!
0083  * Construct with a range indexing into backend storage.
0084  */
0085 template<class T>
0086 CELER_FUNCTION NonuniformGrid<T>::NonuniformGrid(ItemRangeT const& values,
0087                                                  Storage const& storage)
0088     : storage_{storage}, offset_{values}
0089 {
0090     CELER_EXPECT(offset_.size() >= 2);
0091     CELER_EXPECT(*offset_.end() <= storage.size());
0092     CELER_EXPECT(this->front() <= this->back());  // Approximation for "sorted"
0093 }
0094
0095 //---------------------------------------------------------------------------//
0096 /*!
0097  * Get the value at the given grid point.
0098  */
0099 template<class T>
0100 CELER_FUNCTION auto NonuniformGrid<T>::operator[](size_type i) const
0101     -> value_type
0102 {
0103     CELER_EXPECT(i < offset_.size());
0104     return storage_[offset_[i]];
0105 }
0106
0107 //---------------------------------------------------------------------------//
0108 /*!
0109  * Find the value bin such that data[result] <= value < data[result + 1].
0110  *
0111  * The given value *must* be in range, because out-of-bounds values usually
0112  * require different treatment (e.g. clipping to the boundary values rather
0113  * than interpolating). It's easier to test the exceptional cases (final grid
0114  * point) outside of the grid view.
0115  */
0116 template<class T>
0117 CELER_FUNCTION size_type NonuniformGrid<T>::find(value_type value) const
0118 {
0119     CELER_EXPECT(value >= this->front() && value < this->back());
0120
0121     auto iter = celeritas::upper_bound(
0122         offset_.begin() + 1,
0123         offset_.end() - 1,
0124         value,
0125         [&v = storage_](T value, ItemId<T> i) { return value < v[i]; });
0126
0127     if (value < storage_[*iter])
0128     {
0129         // The start point belongs to the previous bin
0130         --iter;
0131     }
0132
0133     return iter - offset_.begin();
0134 }
0135
0136 //---------------------------------------------------------------------------//
0137 /*!
0138  * Construct a span referring to the grid points.
0139  */
0140 template<class T>
0141 CELER_FUNCTION auto NonuniformGrid<T>::values() const -> SpanConstT
0142 {
0143     return storage_[offset_];
0144 }
0145
0146 //---------------------------------------------------------------------------//
0147 }  // namespace celeritas