GridBase.h

 
#pragma once
 
#include <cmath>
#include <exception>
#include <vector>
 
namespace SSAGES
{
 
 
template<typename T>
class GridBase
{
protected:
    std::vector<T> data_;
 
    size_t dimension_;
 
    std::vector<int> numPoints_;
 
    std::pair< std::vector<double>,std::vector<double> > edges_;
 
    std::vector<bool> isPeriodic_;
 
 
    std::vector<int> wrapIndices(const std::vector<int> &indices) const
    {
        std::vector<int> newIndices(indices);
        for (size_t i = 0; i < dimension_; ++i) {
            if (!GetPeriodic(i)) {
                continue;
            }
 
            int index = indices.at(i);
            while (index < 0) { index += numPoints_[i]; }
            while (index >= numPoints_[i]) { index -= numPoints_[i]; }
        if(index == numPoints_[i]-1)
        index = 0;
            newIndices.at(i) = index;
        }
 
        return newIndices;
    }
 
 
    virtual size_t mapTo1d(const std::vector<int> &indices) const = 0;
 
protected:
 
    GridBase(std::vector<int> numPoints,
             std::vector<double> lower,
             std::vector<double> upper,
             std::vector<bool> isPeriodic)
        : dimension_(numPoints.size()),
          numPoints_(numPoints),
          edges_(std::pair< std::vector<double>, std::vector<double> >(lower, upper)),
          isPeriodic_(isPeriodic)
    {
        // Check that vector sizes are correct
        if (edges_.first.size() != dimension_ ||
            edges_.second.size() != dimension_) {
            throw std::invalid_argument("Size of vector containing upper or "
                "lower edges, does not match size of vector containing "
                "number of grid points.");
        }
        if (isPeriodic_.size() == 0) {
            // Default: Non-periodic in all dimensions
            isPeriodic.resize(dimension_, false);
        } else if (isPeriodic_.size() != dimension_) {
            throw std::invalid_argument("Size of vector isPeriodic does not "
                    "match size of vector containing number of grid points.");
        }
        for(size_t i=0 ; i < isPeriodic_.size() ; i++)
        {
            if(isPeriodic_[i])
            {
                if(numPoints_[i] <= 1)
                {
                    throw std::invalid_argument("A periodic grid is incompatible with a grid size of 1.");
                }               
                double spacing = (edges_.second[i] - edges_.first[i]) / (numPoints_[i]-1);
                edges_.first[i] -= spacing/2;
                edges_.second[i] += spacing/2;
            }
        }
        
    }
 
public:
 
    void syncGrid()
    {
        //Convenience
        size_t dim = this->GetDimension();
        
        //Preallocate
        std::vector<int> navigate(dim);
        
        //Loop over all surfaces. Number of surfaces = dim.
        for(size_t i = 0 ; i < dim ; i++)
        {
            //Check if periodic in this dimension.
            if(isPeriodic_[i])
            {
            
                //Calculate surface size. This is equal to number of points in each other dimension multiplied.
                size_t surfsize = 1;
                for(size_t j = 0 ; j < dim ; j++)
                {
                    if(i != j)
                        surfsize*=numPoints_[j];
                }
                
                //Loop over all points of this surface on the 0 side and copy to end.
                for(size_t j = 0 ; j < surfsize ; j++)
                {
                    int runningcount = j;
                    for(size_t k = 0 ; k < dim ; k++)
                    {
                        if(k == i)
                            navigate[k] = 0;
                        else
                        {
                            navigate[k] = runningcount%numPoints_[k];
                            runningcount = runningcount / numPoints_[k];
                        }
                    }       
                    auto temp = this->at(navigate);
                    navigate[i] = numPoints_[i]-1;
                    this->at(navigate) = temp;  
                }
            }
        }
    }
                
 
    size_t GetDimension() const
    {
        return dimension_;
    }
 
 
    const std::vector<int> GetNumPoints() const
    {
        return numPoints_;
    }
 
 
    int GetNumPoints(size_t dim) const
    {
        if (dim >= GetDimension()) {
            std::cerr << "Warning! Grid size requested for a dimension larger "
                         "than the grid dimensionality!\n";
            return 0;
        }
    return numPoints_.at(dim);
    }
 
 
    const std::vector<double> GetLower() const
    {
        std::vector<double> lower(dimension_);
        for(size_t i = 0; i<dimension_;++i) 
            if(GetPeriodic(i)) {lower[i] = edges_.first[i] - ((edges_.first[i] - edges_.second[i]) / (numPoints_[i]))/2;}
            else {lower[i] = edges_.first[i];}
        return lower;
    }
 
 
    double GetLower(size_t dim) const
    {
        if (dim >= GetDimension()) {
            std::cerr << "Warning! Lower edge requested for a dimension larger "
                         "than the grid dimensionality!\n";
            return 0.0;
        }
        if(GetPeriodic(dim)){return edges_.first[dim] - ((edges_.first[dim] - edges_.second[dim]) / (numPoints_[dim]))/2;}
        else{return edges_.first[dim];}
    }
 
 
    const std::vector<double> GetUpper() const
    {
        std::vector<double> upper(dimension_);
        for(size_t i = 0; i<dimension_;++i) 
            if(GetPeriodic(i)) {upper[i] = edges_.second[i] + ((edges_.first[i] - edges_.second[i]) / (numPoints_[i]))/2;}
            else {upper[i] = edges_.second[i];}
        return upper;
    }
 
 
 
