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 NonuniformGrid(ItemRangeT const& values,
0042                                          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     // Construct a span referring to the grid points
0069     inline CELER_FUNCTION SpanConstT values() const;
0070
0071     //! Low-level access to offsets for downstream utilities
0072     CELER_FORCEINLINE_FUNCTION ItemRangeT offset() const { return offset_; }
0073
0074     //! Low-level access to storage for downstream utilities
0075     CELER_FORCEINLINE_FUNCTION Storage const& storage() const
0076     {
0077         return storage_;
0078     }
0079
0080   private:
0081     Storage const storage_;
0082     ItemRangeT offset_;
0083 };
0084
0085 //---------------------------------------------------------------------------//
0086 // INLINE DEFINITIONS
0087 //---------------------------------------------------------------------------//
0088 /*!
0089  * Construct with a range indexing into backend storage.
0090  */
0091 template<class T>
0092 CELER_FUNCTION NonuniformGrid<T>::NonuniformGrid(ItemRangeT const& values,
0093                                                  Storage const& storage)
0094     : storage_{storage}, offset_{values}
0095 {
0096     CELER_EXPECT(offset_.size() >= 2);
0097     CELER_EXPECT(*offset_.end() <= storage.size());
0098     CELER_EXPECT(this->front() <= this->back());  // Approximation for "sorted"
0099 }
0100
0101 //---------------------------------------------------------------------------//
0102 /*!
0103  * Get the value at the given grid point.
0104  */
0105 template<class T>
0106 CELER_FUNCTION auto NonuniformGrid<T>::operator[](size_type i) const
0107     -> value_type
0108 {
0109     CELER_EXPECT(i < offset_.size());
0110     return storage_[offset_[i]];
0111 }
0112
0113 //---------------------------------------------------------------------------//
0114 /*!
0115  * Find the value bin such that data[result] <= value < data[result + 1].
0116  *
0117  * The given value *must* be in range, because out-of-bounds values usually
0118  * require different treatment (e.g. clipping to the boundary values rather
0119  * than interpolating). It's easier to test the exceptional cases (final grid
0120  * point) outside of the grid view.
0121  */
0122 template<class T>
0123 CELER_FUNCTION size_type NonuniformGrid<T>::find(value_type value) const
0124 {
0125     CELER_EXPECT(value >= this->front() && value < this->back());
0126
0127     auto iter
0128         = celeritas::upper_bound(offset_.begin() + 1,
0129                                  offset_.end() - 1,
0130                                  value,
0131                                  [&v = storage_](T lhs, ItemId<T> rhs_id) {
0132                                      return lhs < v[rhs_id].get();
0133                                  });
0134
0135     if (value < storage_[*iter].get())
0136     {
0137         // The start point belongs to the previous bin
0138         --iter;
0139     }
0140
0141     return iter - offset_.begin();
0142 }
0143
0144 //---------------------------------------------------------------------------//
0145 /*!
0146  * Construct a span referring to the grid points.
0147  */
0148 template<class T>
0149 CELER_FUNCTION auto NonuniformGrid<T>::values() const -> SpanConstT
0150 {
0151     return storage_[offset_];
0152 }
0153
0154 //---------------------------------------------------------------------------//
0155 }  // namespace celeritas