Grid.h

Go to the documentation of this file.

#ifndef _Grid_
#define _Grid_
 
#include <algorithm>
#include <iostream>
#include <ostream>
#include <vector>
#include <cassert>
#include <array>
#include <string>
#include <limits>
#include "vapor/VAssert.h"
#include <memory>
#include <vapor/common.h>
 
#ifdef WIN32
    #pragma warning(disable : 4661 4251)    // needed for template class
#endif
 
namespace VAPoR {
 
using CoordType = std::array<double, 3>;
 
using DimsType = std::array<size_t, 3>;
 
//
class VDF_API Grid {
public:
    Grid(const std::vector<size_t> &dims, const std::vector<size_t> &bs, const std::vector<float *> &blks, size_t topology_dimension);
    Grid(const DimsType &dims, const DimsType &bs, const std::vector<float *> &blks, size_t topology_dimension);
 
    Grid();
    virtual ~Grid() = default;
 
    const DimsType &GetDimensions() const { return _dims; }
 
    size_t GetNumDimensions() const { return GetNumDimensions(_dims); }
 
    static size_t GetNumDimensions(DimsType dims);
 
    virtual DimsType GetCoordDimensions(size_t dim) const = 0;
 
    virtual std::string GetType() const = 0;
 
    //
    virtual size_t GetGeometryDim() const = 0;
 
    virtual const DimsType &GetNodeDimensions() const = 0;
    virtual const size_t    GetNumNodeDimensions() const = 0;
    virtual const DimsType &GetCellDimensions() const = 0;
    virtual const size_t    GetNumCellDimensions() const = 0;
 
    //
    const std::vector<size_t> GetDimensionInBlks() const { return (_bdimsDeprecated); }
 
    //
    const std::vector<size_t> &GetBlockSize() const { return (_bsDeprecated); }
 
    const std::vector<float *> &GetBlks() const { return (_blks); };
 
    virtual float GetValueAtIndex(const DimsType &indices) const;
 
    //
    virtual float GetValueAtIndex(const std::vector<size_t> &indices) const
    {
        DimsType a = {0, 0, 0};
        CopyToArr3(indices, a);
        return (GetValueAtIndex(a));
    }
 
    virtual void SetValue(const DimsType &indices, float value);
 
    //
    virtual void SetValue(const size_t indices[3], float value)
    {
        DimsType i3 = {0, 0, 0};
        CopyToArr3(indices, GetNodeDimensions().size(), i3);
        SetValue(i3, value);
    }
 
    //
    virtual float AccessIJK(size_t i, size_t j = 0, size_t k = 0) const;
 
    void SetValueIJK(size_t i, size_t j, size_t k, float v);
    void SetValueIJK(size_t i, size_t j, float v) { SetValueIJK(i, j, 0, v); }
    void SetValueIJK(size_t i, float v) { SetValueIJK(i, 0, 0, v); }
 
    virtual float GetValue(const CoordType &coords) const;
 
    //
    virtual float GetValue(const std::vector<double> &coords) const
    {
        CoordType c3 = {0.0, 0.0, 0.0};
        CopyToArr3(coords, c3);
        return (GetValue(c3));
    };
 
    //
    virtual float GetValue(const double coords[]) const
    {
        CoordType c3 = {0.0, 0.0, 0.0};
        CopyToArr3(coords, GetGeometryDim(), c3);
        return (GetValue(c3));
    }
 
    virtual float GetValue(double x, double y) const
    {
        double coords[] = {x, y, 0.0};
        return (GetValue(coords));
    }
    virtual float GetValue(double x, double y, double z) const
    {
        double coords[] = {x, y, z};
        return (GetValue(coords));
    }
 
    virtual void GetUserExtents(CoordType &minu, CoordType &maxu) const;
 
    //
    virtual void GetUserExtents(double minu[3], double maxu[3]) const
    {
        CoordType minu3 = {0.0, 0.0, 0.0};
        CoordType maxu3 = {0.0, 0.0, 0.0};
        GetUserExtents(minu3, maxu3);
        CopyFromArr3(minu3, minu);
        CopyFromArr3(maxu3, maxu);
    }
 