    double GetUpper(size_t dim) const
    {
        if (dim >= GetDimension()) {
            std::cerr << "Warning! Upper edge requested for a dimension larger "
                         "than the grid dimensionality!\n";
            return 0.0;
        }
        if(GetPeriodic(dim)){return edges_.second[dim] + ((edges_.first[dim] - edges_.second[dim]) / (numPoints_[dim]))/2;}
        else{return edges_.second[dim];}
    }
 
 
    const std::vector<bool>& GetPeriodic() const
    {
        return isPeriodic_;
    }
 
 
    bool GetPeriodic(size_t dim) const
    {
        if (dim >= GetDimension()) {
            std::cerr << "Warning! Periodicity requested for a dimension larger "
                         "than the grid dimensionality!\n";
            return false;
        }
        return GetPeriodic().at(dim);
    }
 
 
    size_t size() const
    {
        return data_.size();
    }
 
 
    T *data()
    {
        return data_.data();
    }
 
 
    T const *data() const
    {
        return data_.data();
    }
 
 
    std::vector<int> GetIndices(const std::vector<double> &x) const
    {
        // Check that input vector has the correct dimensionality
        if (x.size() != dimension_) {
            throw std::invalid_argument("Specified point has a larger "
                    "dimensionality than the grid.");
        }
 
        std::vector<int> indices(dimension_);
        for (size_t i = 0; i < dimension_; ++i) 
    {
            double xpos = x.at(i);
                if (!GetPeriodic(i)) 
        {
                    if (xpos < edges_.first[i]) 
            {
                            indices.at(i) = -1;
                            continue;
                    } 
            else if (xpos > edges_.second[i]) 
            {
                            indices.at(i) = numPoints_[i];
                            continue;
                    }
 
                }
 
                // To make sure, the value is rounded in the correct direction.
        double spacing = (edges_.second[i] - edges_.first[i]) / (numPoints_[i]);
                indices.at(i) = std::floor( (xpos - edges_.first[i]) / spacing);
        }
        
        return wrapIndices(indices);
    }
    
 
 
    int GetIndex(double x) const
    {
        if (dimension_ != 1) {
            throw std::invalid_argument("1d Grid index can only be accessed for "
                   "1d-Grids can be accessed with a.");
        }
        return GetIndices({x}).at(0);
    }
 
 
    /*/!
     *\param coords
     *\param dim
     * To be finalized.
        
    std::vector<double> GetDist(const std::vector<double> &coords)
    {
    
    }*/
 
 
 
    std::vector<double> GetCoordinates(const std::vector<int> &indices)
    {
        if (indices.size() != dimension_) {
            throw std::invalid_argument(
                    "Grid indices specified for GetCoordinates() do not have "
                    "the same dimensionality as the grid.");
        }
 
        std::vector<double> v(dimension_);
 
        for (size_t i = 0; i < dimension_; ++i) {
            double spacing = (edges_.second[i] - edges_.first[i]) / (numPoints_[i]);
            v.at(i) = edges_.first[i] + (indices[i] + 0.5)*spacing;
        }
 
        return v;
    }
 
 
    double GetCoordinate(int index)
    {
        return GetCoordinates({index}).at(0);
    }
 
 
    const T& at(const std::vector<int> &indices) const
    {
        // Check that indices are in bound.
        if (indices.size() != GetDimension()) {
            throw std::invalid_argument("Dimension of indices does not match "
                    "dimension of the grid.");
        }
 
        return data_.at(mapTo1d(indices));
    }
 
 
    T& at(const std::vector<int> &indices)
    {
        return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(indices));
    }
 
 
    template<typename R>
    const T& at(std::initializer_list<R>&& x) const
    {
        return at(static_cast<std::vector<R> >(x));
    }
 
 
    template<typename R>
    T& at(std::initializer_list<R>&& x)
    {
        return at(static_cast<std::vector<R> >(x));
    }
 
 
    const T& at(int index) const
    {
        if (dimension_ != 1) {
            throw std::invalid_argument("Only 1d-Grids can be accessed with a "
                    "single integer as the index.");
        }
        return at({index});
    }
 
 
    T& at(int index)
    {
        return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(index));
    }
 
 
    const T& at(const std::vector<double> &x) const
    {
        return at(GetIndices(x));
    }
 
 
    T& at(const std::vector<double> &x)
    {
        return at(GetIndices(x));
    }
 
 
    const T& at(double x) const
    {
        if (dimension_ != 1) {
            throw std::invalid_argument("Only 1d-Grids can be accessed with a "
                    "single float as the specified point.");
        }
        return at({x});
    }
 
 
    T& at(double x)
    {
        return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(x));
    }
 
 
    const T& operator[](const std::vector<int> &indices) const
    {
        return at(indices);
    }
 
 
    T& operator[](const std::vector<int> &indices)
    {
        return at(indices);
    }
 
 
    template<typename R>
    const T& operator[](std::initializer_list<R>&& x) const
    {
        return at(static_cast<std::vector<R> >(x));
    }
 
 
    template<typename R>
    T& operator[](std::initializer_list<R>&& x)
    {
        return at(static_cast<std::vector<R> >(x));
    }
 
 
    const T& operator[](int index) const
    {
        return at(index);
    }
 
 
    T& operator[](int index)
    {
        return at(index);
    }
 
 
    const T& operator[](const std::vector<double> &x) const
    {
        return at(x);
    }
 
 
    T& operator[](const std::vector<double> &x)
    {
        return at(x);
    }
 
 
    const T& operator[](double x) const
    {
        return at(x);
    }
 
 
    T& operator[](double x)
    {
        return at(x);
    }
 
    virtual ~GridBase() {}
};
 
} // End namespace SSAGES