.. _program_listing_file_src_interpolation_Mesh3DInterpolator.hpp:

Program Listing for File Mesh3DInterpolator.hpp
===============================================

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

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

.. code-block:: cpp

   
   #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