Program Listing for File Mesh2DGrid.hpp
↰ Return to documentation for file (src/interpolation/Mesh2DGrid.hpp)
#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