    //
    virtual void GetUserExtents(std::vector<double> &minu, std::vector<double> &maxu) const
    {
        CoordType minu3 = {0.0, 0.0, 0.0};
        CoordType maxu3 = {0.0, 0.0, 0.0};
        GetUserExtents(minu3, maxu3);
        CopyFromArr3(minu3, minu);
        CopyFromArr3(maxu3, maxu);
 
        // Much of the use of this method assumes that the size of the minu,maxu
        // vectors can be used to determine the number of coordinates. So we
        // maintain this property for now.
        //
        minu.resize(GetGeometryDim());
        maxu.resize(GetGeometryDim());
    }
 
    virtual void GetBoundingBox(const DimsType &min, const DimsType &max, CoordType &minu, CoordType &maxu) const = 0;
 
    virtual void GetBoundingBox(const std::vector<size_t> &min, const std::vector<size_t> &max, std::vector<double> &minu, std::vector<double> &maxu) const
    {
        VAssert(min.size() == max.size());
        DimsType  min3 = {0, 0, 0};
        DimsType  max3 = {0, 0, 0};
        CoordType minu3 = {0.0, 0.0, 0.0};
        CoordType maxu3 = {0.0, 0.0, 0.0};
 
        CopyToArr3(min, min3);
        CopyToArr3(max, max3);
        GetBoundingBox(min3, max3, minu3, maxu3);
        CopyFromArr3(minu3, minu);
        CopyFromArr3(maxu3, maxu);
    }
 
    virtual bool GetEnclosingRegion(const CoordType &minu, const CoordType &maxu, DimsType &min, DimsType &max) const = 0;
 
    virtual float GetMissingValue() const;
 
    void SetMissingValue(float missing_value) { _missingValue = missing_value; };
 
    //
    void SetHasMissingValues(bool flag) { _hasMissing = flag; }
 
    //
    bool HasMissingData() const { return (_hasMissing); };
 
    //
    virtual bool HasInvertedCoordinateSystemHandiness() const { return (true); }
 
    virtual int GetInterpolationOrder() const { return _interpolationOrder; };
 
    virtual void SetInterpolationOrder(int order);
 
    //
    static DimsType Dims(const DimsType &min, const DimsType &max);
 
    virtual void GetUserCoordinates(const DimsType &indices, CoordType &coords) const = 0;
 
    virtual void GetUserCoordinates(const size_t indices[], double coords[]) const
    {
        DimsType  indices3 = {0, 0, 0};
        CoordType coords3 = {0.0, 0.0, 0.0};
        CopyToArr3(indices, GetNodeDimensions().size(), indices3);
        GetUserCoordinates(indices3, coords3);
        CopyFromArr3(coords3, coords);
    }
 
    virtual void GetUserCoordinates(const std::vector<size_t> &indices, std::vector<double> &coords) const
    {
        DimsType  indices3 = {0, 0, 0};
        CoordType coords3 = {0.0, 0.0, 0.0};
        CopyToArr3(indices, indices3);
        GetUserCoordinates(indices3, coords3);
        CopyFromArr3(coords3, coords);
    }
 
    virtual void GetUserCoordinates(size_t i, double &x, double &y, double &z) const;
    virtual void GetUserCoordinates(size_t i, size_t j, double &x, double &y, double &z) const;
    virtual void GetUserCoordinates(size_t i, size_t j, size_t k, double &x, double &y, double &z) const;
 
    virtual bool GetIndicesCell(const CoordType &coords, DimsType &indices) const = 0;
 
    //
    virtual bool GetIndicesCell(const double coords[3], size_t indices[3]) const
    {
        CoordType c3 = {0.0, 0.0, 0.0};
        DimsType  i3 = {0, 0, 0};
        CopyToArr3(coords, GetGeometryDim(), c3);
        bool status = GetIndicesCell(c3, i3);
        CopyFromArr3(i3, indices);
        return (status);
    }
 
    //
    virtual bool GetIndicesCell(const std::vector<double> &coords, std::vector<size_t> &indices) const
    {
        CoordType c3 = {0.0, 0.0, 0.0};
        DimsType  i3 = {0, 0, 0};
        CopyToArr3(coords, c3);
        bool status = GetIndicesCell(c3, i3);
        CopyFromArr3(i3, indices);
        return (status);
    }
 
    virtual void GetRange(float range[2]) const;
 
    virtual void GetRange(const DimsType &min, const DimsType &max, float range[2]) const;
 
