.. _program_listing_file_src_interpolation_Mesh2DGrid.hpp:

Program Listing for File Mesh2DGrid.hpp
=======================================

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

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

.. code-block:: cpp

   
   #ifndef MESH2DGRID_HPP
   #define MESH2DGRID_HPP
   
   #include "../utils/ThreadVec.hpp"
   #include <cassert>
   #include <cmath>
   #include <cstddef>
   #include <cstdint>
   #include <cstdlib>
   #include <limits>
   #include <mdspan> // NOLINT(misc-include-cleaner)
   #include <span>
   #include <tuple>
   #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
   {
   
       class Mesh2DGrid
       {
       public:
   
           // Shorten array types
           using Array2D = std::mdspan<double, std::dextents<size_t, 2>>; // NOLINT(misc-include-cleaner)
           using Array1D = std::mdspan<double, std::dextents<size_t, 1>>; // NOLINT(misc-include-cleaner) -- <mdspan> is included; see Array2D above
   
           // Shorthand for numeric limits
           static constexpr double bigNum = std::numeric_limits<double>::max(); 
           class XView
           {
           public:
               XView() = default;
   
               XView(std::span<const double> data, const size_t nx, const size_t ny) :
                   p0_(data), strideI_(ny), strideJ_(1), nx_(nx), ny_(ny) {}
   
               XView(std::span<const double> p0, std::span<const double> p1,
                   const size_t strideI, const size_t strideJ,
                   const size_t nx, const size_t ny, const double t) :
                   p0_(p0), p1_(p1), strideI_(strideI), strideJ_(strideJ),
                   nx_(nx), ny_(ny), t_(t) {}
   
               [[nodiscard]] auto operator[](const size_t i, const size_t j) const -> double
               {
                   const size_t off = (i * strideI_) + (j * strideJ_);
                   const double v0 = p0_[off];
                   if (p1_.empty()) { return v0; }
                   return v0 + (t_ * (p1_[off] - v0));
               }
   
               [[nodiscard]] auto extent(const size_t r) const -> size_t { return r == 0 ? nx_ : ny_; }
   
               [[nodiscard]] auto size() const -> size_t { return nx_ * ny_; }
   
               [[nodiscard]] auto lazy() const -> bool { return !p1_.empty(); }
   
           private:
               std::span<const double> p0_; 
               std::span<const double> p1_; 
               size_t strideI_ = 0;         
               size_t strideJ_ = 0;         
               size_t nx_ = 0;              
               size_t ny_ = 0;              
               double t_ = 0.0;             
           };
   
           class MView
           {
           public:
               MView() = default;
   
               MView(std::span<const double> data, const size_t ny) :
                   data_(data), ny_(ny) {}
   
               MView(const XView& x, std::span<const double> y) :
                   x_(x), y_(y), ny_(y.size()), lazy_(true) {}
   
               [[nodiscard]] auto operator[](const size_t i, const size_t j) const -> double
               {
                   if (!lazy_) { return data_[(i * (ny_ - 1)) + j]; }
                   const double x0 = x_[i,j];
                   const double x1 = x_[i,j+1];
                   if (x1 != x0) { return (y_[j+1] - y_[j]) / (x1 - x0); }
                   return bigNum;
               }
   
           private:
               std::span<const double> data_; 
               XView x_;                      
               std::span<const double> y_;    
               size_t ny_ = 0;                
               bool lazy_ = false;            
           };
   
           class SView
           {
           public:
               SView() = default;
   
               SView(std::span<const double> data, const size_t ny) :
                   data_(data), ny_(ny) {}
   
               SView(const XView& x, std::span<const double> y) :
                   x_(x), y_(y), ny_(y.size()), lazy_(true) {}
   
               [[nodiscard]] auto operator[](const size_t i, const size_t j) const -> double
               {
                   if (!lazy_) { return data_[(i * ny_) + j]; }
                   double s = 0.0;
                   for (size_t jj = 1; jj <= j; ++jj)
                   {
                       s += std::sqrt(
                           std::pow(x_[i,jj] - x_[i,jj-1], 2) +
                           std::pow(y_[jj] - y_[jj-1], 2));
                   }
                   return s;
               }
   
           private:
               std::span<const double> data_; 
               XView x_;                      
               std::span<const double> y_;    
               size_t ny_ = 0;                
               bool lazy_ = false;            
           };
   
           enum class IntersectionType : std::uint8_t {
               spine, 
               rib,   
               none   
           };
   
           // NOLINTBEGIN(readability-identifier-naming)
           struct xIntersectionDescriptor {
               double y;  
               double xs; 
               IntersectionType t; 
               size_t idx; 
               bool meshExit; 
           };
           // NOLINTEND(readability-identifier-naming)
   
           // NOLINTBEGIN(readability-identifier-naming)
           struct yIntersectionDescriptor {
               double x;  
               double s; 
               size_t idx; 
           };
           // NOLINTEND(readability-identifier-naming)
           
           // Constructors and destructor
           Mesh2DGrid(const Array2D& x,
               const Array1D& y,
               bool copyData = true);
   
           Mesh2DGrid(const XView& x, const Array1D& y);
           virtual ~Mesh2DGrid() = default;
   
           // Copy and move constructors are written out by hand (rather
           // than defaulted) because iSave_ and jSave_ are ThreadVec's,
           // which are themselves non-copyable and non-movable so that
           // each thread's cached search index stays private to that
           // thread; the copy/move constructed object therefore gets its
           // own freshly-initialized search cache rather than one copied
           // or moved from the source, which is correct since the cache
           // is just a performance optimization, not semantic state.
           // Assignment remains disallowed, as it always has been, since
           // x_/y_/m_/s_ are non-owning views into xData_/yData_/mData_/
           // sData_ and a blind member-wise assignment would not repoint
           // them at the assigned-to object's own data.
           Mesh2DGrid(const Mesh2DGrid& other);
           Mesh2DGrid(Mesh2DGrid&& other) noexcept; // NOLINT(bugprone-exception-escape) -- see definition
           auto operator=(const Mesh2DGrid&) -> Mesh2DGrid& = delete;
           auto operator=(Mesh2DGrid&&) -> Mesh2DGrid& = delete;
   
           // Observers
           auto size() const { return x_.size(); }
   
           auto nx() const { return x_.extent(0); }
   
           auto ny() const { return x_.extent(1); }
   
           auto convex() const { return convex_; }
   
           auto xData() const -> const auto& { return x_; }
   
           auto yData() const -> const auto& { return y_; }
   
           auto sData() const -> const auto& { return s_; }
           
           auto xMin() const { return xMin_; }
   
           auto xMin(const double y) const
           {
               assert(y >= yMin_ && y <= yMax_);
               if (y == yMin_) { return x_[0,0]; }
               if (y == yMax_) { return x_[0,ny()-1]; }
               auto j = yIdx(y);
               return x_[0,j] + ((y - y_[j]) / m_[0,j]);
           }
   
           auto xMax() const { return xMax_; }
   
           auto xMax(const double y) const
           {
               assert(y >= yMin_ && y <= yMax_);
               if (y == yMin_) { return x_[nx()-1,0]; }
               if (y == yMax_) { return x_[nx()-1,ny()-1]; }
               auto j = yIdx(y);
               return x_[nx()-1,j] + ((y - y_[j]) / m_[nx()-1,j]);
           }
   
           auto yMin() const { return yMin_; }
   
           auto yMax() const { return yMax_; }
   
           auto yLim(double x)
           const -> std::vector<std::pair<double, double>>;
   
           auto yEdgeSlope(double x, bool leftEdge)
           const -> std::vector<std::pair<double, double>>;
   
           auto contains(const double x, const double y) const
           {
               if (y < yMin_ || y > yMax_) { return false; }
               return (x >= xMin(y) && x <= xMax(y));
           }
   
           // Mesh indexing functions
           auto xIdx(const double x, const size_t j) const -> size_t
           {
               auto& iSaveLoc = iSave_();
               auto& jSaveLoc = jSave_();
               jSaveLoc = j;  // Update cached position
               assert((x >= x_[0,j] && x <= x_[nx()-1,j]));  // Safety check
               if (x < x_[iSaveLoc,j]) { // Check if we need to update cached index
                   iSaveLoc = xSearch(x, 0, j, 0, iSaveLoc); // Below cached position
               } else if (x >= x_[iSaveLoc+1,j]) {
                   // Above cached position, or exactly on its upper edge; in
                   // the latter case iSaveLoc must still advance, since leaving
                   // it unchanged would compute the offset within this cell
                   // as the full cell width instead of zero, which produces
                   // spurious floating-point error in downstream traversals
                   iSaveLoc = xSearch(x, 0, j, iSaveLoc, nx()-1);
               }
               if (iSaveLoc == nx()-1) { --iSaveLoc; } // Handle edge case
               return iSaveLoc;
           }
   
           auto yIdx(const double y) const -> size_t
           {
               auto& jSaveLoc = jSave_();
               assert(y >= yMin_ && y <= yMax_);  // Safety check
               if (y < y_[jSaveLoc]) { // Check if we need to update cached position
                   jSaveLoc = ySearch(y, 0, jSaveLoc);   // Below cached position, search left
               } else if (y >= y_[jSaveLoc + 1]) {
                   // Above cached position, or exactly on its upper edge; in
                   // the latter case jSaveLoc must still advance, since leaving
                   // it unchanged would compute the offset within this cell
                   // as the full cell width instead of zero, which produces
                   // spurious floating-point error in downstream traversals
                   jSaveLoc = ySearch(y, jSaveLoc, ny()-1);
               }
               if (jSaveLoc == ny()-1) { --jSaveLoc; } // Handle edge case
               return jSaveLoc;
           }
   
           auto xyIdx(const double x, const double y) const
           {
               auto& iSaveLoc = iSave_();
               assert(contains(x,y));  // Safety assertion
               auto j = yIdx(y);       // Get y index
               const double dy = y - y_[j];  // Offset in y
               size_t i = 0;  // Output holder
               if (x < x_[iSaveLoc,j] + (dy/m_[iSaveLoc,j])) { // Check if we need to update cached index
                   iSaveLoc = xSearch(x, dy, j, 0, iSaveLoc);  // Too low
                   i = iSaveLoc;
               } else if (x >= x_[iSaveLoc+1,j] + (dy/m_[iSaveLoc+1,j])) {
                   iSaveLoc = xSearch(x, dy, j, iSaveLoc, nx()-1); // Too high
                   i = iSaveLoc;
               } else {
                   i = iSaveLoc; // Cache hit
               }
               if (i == nx()-1) { --iSaveLoc; } // Handle dge case
               return std::pair(i,j);
           }
   
           auto onMesh(double x, double y)
           const -> std::pair<bool, xIntersectionDescriptor>;
   
           // Mesh traversal functions
           auto xIntersect(
               double x,
               double yLo = std::numeric_limits<double>::lowest(),
               double yHi = std::numeric_limits<double>::max()
           ) const -> std::vector<xIntersectionDescriptor>;
   
           auto xIntersectN(
               double x,
               double y,
               size_t n
           ) const -> std::vector<xIntersectionDescriptor>;
   
           auto yIntersect(
               double y,
               double xLo = std::numeric_limits<double>::lowest(),
               double xHi = std::numeric_limits<double>::max()
           ) const -> std::vector<yIntersectionDescriptor>;
           
       private:
   
           // Construction helpers
           static void safetyCheckInputs(const Array2D& x, 
               const Array1D& y);
           void computeSlopesAndLengths();
           void setConvexity();
   
           // Indexing / search helpers
           auto ySearch(const double y,
               const size_t idxLo,
               const size_t idxHi
           ) const -> size_t
           {
               auto jLo = idxLo;
               auto jHi = idxHi;
               while (jHi > jLo + 1) {
                   auto j = (jLo + jHi) / 2;
                   if (y > y_[j]) {
                       jLo = j;
                   } else {
                       jHi = j;
                   }
               }
               if (y == y_[jHi]) { return jHi; }
               return jLo;
           }
           /* @brief Search for index of x position at given y
            * @param x x position for which to search
            * @param dy Distance from y[j]
            * @param j Index of y position
            * @param idxLo Lower index for search
            * @param idxHi Upper index for search
            * @details
            * This routine carries out a binary search for the
            * index at which position y can be found, between
            * idxLo and idxHi. This is a helper function for
            * other functions, so no checking is done to ensure
            * that y is in fact between these indices; the calling
            * routine is responsible for checking that.
            */
           auto xSearch(
               const double x,
               const double dy,
               const size_t j,
               const size_t idxLo,
               const size_t idxHi
           ) const -> size_t
           {
               size_t iLo = idxLo;
               size_t iHi = idxHi;
               while (iHi > iLo + 1) {
                   const auto i = (iLo + iHi) / 2;
                   if (x > x_[i,j] + (dy/m_[i,j])) {
                       iLo = i;
                   } else {
                       iHi = i;
                   }
               }
               if (x == x_[iHi,j] + (dy/m_[iHi,j])) {
                   return iHi; // Edge case
               }
               return iLo;
           }
   
           // Search helpers
           auto yLimTraverse(double x, std::vector<double>& yL) const -> bool;
   
           auto yEdgeSlopeTraverse(double x, std::vector<double>& yL,
               std::vector<std::pair<double, double>>& yEdge, bool leftEdge) const -> bool;
   
           auto xIntersectSeg(
               double x,
               double yMin,
               double yMax,
               bool startInterior,
               bool endInterior
           ) const -> std::vector<xIntersectionDescriptor>;
   
           auto yIntersectBoundaryColumn(
               double y,
               double xLo,
               double xHi
           ) const -> std::vector<yIntersectionDescriptor>;
   
           auto xIntersectBoundaryRow(
               double x,
               double yLo,
               double yHi
           ) const -> std::vector<xIntersectionDescriptor>;
   
           auto xIntersectSegStartInterior(
               double x,
               double& y,
               std::vector<xIntersectionDescriptor>& intList,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> bool;
   
           auto xIntersectSegStartEdge(
               double x,
               double& y,
               double yMin,
               std::vector<xIntersectionDescriptor>& intList,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> bool;
           
           auto findNextIntersect(
               double x,
               double& y,
               double yStop,
               bool endInterior,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> std::pair<bool, xIntersectionDescriptor>;
   
           auto findIntersectDistance(
               double x,
               double y,
               bool searchUp,
               bool lastIntersectLeft,
               bool lastIntersectRight
           ) const -> std::tuple<double, double, double>;
   
           auto handleNextIntersectTopRib(
               double& y,
               double x,
               double yStop,
               bool endInterior,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> std::pair<bool, xIntersectionDescriptor>;
   
           auto handleNextIntersectBottomRib(
               double& y,
               double x,
               double yStop,
               bool endInterior,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> std::pair<bool, xIntersectionDescriptor>;
   
           auto handleNextIntersectLeftSpine(
               double& y,
               double x,
               double yStop,
               bool endInterior,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> std::pair<bool, xIntersectionDescriptor>;
   
           auto handleNextIntersectRightSpine(
               double& y,
               double x,
               double yStop,
               bool endInterior,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> std::pair<bool, xIntersectionDescriptor>;
   
           auto handleTopRibCornerCases(
               double x,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> bool;
    
           auto handleBottomRibCornerCases(
               double x,
               bool& lastIntersectLeft,
               bool& lastIntersectRight
           ) const -> bool;
           
           void rebindViews();
   
           // Input data
           XView x_;                    
           Array1D y_;                  
           std::vector<double> xData_;  
           std::vector<double> yData_;  
           // Descriptors
           double xMin_ = 0.0;           
           double xMax_ = 0.0;           
           double yMin_ = 0.0;           
           double yMax_ = 0.0;           
           bool convex_ = true;          
           // Derived data
           MView m_;                       
           SView s_;                       
           std::vector<double> mData_;     
           std::vector<double> sData_;     
           // Mutables
           //
           // Cached x and y indices for search acceleration. These are
           // ThreadVec's rather than plain size_t's so that each thread
           // gets its own private cached index; without this, concurrent
           // calls from different threads would race to update a single
           // shared cache. Every method that uses these must start by
           // binding a local reference to the element private to the
           // calling thread, e.g. "auto& iSaveLoc = iSave_();", and use
           // that local reference throughout rather than the member
           // directly.
           mutable utils::ThreadVec<size_t> iSave_;    
           mutable utils::ThreadVec<size_t> jSave_;    
       };
   
   } // namespace interp
   
   // NOLINTEND(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
   
   #endif // MESH2DGRID_HPP