ThreeDGrid.h

Go to the documentation of this file.

/* 
 *  Copyright (c) 2012, 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_THREEDGRID_HH_
#define _CLASSIC_FIELDS_THREEDGRID_HH_

#include <math.h>

#include <algorithm>
#include <vector>

#include "Fields/Interpolation/Mesh.h"

namespace interpolation {

class ThreeDGrid : public Mesh {
  public:
    class Iterator;

    ThreeDGrid * clone() {return new ThreeDGrid(*this);}

    Mesh * dual () const;

    ThreeDGrid();

    ThreeDGrid(double dX, double dY, double dZ,
               double minX, double minY, double minZ,
               int numberOfXCoords, int numberOfYCoords, int numberOfZCoords);

    ThreeDGrid(int xSize, const double *x,
               int ySize, const double *y,
               int zSize, const double *z);

    ThreeDGrid(std::vector<double> x,
               std::vector<double> y,
               std::vector<double> z);

    ~ThreeDGrid();

    inline double& x(const int& i) {return x_m[i-1];}

    inline double& y(const int& j) {return y_m[j-1];}

    inline double& z(const int& k) {return z_m[k-1];}

    inline const double& x(const int& i) const {return x_m[i-1];}

    inline const double& y(const int& j) const {return y_m[j-1];}

    inline const double& z(const int& j) const {return z_m[j-1];}

    inline int     xSize() const {return static_cast<int>(x_m.size());}

    inline int     ySize() const {return static_cast<int>(y_m.size());}

    inline int     zSize() const {return static_cast<int>(z_m.size());}

    std::vector<double> xVector() {return std::vector<double>(x_m);}

    std::vector<double> yVector() {return std::vector<double>(y_m);}

    std::vector<double> zVector() {return std::vector<double>(z_m);}

    inline double* newXArray();

    inline double* newYArray();

    inline double* newZArray();

    inline void xLowerBound(const double& x, int& xIndex) const;

    inline void yLowerBound(const double& y, int& yIndex) const;

    inline void zLowerBound(const double& z, int& zIndex) const;

    inline void lowerBound(const double& x, int& xIndex,
                           const double& y, int& yIndex,
                           const double& z, int& zIndex) const;

    inline void lowerBound(const double& x, const double& y, const double& z,
                           Mesh::Iterator& it) const;

    inline double minX() const;

    inline double maxX() const;

    inline double minY() const;

    inline double maxY() const;

    inline double minZ() const;

    inline double maxZ() const;

    inline void setX(int nXCoords, double * x);

    inline void setY(int nYCoords, double * y);

    inline void setZ(int nZCoords, double * z);

    void add(VectorMap* map);

    void remove(VectorMap* map);

    Mesh::Iterator begin() const;

    Mesh::Iterator end()   const;

    virtual void getPosition(const Mesh::Iterator& it, double * position) const;

    inline int getPositionDimension() const;

    inline int toInteger(const Mesh::Iterator& lhs) const;

    inline void setConstantSpacing(bool spacing);

    void setConstantSpacing();

    inline bool getConstantSpacing() const;

    Mesh::Iterator getNearest(const double* position) const;

    static void vectorLowerBound(std::vector<double> vec, double x, int& index);

  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<double>     x_m;
    std::vector<double>     y_m;
    std::vector<double>     z_m;
    int                     xSize_m;
    int                     ySize_m;
    int                     zSize_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* ThreeDGrid::newXArray() {
    double *x = new double[x_m.size()];
    for (unsigned int i = 0; i < x_m.size(); i++)
        x[i] = x_m[i];
    return x;
}

double* ThreeDGrid::newYArray() {
    double *y = new double[y_m.size()];
    for (unsigned int i = 0; i < y_m.size(); i++)
        y[i] = y_m[i];
    return y;
}

double* ThreeDGrid::newZArray() {
    double *z = new double[z_m.size()];
    for (unsigned int i = 0; i < z_m.size(); i++)
        z[i] = z_m[i];
    return z;
}

void ThreeDGrid::lowerBound(const double& x, int& xIndex,
                            const double& y, int& yIndex,
                            const double& z, int& zIndex) const {
    xLowerBound(x, xIndex);
    yLowerBound(y, yIndex);
    zLowerBound(z, zIndex);
}

void ThreeDGrid::lowerBound(const double& x,
                            const double& y,
                            const double& z,
                            Mesh::Iterator& it) const {
    xLowerBound(x, it[0]);
    yLowerBound(y, it[1]);
    zLowerBound(z, it[2]);
    it[0]++;
    it[1]++;
    it[2]++;
}

void ThreeDGrid::xLowerBound(const double& x, int& xIndex) const {
    if (constantSpacing_m)
        xIndex = static_cast<int>(floor((x - x_m[0])/(x_m[1]-x_m[0]) ));
    else {
        vectorLowerBound(x_m, x, xIndex);
    }
}


void ThreeDGrid::yLowerBound(const double& y, int& yIndex) const {
    if (constantSpacing_m)
        yIndex = static_cast<int>(floor((y - y_m[0])/(y_m[1]-y_m[0])));
    else
        vectorLowerBound(y_m, y, yIndex);
}

void ThreeDGrid::zLowerBound(const double& z, int& zIndex) const {
    if (constantSpacing_m)
        zIndex = static_cast<int>(floor((z - z_m[0])/(z_m[1]-z_m[0])));
    else
        vectorLowerBound(z_m, z, zIndex);
}

double ThreeDGrid::minX() const {
    return x_m[0];
}

double ThreeDGrid::maxX() const {
    return x_m[xSize_m-1];
}

double ThreeDGrid::minY() const {
    return y_m[0];
}

double ThreeDGrid::maxY() const {
    return y_m[ySize_m-1];
}

double ThreeDGrid::minZ() const {
    return z_m[0];
}

double ThreeDGrid::maxZ() const {
    return z_m[zSize_m-1];
}

void ThreeDGrid::setX(int nXCoords, double * x) {
    x_m = std::vector<double>(x, x+nXCoords);
}

void ThreeDGrid::setY(int nYCoords, double * y) {
    y_m = std::vector<double>(y, y+nYCoords);
}

void ThreeDGrid::setZ(int nZCoords, double * z) {
    z_m = std::vector<double>(z, z+nZCoords);
}

int ThreeDGrid::getPositionDimension() const {
    return 3;
}

int ThreeDGrid::toInteger(const Mesh::Iterator& lhs) const {
    return (lhs.getState()[0]-1)*zSize_m*ySize_m+
           (lhs.getState()[1]-1)*zSize_m+
           (lhs.getState()[2]-1);
}

void ThreeDGrid::setConstantSpacing(bool spacing) {
    constantSpacing_m = spacing;
}

bool ThreeDGrid::getConstantSpacing() const {
    return constantSpacing_m;
}
}
#endif