    //
    virtual void GetRange(std::vector<size_t> min, std::vector<size_t> max, float range[2]) const
    {
        DimsType min3 = {0, 0, 0};
        DimsType max3 = {0, 0, 0};
        CopyToArr3(min, min3);
        CopyToArr3(max, max3);
        GetRange(min3, max3, range);
    }
 
    virtual bool InsideGrid(const CoordType &coords) const = 0;
 
    //
    virtual bool InsideGrid(const double coords[3]) const
    {
        CoordType c3 = {0.0, 0.0, 0.0};
        CopyToArr3(coords, GetGeometryDim(), c3);
        return (InsideGrid(c3));
    }
 
    //
    virtual bool InsideGrid(const std::vector<double> &coords) const
    {
        CoordType c3 = {0.0, 0.0, 0.0};
        CopyToArr3(coords, c3);
        return (InsideGrid(c3));
    }
 
    virtual bool GetCellNodes(const DimsType &cindices, std::vector<DimsType> &nodes) const = 0;
 
    //
    virtual bool GetCellNodes(const size_t cindices[], std::vector<DimsType> &nodes) const
    {
        DimsType i3 = {0, 0, 0};
        CopyToArr3(cindices, GetCellDimensions().size(), i3);
        return (GetCellNodes(i3, nodes));
    }
 
    virtual bool GetCellNeighbors(const DimsType &cindices, std::vector<DimsType> &cells) const = 0;
 
    virtual bool GetNodeCells(const DimsType &indices, std::vector<DimsType> &cells) const = 0;
 
    virtual size_t GetMaxVertexPerFace() const = 0;
 
    virtual size_t GetMaxVertexPerCell() const = 0;
 
    //
    virtual void ClampCoord(const CoordType &coords, CoordType &cCoords) const = 0;
 
    //
    virtual void ClampCoord(const double coords[3], double cCoords[3]) const
    {
        CoordType c3 = {coords[0], coords[1], coords[2]};
        CoordType cC3;
        ClampCoord(c3, cC3);
        cCoords[0] = cC3[0];
        cCoords[1] = cC3[1];
        cCoords[2] = cC3[2];
    }
 
    //
    virtual void ClampIndex(const DimsType &indices, DimsType &cIndices) const { ClampIndex(GetNodeDimensions(), indices, cIndices); }
 
    //
    virtual void ClampCellIndex(const DimsType &indices, DimsType &cIndices) const { ClampIndex(GetCellDimensions(), indices, cIndices); }
 
    //
    virtual void SetPeriodic(const std::vector<bool> &periodic)
    {
        _periodic.clear();
        int i = 0;
        for (; i < periodic.size() && i < GetGeometryDim(); i++) { _periodic.push_back(periodic[i]); }
        for (; i < GetGeometryDim(); i++) { _periodic.push_back(false); }
    }
 
    //
    virtual const std::vector<bool> &GetPeriodic() const { return (_periodic); }
 
    //
    virtual size_t GetTopologyDim() const { return (_topologyDimension); }
 
    //
    virtual long GetNodeOffset() const { return (_nodeIDOffset); }
    virtual void SetNodeOffset(long offset) { _nodeIDOffset = offset; }
 
    //
    virtual long GetCellOffset() const { return (_cellIDOffset); }
    virtual void SetCellOffset(long offset) { _cellIDOffset = offset; }
 
    virtual DimsType GetMinAbs() const { return (_minAbs); }
 
    //
    virtual void SetMinAbs(const DimsType &minAbs) { _minAbs = minAbs; }
 
    //
    static bool PointInsideBoundingRectangle(const double pt[], const double verts[], int n)
    {
        VAssert(n > 2);
 
        double left = verts[0];
        double right = verts[0];
        double top = verts[1];
        double bottom = verts[1];
 
        for (int i = 1; i < n; i++) {
            if (verts[i * 2 + 0] < left) left = verts[i * 2 + 0];
            if (verts[i * 2 + 0] > right) right = verts[i * 2 + 0];
            if (verts[i * 2 + 1] < top) top = verts[i * 2 + 1];
            if (verts[i * 2 + 1] > bottom) bottom = verts[i * 2 + 1];
        }
 
        return ((left <= pt[0]) && (right >= pt[0]) && (top <= pt[1]) && (bottom >= pt[1]));
    }
 
 
    VDF_API friend std::ostream &operator<<(std::ostream &o, const Grid &g);
 
