Program Listing for File Mesh3DInterpolator.hpp

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

#ifndef MESH3DINTERPOLATOR_HPP
#define MESH3DINTERPOLATOR_HPP

#include "Interpolator1D.hpp"
#include "Mesh2DGrid.hpp"
#include "Mesh2DInterpolator.hpp"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <gsl/gsl_interp.h>
#include <mdspan> // NOLINT(misc-include-cleaner)
#include <memory>
#include <span>
#include <stdexcept>
#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 = 1>
    class Mesh3DInterpolator
    {
    public:

        // Shorten array types
        using Array4D = std::mdspan<double, std::dextents<size_t, 4>>; // NOLINT(misc-include-cleaner)
        using Array3D = std::mdspan<double, std::dextents<size_t, 3>>;
        using Array1D = std::mdspan<double, std::dextents<size_t, 1>>;

        // Constructors and destructor
        Mesh3DInterpolator(
            const Array3D& x,
            const Array1D& y,
            const Array1D& z,
            const Array4D& f,
            const gsl_interp_type* interpType = gsl_interp_steffen,
            const bool monotonic = false
        ) :
        monotonic_(monotonic),
        interpType_(monotonic ? gsl_interp_linear : interpType),
        nx_(x.extent(0)),
        ny_(x.extent(1)),
        nz_(x.extent(2))
        {
            if (x.extent(1) != y.extent(0) || x.extent(2) != z.extent(0))
            {
                throw std::runtime_error(
                    "Mesh3DInterpolator: y and z sizes must match the "
                    "second and third dimensions of x");
            }
            if (f.extent(0) != x.extent(0) ||
                f.extent(1) != x.extent(1) ||
                f.extent(2) != x.extent(2) ||
                f.extent(3) != NF)
            {
                throw std::runtime_error(
                    "Mesh3DInterpolator: x and f must have the same leading "
                    "three dimensions, and the fourth dimension of f must "
                    "match the number of interpolated quantities");
            }
            if (x.extent(0) < 2 ||
                (x.extent(1) < 2 && x.extent(2) < 2))
            {
                throw std::runtime_error(
                    "Mesh3DInterpolator: mesh size nx x ny x nz must obey "
                    "nx > 1 and either ny > 1 or nz > 1"
                );
            }
            copyMeshData(x, y, z, f);
            computeLengths();
            computeInterpolators();
        }

        Mesh3DInterpolator(
            const Array3D& x,
            const Array1D& y,
            const Array1D& z,
            const Array3D& f,
            const gsl_interp_type* interpType = gsl_interp_steffen,
            const bool monotonic = false
        ) requires (NF == 1) :
        Mesh3DInterpolator(x, y, z,
            std::mdspan(f.data_handle(), f.extent(0), f.extent(1), f.extent(2), 1), // NOLINT(misc-include-cleaner)
            interpType,
            monotonic)
        { }

        virtual ~Mesh3DInterpolator() = default;

        // Disallow copy constructor due to the unique_ptrs we use
        // internally, but allow move constructor
        Mesh3DInterpolator(const Mesh3DInterpolator&) = delete;
        Mesh3DInterpolator(Mesh3DInterpolator&&) = default;
        auto operator=(const Mesh3DInterpolator&) -> Mesh3DInterpolator& = delete;
        auto operator=(Mesh3DInterpolator&&) noexcept -> Mesh3DInterpolator& = default;

        // Observers
        [[nodiscard]] auto size() const { return nx_ * ny_ * nz_; }

        [[nodiscard]] auto nx() const { return nx_; }

        [[nodiscard]] auto ny() const { return ny_; }

        [[nodiscard]] auto nz() const { return nz_; }

        [[nodiscard]] auto x() const
        {
            return std::mdspan<const double, std::dextents<size_t, 3>>(
                xData_.data(), nx_, ny_, nz_);
        }

        [[nodiscard]] auto y() const
        {
            return std::mdspan<const double, std::dextents<size_t, 1>>(
                yData_.data(), ny_);
        }

        [[nodiscard]] auto z() const
        {
            return std::mdspan<const double, std::dextents<size_t, 1>>(
                zData_.data(), nz_);
        }

        [[nodiscard]] auto f() const
        {
            return std::mdspan<const double, std::dextents<size_t, 4>>(
                fData_.data(), nx_, ny_, nz_, NF);
        }

        [[nodiscard]] auto xMin() const { return *std::ranges::min_element(xData_); }

        [[nodiscard]] auto xMax() const { return *std::ranges::max_element(xData_); }

        [[nodiscard]] auto yMin() const { return yData_.front(); }

        [[nodiscard]] auto yMax() const { return yData_.back(); }

        [[nodiscard]] auto zMin() const { return zData_.front(); }

        [[nodiscard]] auto zMax() const { return zData_.back(); }

        [[nodiscard]] auto sy() const
        {
            return std::mdspan<const double, std::dextents<size_t, 3>>(
                syData_.data(), nx_, nz_, ny_);
        }

        [[nodiscard]] auto sz() const
        {
            return std::mdspan<const double, std::dextents<size_t, 3>>(
                szData_.data(), nx_, ny_, nz_);
        }

        [[nodiscard]] auto yInterp() const
        {
            return std::mdspan<const std::unique_ptr<Interpolator1D<NF>>, std::dextents<size_t, 2>>(
                yInterpData_.data(), nx_, nz_);
        }

        [[nodiscard]] auto zInterp() const
        {
            return std::mdspan<const std::unique_ptr<Interpolator1D<NF>>, std::dextents<size_t, 2>>(
                zInterpData_.data(), nx_, ny_);
        }

        [[nodiscard]] auto sliceConstYCopy(double y0) const -> Mesh2DInterpolator<NF>
        {
            return buildSliceY(y0);
        }

        [[nodiscard]] auto sliceConstZCopy(double z0) const -> Mesh2DInterpolator<NF>
        {
            return buildSliceZ(z0);
        }

        [[nodiscard]] auto lazySliceConstZ(double z0) const -> Mesh2DInterpolator<NF>
        {
            assert(z0 >= zMin() && z0 <= zMax());
            if (nz_ > 1 && ny_ == 1)
            {
                throw std::runtime_error(
                    "Mesh3DInterpolator::lazySliceConstZ: cannot slice at "
                    "fixed z because the mesh has only one point in the "
                    "y direction");
            }

            // Bracketing planes, and interpolation weight between them;
            // an exact match with a mesh z value (including the
            // single-plane case) uses that plane alone, as
            // buildSliceZ() does
            size_t k0 = 0;
            bool exact = (nz_ == 1);
            for (size_t k = 0; k < nz_ && !exact; ++k)
            {
                if (zData_[k] == z0) { k0 = k; exact = true; }
            }
            double t = 0.0;
            if (!exact)
            {
                k0 = static_cast<size_t>(
                    std::ranges::upper_bound(zData_, z0) - zData_.begin()) - 1; // NOLINT(misc-include-cleaner) -- <algorithm> is included
                t = (z0 - zData_[k0]) / (zData_[k0 + 1] - zData_[k0]);
            }

            const std::span<const double> xs(xData_);
            const Mesh2DGrid::XView xView = exact ?
                Mesh2DGrid::XView(xs.subspan(k0), std::span<const double>(),
                    ny_ * nz_, nz_, nx_, ny_, 0.0) :
                Mesh2DGrid::XView(xs.subspan(k0), xs.subspan(k0 + 1),
                    ny_ * nz_, nz_, nx_, ny_, t);
            std::vector<double> yCopy(yData_.begin(), yData_.end());
            const auto yView = std::mdspan<double, std::dextents<size_t, 1>>(yCopy.data(), ny_);
            return Mesh2DInterpolator<NF>(xView, yView,
                std::make_shared<const ZSliceSource>(f(), sz(), zInterp(), k0, t, exact),
                interpType_, monotonic_);
        }

    private:

        class ZSliceSource : public LazyValueSource<NF>
        {
        public:
            using FView = std::mdspan<const double, std::dextents<size_t, 4>>;
            using SzView = std::mdspan<const double, std::dextents<size_t, 3>>;
            using ZInterpView = std::mdspan<const std::unique_ptr<Interpolator1D<NF>>,
                std::dextents<size_t, 2>>;
            ZSliceSource(const FView f, const SzView sz, const ZInterpView zInterp,
                const size_t k0, const double t, const bool exact) :
                f_(f), sz_(sz), zInterp_(zInterp), k0_(k0), t_(t), exact_(exact) {}

            [[nodiscard]] auto value(const size_t i, const size_t j) const
                -> std::array<double, NF> override
            {
                std::array<double, NF> result{};
                if (exact_)
                {
                    for (size_t n = 0; n < NF; ++n) { result[n] = f_[i, j, k0_, n]; } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- n is a loop index bounded by compile-time constant NF
                    return result;
                }
                const double szTarget = sz_[i, j, k0_] +
                    (t_ * (sz_[i, j, k0_ + 1] - sz_[i, j, k0_]));
                const auto fInterp = (*(zInterp_[i, j]))(szTarget);
                if constexpr (NF == 1) { result[0] = fInterp; }
                else
                {
                    for (size_t n = 0; n < NF; ++n) { result[n] = fInterp[n]; } // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- n is a loop index bounded by compile-time constant NF
                }
                return result;
            }

        private:
            FView f_;
            SzView sz_;
            ZInterpView zInterp_;
            size_t k0_;
            double t_;
            bool exact_;
        };

        bool monotonic_;
        const gsl_interp_type* interpType_;
        size_t nx_;
        size_t ny_;
        size_t nz_;
        std::vector<double> xData_;
        std::vector<double> yData_;
        std::vector<double> zData_;
        std::vector<double> fData_;
        std::vector<double> syData_;
        std::vector<double> szData_;
        // Interpolators are stored 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 issues with gsl opaque objects
        std::vector<
            std::unique_ptr<Interpolator1D<NF>>
        > yInterpData_;
        std::vector<
            std::unique_ptr<Interpolator1D<NF>>
        > zInterpData_;
        auto copyMeshData(const Array3D& x, const Array1D& y,
            const Array1D& z, const Array4D& f) -> void
        {
            xData_.resize(nx_ * ny_ * nz_);
            auto xView = Array3D(xData_.data(), nx_, ny_, nz_);
            for (size_t i = 0; i < nx_; ++i)
            {
                for (size_t j = 0; j < ny_; ++j)
                {
                    for (size_t k = 0; k < nz_; ++k) { xView[i, j, k] = x[i, j, k]; }
                }
            }

            yData_.resize(ny_);
            auto yView = Array1D(yData_.data(), ny_);
            for (size_t j = 0; j < ny_; ++j) { yView[j] = y[j]; }

            zData_.resize(nz_);
            auto zView = Array1D(zData_.data(), nz_);
            for (size_t k = 0; k < nz_; ++k) { zView[k] = z[k]; }

            fData_.resize(nx_ * ny_ * nz_ * NF);
            auto fView = Array4D(fData_.data(), nx_, ny_, nz_, NF);
            for (size_t i = 0; i < nx_; ++i)
            {
                for (size_t j = 0; j < ny_; ++j)
                {
                    for (size_t k = 0; k < nz_; ++k)
                    {
                        for (size_t n = 0; n < NF; ++n) { fView[i, j, k, n] = f[i, j, k, n]; }
                    }
                }
            }
        }

        auto computeLengths() -> void
        {
            auto xView = Array3D(xData_.data(), nx_, ny_, nz_);
            auto yView = Array1D(yData_.data(), ny_);
            auto zView = Array1D(zData_.data(), nz_);

            // sy: shape (Nx, Nz, Ny)
            if (ny_ > 1)
            {
                syData_.resize(nx_ * nz_ * ny_);
                auto syView = std::mdspan<double, std::dextents<size_t, 3>>(
                    syData_.data(), nx_, nz_, ny_);
                for (size_t i = 0; i < nx_; ++i)
                {
                    for (size_t k = 0; k < nz_; ++k)
                    {
                        syView[i, k, 0] = 0.0;
                        for (size_t j = 1; j < ny_; ++j)
                        {
                            syView[i, k, j] = syView[i, k, j - 1] + std::sqrt(
                                std::pow(xView[i, j, k] - xView[i, j - 1, k], 2) +
                                std::pow(yView[j] - yView[j - 1], 2));
                        }
                    }
                }
            }

            // sz: shape (Nx, Ny, Nz)
            if (nz_ > 1)
            {
                szData_.resize(nx_ * ny_ * nz_);
                auto szView = std::mdspan<double, std::dextents<size_t, 3>>(
                    szData_.data(), nx_, ny_, nz_);
                for (size_t i = 0; i < nx_; ++i)
                {
                    for (size_t j = 0; j < ny_; ++j)
                    {
                        szView[i, j, 0] = 0.0;
                        for (size_t k = 1; k < nz_; ++k)
                        {
                            szView[i, j, k] = szView[i, j, k - 1] + std::sqrt(
                                std::pow(xView[i, j, k] - xView[i, j, k - 1], 2) +
                                std::pow(zView[k] - zView[k - 1], 2));
                        }
                    }
                }
            }
        }

        template <typename SView, typename FAt>
        auto buildColumnInterpolators(size_t n1, size_t n2, size_t nPrimary,
            SView sView, FAt fAt,
            std::vector<std::unique_ptr<Interpolator1D<NF>>>& out) -> void
        {
            out.resize(n1 * n2);
            auto outView = std::mdspan<std::unique_ptr<Interpolator1D<NF>>, std::dextents<size_t, 2>>(
                out.data(), n1, n2);
            for (size_t a = 0; a < n1; ++a)
            {
                for (size_t b = 0; b < n2; ++b)
                {
                    std::vector<double> col(nPrimary);
                    std::array<std::vector<double>, NF> fCol;
                    // NOLINTBEGIN(cppcoreguidelines-pro-bounds-constant-array-index)
                    for (size_t n = 0; n < NF; ++n) { fCol[n].resize(nPrimary); }
                    for (size_t p = 0; p < nPrimary; ++p)
                    {
                        col[p] = sView[a, b, p];
                        for (size_t n = 0; n < NF; ++n) { fCol[n][p] = fAt(a, b, p, n); }
                    }
                    // NOLINTEND(cppcoreguidelines-pro-bounds-constant-array-index)
                    outView[a, b] = std::make_unique<Interpolator1D<NF>>(col, fCol, interpType_);
                }
            }
        }

        auto computeInterpolators() -> void
        {
            auto fView = Array4D(fData_.data(), nx_, ny_, nz_, NF);

            // y-direction interpolators, one per (i, k) pair
            if (ny_ > 1)
            {
                buildColumnInterpolators(nx_, nz_, ny_, sy(),
                    [&fView](size_t i, size_t k, size_t j, size_t n) -> double { return fView[i, j, k, n]; },
                    yInterpData_);
            }

            // z-direction interpolators, one per (i, j) pair
            if (nz_ > 1)
            {
                buildColumnInterpolators(nx_, ny_, nz_, sz(),
                    [&fView](size_t i, size_t j, size_t k, size_t n) -> double { return fView[i, j, k, n]; },
                    zInterpData_);
            }
        }

        [[nodiscard]] auto buildSliceY(double y0) const -> Mesh2DInterpolator<NF>
        {
            assert(y0 >= yMin() && y0 <= yMax());

            // If there is only one y-plane, the mesh is already entirely
            // at y = y0 (the assert above guarantees y0 matches it), so
            // the slice is just this mesh's (x, z, f) data with the
            // degenerate y dimension dropped -- no interpolation needed.
            // Since ny_ == 1 here, the constructor guard guarantees
            // nz_ >= 2, so the resulting Mesh2DInterpolator is valid.
            // Mesh2DInterpolator requires a mutable view even though it
            // only reads from it to build its own internal copies, so
            // casting away constness here is safe.
            if (ny_ == 1)
            {
                auto* xPtr = const_cast<double*>(xData_.data()); // NOLINT(cppcoreguidelines-pro-type-const-cast) -- see comment above: casting away constness here is safe
                auto* zPtr = const_cast<double*>(zData_.data()); // NOLINT(cppcoreguidelines-pro-type-const-cast) -- see comment above: casting away constness here is safe
                auto* fPtr = const_cast<double*>(fData_.data()); // NOLINT(cppcoreguidelines-pro-type-const-cast) -- see comment above: casting away constness here is safe
                auto xView = std::mdspan<double, std::dextents<size_t, 2>>(xPtr, nx_, nz_);
                auto zView = std::mdspan<double, std::dextents<size_t, 1>>(zPtr, nz_);
                auto fView = std::mdspan<double, std::dextents<size_t, 3>>(fPtr, nx_, nz_, NF);
                return Mesh2DInterpolator<NF>(xView, zView, fView, interpType_, monotonic_);
            }

            // If there is only one z-plane, the (x, z) slice this
            // function is meant to build would have a degenerate z axis
            // of size 1, which Mesh2DInterpolator cannot represent
            if (nz_ == 1)
            {
                throw std::runtime_error(
                    "Mesh3DInterpolator::sliceConstYCopy: cannot slice at "
                    "fixed y because the mesh has only one point in the "
                    "z direction");
            }

            auto xv = x();
            auto fv = f();

            // The Mesh2DInterpolator constructor takes mutable mdspan
            // views, so the z coordinate array must be copied
            std::vector<double> zCopy(zData_.begin(), zData_.end());
            auto zSlice = std::mdspan<double, std::dextents<size_t, 1>>(zCopy.data(), nz_);

            std::vector<double> xSliceData(nx_ * nz_);
            std::vector<double> fSliceData(nx_ * nz_ * NF);
            auto xSlice = std::mdspan<double, std::dextents<size_t, 2>>(xSliceData.data(), nx_, nz_);
            auto fSlice = std::mdspan<double, std::dextents<size_t, 3>>(fSliceData.data(), nx_, nz_, NF);

            // Check for an exact match with one of the mesh y values
            for (size_t j = 0; j < ny_; ++j)
            {
                if (yData_[j] == y0)
                {
                    for (size_t i = 0; i < nx_; ++i)
                    {
                        for (size_t k = 0; k < nz_; ++k)
                        {
                            xSlice[i, k] = xv[i, j, k];
                            for (size_t n = 0; n < NF; ++n) { fSlice[i, k, n] = fv[i, j, k, n]; }
                        }
                    }
                    return Mesh2DInterpolator<NF>(xSlice, zSlice, fSlice, interpType_, monotonic_);
                }
            }

            // No exact match; find the bracketing y planes via binary
            // search and interpolate
            const size_t j0 = static_cast<size_t>(
                std::ranges::upper_bound(yData_, y0) - yData_.begin()) - 1; // NOLINT(misc-include-cleaner)
            const double t = (y0 - yData_[j0]) / (yData_[j0 + 1] - yData_[j0]);

            auto syView = sy();
            auto yInterpView = yInterp();
            for (size_t i = 0; i < nx_; ++i)
            {
                for (size_t k = 0; k < nz_; ++k)
                {
                    xSlice[i, k] = xv[i, j0, k] +
                        (t * (xv[i, j0 + 1, k] - xv[i, j0, k]));

                    const double syTarget = syView[i, k, j0] +
                        (t * (syView[i, k, j0 + 1] - syView[i, k, j0]));
                    const auto fInterp = (*(yInterpView[i, k]))(syTarget);
                    storeSliceResult(fSlice, i, k, fInterp);
                }
            }

            return Mesh2DInterpolator<NF>(xSlice, zSlice, fSlice, interpType_, monotonic_);
        }

        [[nodiscard]] auto buildSliceZ(double z0) const -> Mesh2DInterpolator<NF>
        {
            assert(z0 >= zMin() && z0 <= zMax());

            // If there is only one z-plane, the mesh is already entirely
            // at z = z0 (the assert above guarantees z0 matches it), so
            // the slice is just this mesh's (x, y, f) data with the
            // degenerate z dimension dropped -- no interpolation needed.
            // Since nz_ == 1 here, the constructor guard guarantees
            // ny_ >= 2, so the resulting Mesh2DInterpolator is valid.
            // Mesh2DInterpolator requires a mutable view even though it
            // only reads from it to build its own internal copies, so
            // casting away constness here is safe.
            if (nz_ == 1)
            {
                auto* xPtr = const_cast<double*>(xData_.data()); // NOLINT(cppcoreguidelines-pro-type-const-cast) -- see comment above: casting away constness here is safe
                auto* yPtr = const_cast<double*>(yData_.data()); // NOLINT(cppcoreguidelines-pro-type-const-cast) -- see comment above: casting away constness here is safe
                auto* fPtr = const_cast<double*>(fData_.data()); // NOLINT(cppcoreguidelines-pro-type-const-cast) -- see comment above: casting away constness here is safe
                auto xView = std::mdspan<double, std::dextents<size_t, 2>>(xPtr, nx_, ny_);
                auto yView = std::mdspan<double, std::dextents<size_t, 1>>(yPtr, ny_);
                auto fView = std::mdspan<double, std::dextents<size_t, 3>>(fPtr, nx_, ny_, NF);
                return Mesh2DInterpolator<NF>(xView, yView, fView, interpType_, monotonic_);
            }

            // If there is only one y-plane, the (x, y) slice this
            // function is meant to build would have a degenerate y axis
            // of size 1, which Mesh2DInterpolator cannot represent
            if (ny_ == 1)
            {
                throw std::runtime_error(
                    "Mesh3DInterpolator::sliceConstZCopy: cannot slice at "
                    "fixed z because the mesh has only one point in the "
                    "y direction");
            }

            auto xv = x();
            auto fv = f();

            // The Mesh2DInterpolator constructor takes mutable mdspan
            // views, so the y coordinate array must be copied
            std::vector<double> yCopy(yData_.begin(), yData_.end());
            auto ySlice = std::mdspan<double, std::dextents<size_t, 1>>(yCopy.data(), ny_);

            std::vector<double> xSliceData(nx_ * ny_);
            std::vector<double> fSliceData(nx_ * ny_ * NF);
            auto xSlice = std::mdspan<double, std::dextents<size_t, 2>>(xSliceData.data(), nx_, ny_);
            auto fSlice = std::mdspan<double, std::dextents<size_t, 3>>(fSliceData.data(), nx_, ny_, NF);

            // Check for an exact match with one of the mesh z values
            for (size_t k = 0; k < nz_; ++k)
            {
                if (zData_[k] == z0)
                {
                    for (size_t i = 0; i < nx_; ++i)
                    {
                        for (size_t j = 0; j < ny_; ++j)
                        {
                            xSlice[i, j] = xv[i, j, k];
                            for (size_t n = 0; n < NF; ++n) { fSlice[i, j, n] = fv[i, j, k, n]; }
                        }
                    }
                    return Mesh2DInterpolator<NF>(xSlice, ySlice, fSlice, interpType_, monotonic_);
                }
            }

            // No exact match; find the bracketing z planes via binary
            // search and interpolate
            const size_t k0 = static_cast<size_t>(
                std::ranges::upper_bound(zData_, z0) - zData_.begin()) - 1;
            const double t = (z0 - zData_[k0]) / (zData_[k0 + 1] - zData_[k0]);

            auto szView = sz();
            auto zInterpView = zInterp();
            for (size_t i = 0; i < nx_; ++i)
            {
                for (size_t j = 0; j < ny_; ++j)
                {
                    xSlice[i, j] = xv[i, j, k0] +
                        (t * (xv[i, j, k0 + 1] - xv[i, j, k0]));

                    const double szTarget = szView[i, j, k0] +
                        (t * (szView[i, j, k0 + 1] - szView[i, j, k0]));
                    const auto fInterp = (*(zInterpView[i, j]))(szTarget);
                    storeSliceResult(fSlice, i, j, fInterp);
                }
            }

            return Mesh2DInterpolator<NF>(xSlice, ySlice, fSlice, interpType_, monotonic_);
        }

        static auto storeSliceResult(
            std::mdspan<double, std::dextents<size_t, 3>> fSlice,
            size_t d0, size_t d1, const auto& fInterp) -> void
        {
            if constexpr (NF == 1) { fSlice[d0, d1, 0] = fInterp; }
            else
            {
                // NOLINTBEGIN(cppcoreguidelines-pro-bounds-constant-array-index)
                for (size_t n = 0; n < NF; ++n) { fSlice[d0, d1, n] = fInterp[n]; }
                // NOLINTEND(cppcoreguidelines-pro-bounds-constant-array-index)
            }
        }
    };

} // namespace interp

// NOLINTEND(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)

#endif // MESH3DINTERPOLATOR_HPP