Program Listing for File Mesh2DInterpolator.hpp

↰ Return to documentation for file (src/interpolation/Mesh2DInterpolator.hpp)

#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