    //
    // Iterators
    //
 
    //
    class InsideBox {
    public:
        InsideBox(const CoordType &min, const CoordType &max) : _min(min), _max(max) {}
        InsideBox() {}
 
        bool operator()(const CoordType &pt) const
        {
            if (_min == _max) return (true);
            for (int i = 0; i < _min.size() && i < pt.size(); i++) {
                if (pt[i] < _min[i] || pt[i] > _max[i]) return (false);
            }
            return (true);
        }
        bool operator()(const double pt[]) const
        {
            if (_min == _max) return (true);
            for (int i = 0; i < _min.size(); i++) {
                if (pt[i] < _min[i] || pt[i] > _max[i]) return (false);
            }
            return (true);
        }
 
        bool Enabled() const { return (_min != _max); }
 
    private:
        CoordType _min;
        CoordType _max;
    };
 
    //
    // Define polymorphic iterator that can be used with any
    // class derived from this class
    //
    //
 
    // Interface for iterator specializations
    //
    template<typename T> class AbstractIterator {
    public:
        virtual ~AbstractIterator() {}
        virtual void                              next() = 0;
        virtual void                              next(const long &offset) = 0;
        virtual T &                               deref() const = 0;
        virtual const void *                      address() const = 0;
        virtual bool                              equal(const void *other) const = 0;
        virtual std::unique_ptr<AbstractIterator> clone() const = 0;
    };
 
    // Overloaded operators that will act on spealizations of
    // AbstractIterator
    //
    template<typename T> class PolyIterator {
    public:
        PolyIterator(std::unique_ptr<AbstractIterator<T>> it) : _impl(std::move(it)) {}
 
        PolyIterator(PolyIterator const &rhs) : _impl(rhs._impl->clone()) {}
 
        PolyIterator &operator=(PolyIterator const &rhs)
        {
            _impl = rhs._impl->clone();
            return *this;
        }
 
        // PolyIterator has a unique_ptr member so we must provide
        // std::move constructors
        //
        PolyIterator(PolyIterator &&rhs) { _impl = std::move(rhs._impl); }
 
        PolyIterator &operator=(PolyIterator &&rhs)
        {
            if (this != &rhs) { _impl = std::move(rhs._impl); }
            return (*this);
        }
 
        PolyIterator() : _impl(nullptr) {}
 
        PolyIterator &operator++()
        {    // ++prefix
            _impl->next();
            return *this;
        };
 
        PolyIterator operator++(int)
        {    // postfix++
            PolyIterator temp(*this);
            ++(*this);
            return (temp);
        };
 
        PolyIterator &operator+=(const long &offset)
        {
            _impl->next(offset);
            return (*this);
        };
 
        PolyIterator operator+(const long &offset) const
        {
            PolyIterator temp(*this);
            temp += offset;
            return (temp);
        }
 
        const T &operator*() const { return _impl->deref(); }
 
        bool operator==(const PolyIterator &rhs) const { return (_impl->equal(rhs._impl->address())); }
 
        bool operator!=(const PolyIterator &rhs) const { return (!(*this == rhs)); }
 
    private:
        std::unique_ptr<AbstractIterator<T>> _impl;
    };
 
    //
    typedef const CoordType                        ConstCoordType;
    typedef Grid::PolyIterator<ConstCoordType>     ConstCoordItr;
    typedef Grid::AbstractIterator<ConstCoordType> ConstCoordItrAbstract;
 
    //
    virtual ConstCoordItr ConstCoordBegin() const = 0;
    virtual ConstCoordItr ConstCoordEnd() const = 0;
 
    typedef const DimsType                         ConstIndexType;
    typedef Grid::PolyIterator<ConstIndexType>     ConstNodeIterator;
    typedef Grid::AbstractIterator<ConstIndexType> ConstNodeIteratorAbstract;
 
    //
    typedef Grid::PolyIterator<ConstIndexType>     ConstCellIterator;
    typedef Grid::AbstractIterator<ConstIndexType> ConstCellIteratorAbstract;
 
    //
    // Iterators
    //
 
