NDGrid.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2017, Chris Rogers
 *  All rights reserved.
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *  1. Redistributions of source code must retain the above copyright notice,
 *     this list of conditions and the following disclaimer.
 *  2. Redistributions in binary form must reproduce the above copyright notice,
 *     this list of conditions and the following disclaimer in the documentation
 *     and/or other materials provided with the distribution.
 *  3. Neither the name of STFC nor the names of its contributors may be used to
 *     endorse or promote products derived from this software without specific
 *     prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 */
 
#ifndef _OPALX_FIELDS_NDGRID_HH_
#define _OPALX_FIELDS_NDGRID_HH_
 
#include <algorithm>
#include <cmath>
#include <ostream>
#include <vector>
 
#include "Fields/Interpolation/Mesh.h"
 
namespace interpolation {
 
    class NDGrid : public Mesh {
    public:
        class Iterator;
 
        inline NDGrid* clone();
 
        Mesh* dual() const;
 
        NDGrid();
 
        NDGrid(int nDimensions, int* size, double* spacing, double* min);
 
        NDGrid(std::vector<int> size, std::vector<const double*> gridCoordinates);
 
        NDGrid(std::vector<std::vector<double> > gridCoordinates);
 
        NDGrid(const NDGrid& grid);
 
        ~NDGrid() { ; }
 
        inline double& coord(const int& index, const int& dimension);
 
        inline const double& coord(const int& index, const int& dimension) const;
 
        inline int size(const int& dimension) const;
 
        inline std::vector<double> coordVector(const int& dimension) const;
 
        double* newCoordArray(const int& dimension) const;
 
        inline void coordLowerBound(const double& x, const int& dimension, int& xIndex) const;
 
        inline void lowerBound(const std::vector<double>& pos, std::vector<int>& xIndex) const;
 
        inline double min(const int& dimension) const;
 
        inline double max(const int& dimension) const;
 
        inline void setCoord(int dimension, int nCoords, double* x);
 
        inline Mesh::Iterator begin() const;
 
        inline Mesh::Iterator end() const;
 
        inline void getPosition(const Mesh::Iterator& it, double* position) const;
 
        inline int getPositionDimension() const;
 
        inline bool getConstantSpacing() const;
 
        inline void setConstantSpacing(bool spacing);
 
        void setConstantSpacing(double tolerance_m = 1e-9);
 
        int toInteger(const Mesh::Iterator& lhs) const;
 
        Mesh::Iterator getNearest(const double* position) const;
 
    protected:
        // Change position
        virtual Mesh::Iterator& addEquals(Mesh::Iterator& lhs, int difference) const;
        virtual Mesh::Iterator& subEquals(Mesh::Iterator& lhs, int difference) const;
        virtual Mesh::Iterator& addEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
        virtual Mesh::Iterator& subEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
        virtual Mesh::Iterator& addOne(Mesh::Iterator& lhs) const;
        virtual Mesh::Iterator& subOne(Mesh::Iterator& lhs) const;
        // Check relative position
        virtual bool isGreater(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
 
    private:
        std::vector<std::vector<double> > coord_m;
        std::vector<VectorMap*> maps_m;
        bool constantSpacing_m;
 
        friend Mesh::Iterator operator++(Mesh::Iterator& lhs, int);
        friend Mesh::Iterator operator--(Mesh::Iterator& lhs, int);
        friend Mesh::Iterator& operator++(Mesh::Iterator& lhs);
        friend Mesh::Iterator& operator--(Mesh::Iterator& lhs);
        friend Mesh::Iterator operator-(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
        friend Mesh::Iterator operator+(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
        friend Mesh::Iterator& operator-=(Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
        friend Mesh::Iterator& operator+=(Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
 
        friend bool operator==(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
        friend bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
        friend bool operator>=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
        friend bool operator<=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
        friend bool operator<(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
        friend bool operator>(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
    };
 
    double& NDGrid::coord(const int& index, const int& dimension) {
        return coord_m[dimension][index - 1];
    }
 
    const double& NDGrid::coord(const int& index, const int& dimension) const {
        return coord_m[dimension][index - 1];
    }
 
    int NDGrid::size(const int& dimension) const { return coord_m[dimension].size(); }
 
    std::vector<double> NDGrid::coordVector(const int& dimension) const {
        return coord_m[dimension];
    }
 
    void NDGrid::coordLowerBound(const double& x, const int& dimension, int& xIndex) const {
        if (constantSpacing_m) {
            double x0     = coord_m[dimension][0];
            double x1     = coord_m[dimension][1];
            xIndex        = static_cast<int>(std::floor((x - x0) / (x1 - x0)));
            int coordSize = static_cast<int>(coord_m[dimension].size());
            if (xIndex >= coordSize) {
                xIndex = coordSize - 1;
            } else if (xIndex < -1) {
                xIndex = -1;
            }
        } else {
            const std::vector<double>& c_t(coord_m[dimension]);
            xIndex = std::lower_bound(c_t.begin(), c_t.end(), x) - c_t.begin() - 1;
        }
    }
 
    void NDGrid::lowerBound(const std::vector<double>& pos, std::vector<int>& xIndex) const {
        xIndex = std::vector<int>(pos.size());
        for (unsigned int i = 0; i < pos.size(); i++) {
            coordLowerBound(pos[i], i, xIndex[i]);
        }
    }
 
    double NDGrid::min(const int& dimension) const { return coord_m[dimension][0]; }
 
    double NDGrid::max(const int& dimension) const {
        return coord_m[dimension][coord_m[dimension].size() - 1];
    }
 
    NDGrid* NDGrid::clone() { return new NDGrid(*this); }
 
    void NDGrid::setCoord(int dimension, int nCoords, double* x) {
        coord_m[dimension] = std::vector<double>(x, x + nCoords);
    }
 
    Mesh::Iterator NDGrid::begin() const {
        return Mesh::Iterator(std::vector<int>(coord_m.size(), 1), this);
    }
 
    Mesh::Iterator NDGrid::end() const {
        if (coord_m.size() == 0) {
            return Mesh::Iterator(std::vector<int>(), this);
        }
        std::vector<int> end(coord_m.size(), 1);
        end[0] = coord_m[0].size() + 1;
        return Mesh::Iterator(end, this);
    }
 
    void NDGrid::getPosition(const Mesh::Iterator& it, double* position) const {
        for (unsigned int i = 0; i < it.getState().size(); i++)
            position[i] = coord_m[i][it[i] - 1];
    }
 
    int NDGrid::getPositionDimension() const { return coord_m.size(); }
 
    bool NDGrid::getConstantSpacing() const { return constantSpacing_m; }
 
    void NDGrid::setConstantSpacing(bool spacing) { constantSpacing_m = spacing; }
 
}  // namespace interpolation
 
#endif  // _OPALX_FIELDS_NDGRID_HH_