.. _program_listing_file_src_interpolation_Mesh2DInterpolator.hpp:

Program Listing for File Mesh2DInterpolator.hpp
===============================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_interpolation_Mesh2DInterpolator.hpp>` (``src/interpolation/Mesh2DInterpolator.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #ifndef MESH2DINTERPOLATOR_HPP
   #define MESH2DINTERPOLATOR_HPP
   
   #include "Interpolator1D.hpp"
   #include "Mesh2DGrid.hpp"
   #include <algorithm>
   #include <array>
   #include <cmath>
   #include <cstddef>
   #include <gsl/gsl_interp.h>
   #include <mdspan> // NOLINT(misc-include-cleaner)
   #include <memory>
   #include <stdexcept>
   #include <utility>
   #include <vector>
   
   // Disable linting for array bounds checking in this
   // file, since the overhead associated with enforcing
   // such checks severely interferes with performance
   // NOLINTBEGIN(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
   
   
   namespace interp
   {
   
       template <size_t NF>
       class LazyValueSource
       {
       public:
           LazyValueSource() = default;
           virtual ~LazyValueSource() = default;
           LazyValueSource(const LazyValueSource&) = default;
           LazyValueSource(LazyValueSource&&) = default;
           auto operator=(const LazyValueSource&) -> LazyValueSource& = default;
           auto operator=(LazyValueSource&&) -> LazyValueSource& = default;
   
           [[nodiscard]] virtual auto value(size_t i, size_t j) const -> std::array<double, NF> = 0;
       };
   
       template <size_t NF = 1>
       class Mesh2DInterpolator
       {
       public:
   
           // Shorten array types
           using Array3D = std::mdspan<double, std::dextents<size_t, 3>>; // NOLINT(misc-include-cleaner)
           using Array2D = std::mdspan<double, std::dextents<size_t, 2>>;
           using Array1D = std::mdspan<double, std::dextents<size_t, 1>>;
   
           // Constructors and destructor
           Mesh2DInterpolator(
               const Array2D& x,
               const Array1D& y,
               const Array3D& f,
               const gsl_interp_type* interpType = gsl_interp_steffen,
               const bool monotonic = false
           ) :
           interpType_(monotonic ? gsl_interp_linear : interpType),
           monotonic_(monotonic),
           mesh_(x, y, true)
           {
               // Safety check
               if (f.extent(0) != x.extent(0) || 
                   f.extent(1) != x.extent(1) ||
                   f.extent(2) != NF)
               {
                   throw std::runtime_error(
                       "Mesh2DInterpolator: x and f must have same leading "
                       "two dimensions, and the third dimension of f must "
                       "match the number of interpolated quantities");
               }
   
               // Build rib interpolators
               for (size_t j = 0; j < ny(); ++j)
               {
                   std::vector<double> xRib(nx());
                   std::array<std::vector<double>, NF> fRib;
                   for (size_t k = 0; k < NF; k++) { fRib[k].resize(nx()); } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k is a loop index bounded by compile-time constant NF
                   for (size_t i = 0; i < nx(); ++i)
                   {
                       xRib[i] = x[i,j];
                       for (size_t k = 0; k < NF; ++k) { fRib[k][i] = f[i,j,k]; } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k is a loop index bounded by compile-time constant NF
                   }
                   auto ri = std::make_unique<Interpolator1D<NF>>(xRib, fRib, interpType_);
                   ribInterp_.push_back(std::move(ri));
               }
   
               // Build the spine interpolators
               for (size_t i = 0; i < nx(); ++i)
               {
                   std::vector<double> xSpine(ny());
                   std::array<std::vector<double>, NF> fSpine;
                   for (size_t k = 0; k < NF; k++) { fSpine[k].resize(ny()); } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k is a loop index bounded by compile-time constant NF
                   for (size_t j = 0; j < ny(); ++j)
                   {
                       if (!monotonic_) { xSpine[j] = s()[i,j]; }
                       else { xSpine[j] = y[j]; }
                       for (size_t k = 0; k < NF; ++k) { fSpine[k][j] = f[i,j,k]; } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k is a loop index bounded by compile-time constant NF
                   }
                   auto si = std::make_unique<Interpolator1D<NF>>(xSpine, fSpine, interpType_);
                   spineInterp_.push_back(std::move(si));
               }
           }
   
           Mesh2DInterpolator(
               const Array2D& x,
               const Array1D& y,
               const Array2D& f,
               const gsl_interp_type* interpType = gsl_interp_steffen,
               const bool monotonic = false
           ) requires (NF == 1) :
           Mesh2DInterpolator(x, y,
               std::mdspan(f.data_handle(), f.extent(0), f.extent(1), 1), // NOLINT(misc-include-cleaner)
               interpType,
               monotonic)
           { }
   
           Mesh2DInterpolator(
               const Mesh2DGrid::XView& x,
               const Array1D& y,
               std::shared_ptr<const LazyValueSource<NF>> source,
               const gsl_interp_type* interpType = gsl_interp_steffen,
               const bool monotonic = false
           ) :
           interpType_(monotonic ? gsl_interp_linear : interpType),
           monotonic_(monotonic),
           mesh_(x, y),
           lazySource_(std::move(source))
           {
               if (interpType_ != gsl_interp_linear &&
                   interpType_ != gsl_interp_steffen &&
                   interpType_ != gsl_interp_akima)
               {
                   throw std::runtime_error(
                       "Mesh2DInterpolator: lazily-evaluated meshes support "
                       "only local interpolation types (linear, steffen, "
                       "akima)");
               }
           }
   
           virtual ~Mesh2DInterpolator() = default;
   
           // Disallow copy constructor due to the unique_ptrs we use
           // internally, but allow move constructor
           Mesh2DInterpolator(const Mesh2DInterpolator&) = delete;
           Mesh2DInterpolator(Mesh2DInterpolator&&) = default;
           auto operator=(const Mesh2DInterpolator&) -> Mesh2DInterpolator& = delete;
           auto operator=(Mesh2DInterpolator&&) noexcept -> Mesh2DInterpolator& = default;
   
           // Observers
           [[nodiscard]] auto size() const { return mesh_.size(); }
   
           [[nodiscard]] auto nx() const { return mesh_.nx(); }
   
           [[nodiscard]] auto ny() const { return mesh_.ny(); }
   
           [[nodiscard]] auto convex() const { return mesh_.convex(); }
   
           [[nodiscard]] auto x() const -> const auto& { return mesh_.xData(); }
   
           [[nodiscard]] auto y() const -> const auto& { return mesh_.yData(); }
   
           [[nodiscard]] auto s() const -> const auto& { return mesh_.sData(); }
   
           [[nodiscard]] auto xMin() const { return mesh_.xMin(); }
   
           [[nodiscard]] auto xMin(const double y) const { return mesh_.xMin(y); }
   
           [[nodiscard]] auto xMax() const { return mesh_.xMax(); }
   
           [[nodiscard]] auto xMax(const double y) const { return mesh_.xMax(y); }
   
          [[nodiscard]] auto yMin() const { return mesh_.yMin(); }
   
           [[nodiscard]] auto yMax() const { return mesh_.yMax(); }
   
           [[nodiscard]] auto yLim(double x) const { return mesh_.yLim(x); }
   
           [[nodiscard]] auto yEdgeSlope(double x, bool leftEdge) const
           { return mesh_.yEdgeSlope(x, leftEdge); }
   
           [[nodiscard]] auto contains(const double x, const double y) const
           { return mesh_.contains(x, y); }
   
   
           // Interpolators
           [[nodiscard]] auto interpConstX(double x) const -> 
           std::vector<std::unique_ptr<Interpolator1D<NF>>>
           {
               // Grab list of intersection points
               const auto intersect = mesh_.xIntersect(x);
   
               // Output holders
               std::vector<std::unique_ptr<Interpolator1D<NF>>> result;
               std::vector<double> y;
               std::array<std::vector<double>, NF> f;
   
               // Build an interpolator from the accumulated segment and
               // push it onto result, unless the segment has too few
               // points to interpolate at all (e.g. a single isolated
               // tangent point, where the line of constant x touches the
               // mesh at only one mass); such a segment carries no
               // interpolable information, so it is silently dropped
               // rather than attempted, exactly as if that portion of
               // the mesh had been missed entirely
               const auto flushSegment = [&]() -> void
               {
                   if (y.size() < 2) { return; }
                   auto newInterp =
                       std::make_unique<Interpolator1D<NF>>(y, f, interpType_);
                   result.push_back(std::move(newInterp));
               };
   
               // Loop over intersection points
               for (const auto& pt : intersect)
               {
   
                   // Add y to accumulator
                   y.push_back(pt.y);
   
                   // Get point at which to evaluate spline
                   const double s =
                       (pt.t == Mesh2DGrid::IntersectionType::spine &&
                       monotonic_) ? pt.y : pt.xs;
   
                   // Evaluate spline and save
                   const auto fInterp = pt.t == Mesh2DGrid::IntersectionType::rib ?
                       evalRib(pt.idx, s) : evalSpine(pt.idx, s);
                   if constexpr (NF == 1) { f[0].push_back(fInterp); }
                   else
                   {
                       for (size_t k = 0; k < NF; k++) { f[k].push_back(fInterp[k]); } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k is a loop index bounded by compile-time constant NF
                   }
   
                   // If this point is a mesh exit, construct new
                   // interpolator here and push onto output list, then
                   // empty the accumulators
                   if (pt.meshExit)
                   {
                       flushSegment();
                       y.clear();
                       for (auto &fi : f) { fi.clear(); }
                   }
               }
   
               // If any points remain in accumulator, create an interpolator from them
               if (!y.empty())
               {
                   flushSegment();
               }
   
               // Return result
               return result;
           }
           
           [[nodiscard]] auto interpConstY(double y) const -> 
           std::unique_ptr<Interpolator1D<NF>>
           {
               // Grab list of intersection points
               auto intersect = mesh_.yIntersect(y);
   
               // Catch case where mesh is missed entirely, so we just return
               // a null pointer
               if (intersect.empty()) { return nullptr; }
   
               // Accumulators
               std::vector<double> x(intersect.size());
               std::array<std::vector<double>, NF> f;
               for (auto& fi : f) { fi.resize(x.size()); }
   
               // Loop over intersection points, evaluating function at each
               for (size_t i = 0; i < x.size(); i++)
               {
                   const auto& pt = intersect[i];
                   x[i] = pt.x;
                   auto fInterp = evalSpine(pt.idx, pt.s);
                   if constexpr (NF == 1) { f[0][i] = fInterp; }
                   else
                   { 
                       for (size_t k = 0; k < NF; k++) { f[k][i] = fInterp[k]; } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k is a loop index bounded by compile-time constant NF
                   }
               }
   
               // Build final interpolation object to return
               return std::make_unique<Interpolator1D<NF>>(x, f, interpType_);
           }
   
           [[nodiscard]] auto operator()(double x, double y, bool linear = false) const
           {
               // Number of points needed to support the requested
               // interpolation type
               const size_t n = linear ?
                   gsl_interp_type_min_size(gsl_interp_linear) :
                   gsl_interp_type_min_size(interpType_);
   
               // Grab list of intersection points
               const auto intersect = mesh_.xIntersectN(x, y, n);
   
               // A single point means (x, y) is an isolated tangent point
               // of the mesh; there is nothing to interpolate, so just
               // evaluate the appropriate rib or spine interpolator there
               // directly and return
               if (intersect.size() == 1)
               {
                   const auto& pt = intersect.front();
                   const double s =
                       (pt.t == Mesh2DGrid::IntersectionType::spine &&
                       monotonic_) ? pt.y : pt.xs;
                   return pt.t == Mesh2DGrid::IntersectionType::rib ?
                       evalRib(pt.idx, s) : evalSpine(pt.idx, s);
               }
   
               // Accumulators
               std::vector<double> yPts;
               std::array<std::vector<double>, NF> f;
   
               // Loop over intersection points
               for (const auto& pt : intersect)
               {
   
                   // Add y to accumulator
                   yPts.push_back(pt.y);
   
                   // Get point at which to evaluate spline
                   const double s =
                       (pt.t == Mesh2DGrid::IntersectionType::spine &&
                       monotonic_) ? pt.y : pt.xs;
   
                   // Evaluate spline and save
                   const auto fInterp = pt.t == Mesh2DGrid::IntersectionType::rib ?
                       evalRib(pt.idx, s) : evalSpine(pt.idx, s);
                   if constexpr (NF == 1) { f[0].push_back(fInterp); }
                   else
                   {
                       for (size_t k = 0; k < NF; k++) { f[k].push_back(fInterp[k]); } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k is a loop index bounded by compile-time constant NF
                   }
               }
   
               // If fewer points were found than the requested type
               // needs, fall back to linear interpolation
               const gsl_interp_type* useType = interpType_;
               if (linear || intersect.size() < gsl_interp_type_min_size(interpType_))
               {
                   useType = gsl_interp_linear;
               }
   
               // Build interpolator and evaluate at y
               const Interpolator1D<NF> pointInterp(yPts, f, useType);
               return pointInterp(y);
           }
   
   
   
       private:
   
           static constexpr size_t stencilHalfWidth = 3;
   
           [[nodiscard]] auto evalLocal(const auto& coord, const auto& value,
               const size_t n, const double q) const
           {
               size_t lo = 0;
               size_t hi = n - 1;
               while (hi - lo > 1)
               {
                   const size_t mid = (lo + hi) / 2;
                   if (coord(mid) <= q) { lo = mid; }
                   else { hi = mid; }
               }
               size_t first = lo;
               for (size_t nDistinct = 1; first > 0 && nDistinct < stencilHalfWidth; --first)
               {
                   if (coord(first - 1) != coord(first)) { ++nDistinct; }
               }
               while (first > 0 && coord(first - 1) == coord(first)) { --first; }
               size_t last = lo;
               for (size_t nDistinct = 0; last < n - 1 && nDistinct < stencilHalfWidth; ++last)
               {
                   if (coord(last + 1) != coord(last)) { ++nDistinct; }
               }
               std::vector<double> c(last - first + 1);
               std::array<std::vector<double>, NF> f;
               for (auto& fk : f) { fk.resize(c.size()); }
               for (size_t k = 0; k < c.size(); ++k)
               {
                   c[k] = coord(first + k);
                   const auto v = value(first + k);
                   for (size_t n2 = 0; n2 < NF; ++n2) { f[n2][k] = v[n2]; } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- n2 is a loop index bounded by compile-time constant NF
               }
               const Interpolator1D<NF> local(c, f, interpType_);
               return local(q);
           }
   
           [[nodiscard]] auto evalRib(const size_t j, const double x) const
           {
               if (!lazySource_) { return (*(ribInterp_[j]))(x); }
               const auto& xv = mesh_.xData();
               return evalLocal(
                   [&xv, j](const size_t i) -> double { return xv[i, j]; },
                   [this, j](const size_t i) { return lazySource_->value(i, j); },
                   mesh_.nx(), x);
           }
   
           [[nodiscard]] auto evalSpine(const size_t i, const double s) const
           {
               if (!lazySource_) { return (*(spineInterp_[i]))(s); }
               const size_t nyLoc = mesh_.ny();
               std::vector<double> coord(nyLoc);
               if (monotonic_)
               {
                   for (size_t j = 0; j < nyLoc; ++j) { coord[j] = mesh_.yData()[j]; }
               }
               else
               {
                   const auto& xv = mesh_.xData();
                   const auto& yv = mesh_.yData();
                   coord[0] = 0.0;
                   for (size_t j = 1; j < nyLoc; ++j)
                   {
                       coord[j] = coord[j-1] + std::sqrt(
                           std::pow(xv[i,j] - xv[i,j-1], 2) +
                           std::pow(yv[j] - yv[j-1], 2));
                   }
               }
               return evalLocal(
                   [&coord](const size_t j) -> double { return coord[j]; },
                   [this, i](const size_t j) { return lazySource_->value(i, j); },
                   nyLoc, s);
           }
   
           // Control parameters
           const gsl_interp_type *interpType_; 
           bool monotonic_;                    
           // Internal storage; note that we handle the rib and spine
           // interpolators as vectors of unique_ptr rather than vectors
           // of Interpolator1D because we need to be sure to invoke the
           // destructor of Interpolator1D objects once only to avoid memory
           // issue with gsl opaque objects
           Mesh2DGrid mesh_;            
           std::vector<
               std::unique_ptr<Interpolator1D<NF>>
           > ribInterp_;                 
           std::vector<
               std::unique_ptr<Interpolator1D<NF>>
           > spineInterp_;               
           std::shared_ptr<const LazyValueSource<NF>>
               lazySource_;              
       };
   
   } // namespace interp
   
   // NOLINTEND(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
   
   #endif // MESH2DINTERPOLATOR_HPP