    //
    // Node index iterator. Iterates over node indices
    //
    class ConstNodeIteratorSG : public Grid::ConstNodeIteratorAbstract {
    public:
        ConstNodeIteratorSG(const Grid *g, bool begin);
        ConstNodeIteratorSG(const ConstNodeIteratorSG &rhs);
        ConstNodeIteratorSG();
 
        virtual ~ConstNodeIteratorSG() {}
 
        virtual void            next();
        virtual void            next(const long &offset);
        virtual ConstIndexType &deref() const { return (_index); }
        virtual const void *    address() const { return this; };
 
        virtual bool equal(const void *rhs) const
        {
            const ConstNodeIteratorSG *itrptr = static_cast<const ConstNodeIteratorSG *>(rhs);
            return (_index == itrptr->_index);
        }
 
        virtual std::unique_ptr<ConstNodeIteratorAbstract> clone() const { return std::unique_ptr<ConstNodeIteratorAbstract>(new ConstNodeIteratorSG(*this)); };
 
    protected:
        DimsType _dims;
        DimsType _index;
        DimsType _lastIndex;
    };
 
    class ConstNodeIteratorBoxSG : public ConstNodeIteratorSG {
    public:
        ConstNodeIteratorBoxSG(const Grid *g, const CoordType &minu, const CoordType &maxu);
        ConstNodeIteratorBoxSG(const ConstNodeIteratorBoxSG &rhs);
        ConstNodeIteratorBoxSG();
 
        virtual ~ConstNodeIteratorBoxSG() {}
 
        virtual void next();
        virtual void next(const long &offset);
 
        virtual std::unique_ptr<ConstNodeIteratorAbstract> clone() const { return std::unique_ptr<ConstNodeIteratorAbstract>(new ConstNodeIteratorBoxSG(*this)); };
 
    private:
        InsideBox     _pred;
        ConstCoordItr _coordItr;
    };
 
    virtual ConstNodeIterator ConstNodeBegin() const { return ConstNodeIterator(std::unique_ptr<ConstNodeIteratorAbstract>(new ConstNodeIteratorSG(this, true))); }
 
    virtual ConstNodeIterator ConstNodeBegin(const CoordType &minu, const CoordType &maxu) const
    {
        return ConstNodeIterator(std::unique_ptr<ConstNodeIteratorAbstract>(new ConstNodeIteratorBoxSG(this, minu, maxu)));
    }
 
    virtual ConstNodeIterator ConstNodeBegin(const std::vector<double> &minu, const std::vector<double> &maxu) const
    {
        CoordType minuCT, maxuCT;
        CopyToArr3(minu, minuCT);
        CopyToArr3(maxu, maxuCT);
        return ConstNodeBegin(minuCT, maxuCT);
    }
 
    virtual ConstNodeIterator ConstNodeEnd() const { return ConstNodeIterator(std::unique_ptr<ConstNodeIteratorAbstract>(new ConstNodeIteratorSG(this, false))); }
 
    //
    // Cell index iterator. Iterates over cell indices
    //
    class ConstCellIteratorSG : public Grid::ConstCellIteratorAbstract {
    public:
        ConstCellIteratorSG(const Grid *g, bool begin);
        ConstCellIteratorSG(const ConstCellIteratorSG &rhs);
        ConstCellIteratorSG();
 
        virtual ~ConstCellIteratorSG() {}
 
        virtual void            next();
        virtual void            next(const long &offset);
        virtual ConstIndexType &deref() const { return (_index); }
        virtual const void *    address() const { return this; };
 
        virtual bool equal(const void *rhs) const
        {
            const ConstCellIteratorSG *itrptr = static_cast<const ConstCellIteratorSG *>(rhs);
 
            return (_index == itrptr->_index);
        }
 
        virtual std::unique_ptr<ConstCellIteratorAbstract> clone() const { return std::unique_ptr<ConstCellIteratorAbstract>(new ConstCellIteratorSG(*this)); };
 
    protected:
        DimsType _dims;
        DimsType _index;
        DimsType _lastIndex;
    };
 
    class ConstCellIteratorBoxSG : public ConstCellIteratorSG {
    public:
        ConstCellIteratorBoxSG(const Grid *g, const CoordType &minu, const CoordType &maxu);
        ConstCellIteratorBoxSG(const ConstCellIteratorBoxSG &rhs);
        ConstCellIteratorBoxSG();
 
        virtual ~ConstCellIteratorBoxSG() {}
 
