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 _CLASSIC_FIELDS_NDGRID_HH_
#define _CLASSIC_FIELDS_NDGRID_HH_

#include <vector>
#include <algorithm>
#include <math.h>
#include <ostream>

#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() {;}

    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>(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 // _CLASSIC_FIELDS_NDGRID_HH_