corecel/grid/TwodSubgridCalculator.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/TwodSubgridCalculator.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/data/Collection.hh"
0010
0011 #include "FindInterp.hh"
0012 #include "NonuniformGrid.hh"
0013 #include "TwodGridData.hh"
0014
0015 namespace celeritas
0016 {
0017 //---------------------------------------------------------------------------//
0018 /*!
0019  * Do bilinear interpolation on a 2D grid with the x value preselected.
0020  *
0021  * This is usually not called directly but rather given as the return result of
0022  * the TwodGridCalculator.
0023  */
0024 class TwodSubgridCalculator
0025 {
0026   public:
0027     //!@{
0028     //! \name Type aliases
0029     using Values
0030         = Collection<real_type, Ownership::const_reference, MemSpace::native>;
0031     using InterpT = FindInterp<real_type>;
0032     //!@}
0033
0034   public:
0035     // Construct with grid data, backend values, and lower X data.
0036     inline CELER_FUNCTION TwodSubgridCalculator(TwodGridData const& grid,
0037                                                 Values const& storage,
0038                                                 InterpT x_loc);
0039
0040     // Calculate the value at the given y coordinate
0041     inline CELER_FUNCTION real_type operator()(real_type y) const;
0042
0043     // Calculate the value at the given y grid point
0044     inline CELER_FUNCTION real_type operator[](size_type i) const;
0045
0046     //! Index of the preselected lower x value
0047     CELER_FORCEINLINE_FUNCTION size_type x_index() const
0048     {
0049         return x_loc_.index;
0050     }
0051
0052     //! Fraction between the lower and upper x grid values
0053     CELER_FORCEINLINE_FUNCTION real_type x_fraction() const
0054     {
0055         return x_loc_.fraction;
0056     }
0057
0058   private:
0059     TwodGridData const& grids_;
0060     Values const& storage_;
0061     InterpT const x_loc_;
0062
0063     inline CELER_FUNCTION real_type at(size_type x_idx, size_type y_idx) const;
0064 };
0065
0066 //---------------------------------------------------------------------------//
0067 // INLINE DEFINITIONS
0068 //---------------------------------------------------------------------------//
0069 /*!
0070  * Construct with grid data, backend values, and lower X data.
0071  *
0072  * This is typically constructed from a TwodGridCalculator. The interpolated x
0073  * location could be extended to allow a fractional value of 1 to support
0074  * interpolating on the highest value of the x grid.
0075  */
0076 CELER_FUNCTION
0077 TwodSubgridCalculator::TwodSubgridCalculator(TwodGridData const& grids,
0078                                              Values const& storage,
0079                                              InterpT x_loc)
0080     : grids_{grids}, storage_(storage), x_loc_(x_loc)
0081 {
0082     CELER_EXPECT(grids);
0083     CELER_EXPECT(grids.values.back() < storage.size());
0084     CELER_EXPECT(x_loc.index + 1 < grids.x.size());
0085     CELER_EXPECT(x_loc.fraction >= 0 && x_loc_.fraction < 1);
0086 }
0087
0088 //---------------------------------------------------------------------------//
0089 /*!
0090  * Calculate the value at the given y coordinate for preselected x.
0091  *
0092  * This uses *bilinear* interpolation and and therefore exactly represents
0093  * functions that are a linear combination of 1, x, y, and xy.
0094  */
0095 CELER_FUNCTION real_type TwodSubgridCalculator::operator()(real_type y) const
0096 {
0097     NonuniformGrid<real_type> const y_grid{grids_.y, storage_};
0098     CELER_EXPECT(y >= y_grid.front() && y < y_grid.back());
0099
0100     InterpT const y_loc = find_interp(y_grid, y);
0101     auto at_corner = [this, y_loc](size_type xo, size_type yo) -> real_type {
0102         return this->at(x_loc_.index + xo, y_loc.index + yo);
0103     };
0104
0105     return (1 - x_loc_.fraction)
0106                * ((1 - y_loc.fraction) * at_corner(0, 0)
0107                   + (y_loc.fraction) * at_corner(0, 1))
0108            + (x_loc_.fraction)
0109                  * ((1 - y_loc.fraction) * at_corner(1, 0)
0110                     + (y_loc.fraction) * at_corner(1, 1));
0111 }
0112
0113 //---------------------------------------------------------------------------//
0114 /*!
0115  * Calculate the value at the given y grid point for preselected x.
0116  */
0117 CELER_FUNCTION real_type TwodSubgridCalculator::operator[](size_type i) const
0118 {
0119     CELER_EXPECT(i < grids_.y.size());
0120     return (1 - x_loc_.fraction) * this->at(x_loc_.index, i)
0121            + x_loc_.fraction * this->at(x_loc_.index + 1, i);
0122 }
0123
0124 //---------------------------------------------------------------------------//
0125 /*!
0126  * Get the value at the specified x/y coordinate.
0127  *
0128  * NOTE: this must match TwodGridData::index.
0129  */
0130 CELER_FUNCTION real_type TwodSubgridCalculator::at(size_type x_idx,
0131                                                    size_type y_idx) const
0132 {
0133     return storage_[grids_.at(x_idx, y_idx)];
0134 }
0135
0136 //---------------------------------------------------------------------------//
0137 }  // namespace celeritas