        virtual void next();
        virtual void next(const long &offset);
 
        virtual std::unique_ptr<ConstCellIteratorAbstract> clone() const { return std::unique_ptr<ConstCellIteratorAbstract>(new ConstCellIteratorBoxSG(*this)); };
 
    private:
        InsideBox _pred;
        const Grid *_g;
        bool        _cellInsideBox(const size_t cindices[]) const;
    };
 
    virtual ConstCellIterator ConstCellBegin() const { return ConstCellIterator(std::unique_ptr<ConstCellIteratorAbstract>(new ConstCellIteratorSG(this, true))); }
 
    virtual ConstCellIterator ConstCellBegin(const CoordType &minu, const CoordType &maxu) const
    {
        return ConstCellIterator(std::unique_ptr<ConstCellIteratorAbstract>(new ConstCellIteratorBoxSG(this, minu, maxu)));
    }
 
    virtual ConstCellIterator ConstCellBegin(const std::vector<double> &minu, const std::vector<double> &maxu) const
    {
        CoordType minuCT, maxuCT;
        CopyToArr3(minu, minuCT);
        CopyToArr3(maxu, maxuCT);
        return ConstCellBegin(minuCT, maxuCT);
    }
 
    virtual ConstCellIterator ConstCellEnd() const { return ConstCellIterator(std::unique_ptr<ConstCellIteratorAbstract>(new ConstCellIteratorSG(this, false))); }
 
    //
    template<class T> class VDF_API ForwardIterator {
    public:
        ForwardIterator(T *rg, bool begin = true, const CoordType &minu = {0.0, 0.0, 0.0}, const CoordType &maxu = {0.0, 0.0, 0.0});
        ForwardIterator();
        ForwardIterator(const ForwardIterator<T> &) = default;
        ForwardIterator(ForwardIterator<T> &&rhs);
        ~ForwardIterator() {}
 
        float &operator*() { return (*_itr); }
 
        ForwardIterator<T> &operator++();    // ++prefix
 
        ForwardIterator<T> operator++(int);    // postfix++
 
        ForwardIterator<T> &operator+=(const long int &offset);
        ForwardIterator<T>  operator+(const long int &offset) const;
 
        ForwardIterator<T> &operator=(ForwardIterator<T> rhs);
        ForwardIterator<T> &operator=(ForwardIterator<T> &rhs) = delete;
 
        bool operator==(const ForwardIterator<T> &rhs) const { return (_index == rhs._index); }
        bool operator!=(const ForwardIterator<T> &rhs) { return (_index != rhs._index); }
 
        const ConstCoordItr &GetCoordItr() { return (_coordItr); }
 
        friend void swap(Grid::ForwardIterator<T> &a, Grid::ForwardIterator<T> &b)
        {
            std::swap(a._blks, b._blks);
            std::swap(a._dims3d, b._dims3d);
            std::swap(a._bdims3d, b._bdims3d);
            std::swap(a._bs3d, b._bs3d);
            std::swap(a._blocksize, b._blocksize);
            std::swap(a._coordItr, b._coordItr);
            std::swap(a._index, b._index);
            std::swap(a._lastIndex, b._lastIndex);
            std::swap(a._xb, b._xb);
            std::swap(a._itr, b._itr);
            std::swap(a._pred, b._pred);
        }
 
    private:
        std::vector<float *> _blks;
        DimsType             _dims3d;
        DimsType             _bdims3d;
        DimsType             _bs3d;
        size_t               _blocksize;
        ConstCoordItr        _coordItr;
        DimsType             _index;        // current index into grid
        DimsType             _lastIndex;    // Last valid index
        size_t               _xb;           // x index within a block
        float *              _itr;
        InsideBox            _pred;
    };
 
    typedef Grid::ForwardIterator<Grid>       Iterator;
    typedef Grid::ForwardIterator<Grid const> ConstIterator;
 
    //
    Iterator begin(const CoordType &minu, const CoordType &maxu) { return (Iterator(this, true, minu, maxu)); }
 
    Iterator begin(const std::vector<double> &minu, const std::vector<double> &maxu)
    {
        CoordType minuCT, maxuCT;
        CopyToArr3(minu, minuCT);
        CopyToArr3(maxu, maxuCT);
        return begin(minuCT, maxuCT);
    }
 
    Iterator begin() { return (Iterator(this, true)); }
 
    Iterator end() { return (Iterator(this, false)); }
 
    ConstIterator cbegin(const CoordType &minu, const CoordType &maxu) const { return (ConstIterator(this, true, minu, maxu)); }
 
    ConstIterator cbegin(const std::vector<double> &minu, const std::vector<double> &maxu)
    {
        CoordType minuCT, maxuCT;
        CopyToArr3(minu, minuCT);
        CopyToArr3(maxu, maxuCT);
        return cbegin(minuCT, maxuCT);
    }
    ConstIterator cbegin() const { return (ConstIterator(this, true)); }
 
    ConstIterator cend() const { return (ConstIterator(this, false)); }
 
    template<typename T> static void CopyToArr3(const std::vector<T> &src, std::array<T, 3> &dst)
    {
        for (int i = 0; i < src.size() && i < dst.size(); i++) { dst[i] = src[i]; }
    }
 
    template<typename T> static void CopyToArr3(const T *src, size_t n, std::array<T, 3> &dst)
    {
        for (int i = 0; i < n && i < dst.size(); i++) { dst[i] = src[i]; }
    }
    template<typename T> static void CopyFromArr3(const std::array<T, 3> &src, std::vector<T> &dst)
    {
        dst.resize(src.size());
        for (int i = 0; i < src.size() && i < dst.size(); i++) { dst[i] = src[i]; }
    }
    template<typename T> static void CopyFromArr3(const std::array<T, 3> &src, T *dst)
    {
        for (int i = 0; i < src.size(); i++) { dst[i] = src[i]; }
    }
 
protected:
    virtual float GetValueNearestNeighbor(const CoordType &coords) const = 0;
 
    virtual float GetValueLinear(const CoordType &coords) const = 0;
 
    virtual void GetUserExtentsHelper(CoordType &minu, CoordType &maxu) const = 0;
 
    virtual float *GetValuePtrAtIndex(const std::vector<float *> &blks, const DimsType &indices) const;
 
    // Deprecated
    virtual void ClampIndex(const std::vector<size_t> &dims, const DimsType indices, DimsType &cIndices) const
    {
        cIndices = {0, 0, 0};
 
        for (int i = 0; i < dims.size(); i++) {
            cIndices[i] = indices[i];
            if (cIndices[i] >= dims[i]) { cIndices[i] = dims[i] - 1; }
        }
    }
 
    virtual void ClampIndex(const DimsType &dims, const DimsType indices, DimsType &cIndices) const
    {
        cIndices = {0, 0, 0};
 
        for (int i = 0; i < dims.size(); i++) {
            cIndices[i] = indices[i];
            if (cIndices[i] >= dims[i]) {
                assert(dims[i] > 0);
                cIndices[i] = dims[i] - 1;
            }
        }
    }
 
    float BilinearInterpolate(size_t i, size_t j, size_t k, const double xwgt, const double ywgt) const;
 
    float TrilinearInterpolate(size_t i, size_t j, size_t k, const double xwgt, const double ywgt, const double zwgt) const;
 
private:
    DimsType             _dims;                   // dimensions of grid arrays
    DimsType             _bs = {{1, 1, 1}};       // dimensions of each block
    DimsType             _bdims = {{1, 1, 1}};    // dimensions (specified in blocks) of ROI
    std::vector<size_t>  _bsDeprecated;           // legacy API
    std::vector<size_t>  _bdimsDeprecated;        // legacy API
    std::vector<float *> _blks;
    std::vector<bool>    _periodic;    // periodicity of boundaries
    DimsType             _minAbs;      // Offset to start of grid
    size_t               _topologyDimension = 0;
    float                _missingValue = std::numeric_limits<float>::infinity();
    bool                 _hasMissing = false;
    int                  _interpolationOrder = 1;    // Order of interpolation
    long                 _nodeIDOffset = 0;
    long                 _cellIDOffset = 0;
    mutable CoordType    _minuCache = {{std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity()}};
    mutable CoordType    _maxuCache = {{std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity()}};
 
    void _grid(const DimsType &dims, const DimsType &bs, const std::vector<float *> &blks, size_t topology_dimension);
};
 
template void Grid::CopyToArr3<size_t>(const std::vector<size_t> &src, std::array<size_t, 3> &dst);
 
};    // namespace VAPoR
#endif