grid_sm.hpp

#ifndef GRID_HPP
#define GRID_HPP

#include "config.h"
#include <boost/shared_array.hpp>
#include <vector>
#include <initializer_list>
#include <array>
#include "memory/memory.hpp"
#include "Space/Shape/Box.hpp"
#include "grid_key.hpp"
#include <iostream>
#include "util/mathutil.hpp"
#include "iterators/stencil_type.hpp"

#define PERIODIC 1
#define NON_PERIODIC 0

// Box need the definition of grid_key_dx_r
#define HARDWARE 1

class NoCheck
{
public:
    static bool valid(size_t v_id, size_t sz)
    {
        return true;
    }
};

class CheckExistence
{
public:
    static bool valid(size_t v_id, size_t sz)
    {
        return v_id < sz;
    }
};

template<unsigned int N, typename T> class grid_sm;

template <unsigned int dim> class print_warning_on_adjustment;

template<unsigned int dim,typename stencil=no_stencil,typename warn=print_warning_on_adjustment<dim>> class grid_key_dx_iterator_sub;

template<unsigned int N, typename T>
class grid_sm
{
    Box<N,size_t> box;

    size_t size_tot;

    size_t sz[N];

    size_t sz_s[N];

    inline void Initialize(const size_t sz)
    {
        sz_s[0] = sz;
        this->sz[0] = sz;

        // set the box
        box.setHigh(0,sz);
        box.setLow(0,0);

        for (size_t i = 1 ;  i < N ; i++)
        {
            /* coverity[dead_error_begin] */
            sz_s[i] = sz*sz_s[i-1];
            this->sz[i] = sz;

            // set the box
            box.setHigh(i,sz);
            box.setLow(i,0);
        }
    }

    inline void Initialize(const size_t (& sz)[N])
    {
        sz_s[0] = sz[0];
        this->sz[0] = sz[0];

        // set the box
        box.setHigh(0,sz[0]);
        box.setLow(0,0);

        for (size_t i = 1 ;  i < N ; i++)
        {
            /* coverity[dead_error_begin] */
            sz_s[i] = sz[i]*sz_s[i-1];
            this->sz[i] = sz[i];

            // set the box
            box.setHigh(i,sz[i]);
            box.setLow(i,0);
        }
    }

    inline void Initialize()
    {
        sz_s[0] = 0;
        this->sz[0] = 0;

        // set the box
        box.setHigh(0,0);
        box.setLow(0,0);

        for (size_t i = 1 ;  i < N ; i++)
        {
            /* coverity[dead_error_begin] */
            sz_s[i] = sz[i]*sz_s[i-1];

            // set the box
            box.setHigh(i,sz[i]);
            box.setLow(i,0);
        }
    }

public:

    inline const Box<N,size_t> & getBox() const
    {
        return box;
    }

    inline const Box<N,size_t> getBoxKey() const
    {
        Box<N,size_t> bx;

        for (size_t i = 0 ; i < N ; i++)
        {
            bx.setLow(i,box.getLow(i));
            bx.setHigh(i,box.getHigh(i) - 1);
        }

        return bx;
    }

    inline void setDimensions(const size_t  (& dims)[N])
    {
        Initialize(dims);
        size_tot = totalSize(dims);
    }

    inline grid_sm()
    :size_tot(0)
    {
        // Initialize sz
        for (size_t i = 0 ; i < N ; i++)
        {sz[i] = 0;}

        Initialize();
    }

    template<typename S> inline grid_sm(const grid_sm<N,S> & g)
    {
        box = g.box;
        size_tot = g.size_tot;

        for (size_t i = 0 ; i < N ; i++)
        {
            sz[i] = g.sz[i];
            sz_s[i] = g.sz_s[i];
        }
    }

    // Static element to calculate total size

    inline size_t totalSize(const size_t sz)
    {
        size_t tSz = 1;

        for (size_t i = 0 ;  i < N ; i++)
        {
            tSz *= sz;
        }

        return tSz;
    }

    // Static element to calculate total size

    inline size_t totalSize(const size_t (& sz)[N])
    {
        size_t tSz = 1;

        for (size_t i = 0 ;  i < N ; i++)
        {
            tSz *= sz[i];
        }

        return tSz;
    }


    inline grid_sm(const size_t & sz)
    : size_tot(totalSize(sz))
    {
        Initialize(sz);
    }

    inline grid_sm(const size_t (& sz)[N])
    : size_tot(totalSize(sz))
    {
        Initialize(sz);
    }

    template<typename check=NoCheck> inline mem_id LinId(const grid_key_dx<N> & gk, const char sum_id[N]) const
    {
        mem_id lid;

        // Check the sum produce a valid key

        if (check::valid(gk.get(0) + sum_id[0],sz[0]) == false)
            return -1;

        lid = gk.get(0) + sum_id[0];


        for (mem_id i = 1 ; i < N ; i++)
        {
            // Check the sum produce a valid key

            if (check::valid(gk.get(i) + sum_id[i],sz[i]) == false)
                return -1;

            lid += (gk.get(i) + sum_id[i]) * sz_s[i-1];
        }

        return lid;
    }

    template<typename check=NoCheck> inline mem_id LinId(const grid_key_dx<N> & gk, const char sum_id[N], const size_t (&bc)[N]) const
    {
        mem_id lid;

        // Check the sum produce a valid key

        if (bc[0] == NON_PERIODIC)
        {
            if (check::valid(gk.get(0) + sum_id[0],sz[0]) == false)
                return -1;

            lid = gk.get(0) + sum_id[0];
        }
        else
        {
            lid = openfpm::math::positive_modulo(gk.get(0) + sum_id[0],sz[0]);
        }

        for (mem_id i = 1 ; i < N ; i++)
        {
            // Check the sum produce a valid key

            /* coverity[dead_error_line] */
            if (bc[i] == NON_PERIODIC)
            {
                if (check::valid(gk.get(i) + sum_id[i],sz[i]) == false)
                    return -1;

                lid += (gk.get(i) + sum_id[i]) * sz_s[i-1];
            }
            else
            {
                lid += (openfpm::math::positive_modulo(gk.get(i) + sum_id[i],sz[i])) * sz_s[i-1];
            }
        }

        return lid;
    }

    inline mem_id LinIdPtr(size_t * k) const
    {
        mem_id lid = k[0];
        for (mem_id i = 1 ; i < N ; i++)
        {
            lid += k[i] * sz_s[i-1];
        }

        return lid;
    }

    inline mem_id LinId(const size_t (& k)[N]) const
    {
        mem_id lid = k[0];
        for (mem_id i = 1 ; i < N ; i++)
        {
            /* coverity[dead_error_line] */
            lid += k[i] * sz_s[i-1];
        }

        return lid;
    }

    inline mem_id LinId(const grid_key_dx<N> & gk) const
    {
        mem_id lid = gk.get(0);
        for (mem_id i = 1 ; i < N ; i++)
        {
            /* coverity[dead_error_begin */
            lid += gk.get(i) * sz_s[i-1];
        }

        return lid;
    }

    template<typename a, typename ...lT> inline mem_id Lin(a v,lT...t) const
    {
#ifdef DEBUG
        if (sizeof...(t)+1 > N)
        {
            std::cerr << "Error incorrect grid cannot linearize more index than its dimensionality" << "\n";
        }
#endif

        return v*sz_s[sizeof...(t)-1] + Lin(t...);
    }

    template<typename a> inline mem_id Lin(a v) const
    {
        return v;
    }


    inline grid_key_dx<N> InvLinId(mem_id id) const
    {
        // Inversion of linearize

        grid_key_dx<N> gk;

        for (mem_id i = 0 ; i < N ; i++)
        {
            gk.set_d(i,id % sz[i]);
            id /= sz[i];
        }

        return gk;
    }

    //#pragma openfpm layout(get)
    template<unsigned int dim, unsigned int p> inline mem_id LinId(const grid_key_d<dim,p> & gk) const
    {
        mem_id lid = gk.k[0];
        for (mem_id i = 1 ; i < dim ; i++)
        {
            lid += gk.k[i] * sz_s[i-1];
        }

        return lid;
    }

    //#pragma openfpm layout(get)
    inline mem_id LinId(mem_id * id) const
    {
        mem_id lid = 0;
        lid += id[0];
        for (mem_id i = 1 ; i < N ; i++)
        {
            lid += id[i] * sz_s[i-1];
        }

        return lid;
    }

    ~grid_sm() {};

    inline size_t size() const
    {
        return size_tot;
    };

    inline grid_sm<N,T> & operator=(const grid_sm<N,T> & g)
    {
        box = g.box;
        size_tot = g.size_tot;

        for (size_t i = 0 ; i < N ; i++)
        {
            sz[i] = g.sz[i];
            sz_s[i] = g.sz_s[i];
        }

        return *this;
    }

    inline bool operator==(const grid_sm<N,T> & g)
    {
        for (size_t i = 0 ; i < N ; i++)
        {
            if (sz[i] != g.sz[i])
                return false;
        }

        if (box != g.box)
            return false;

#ifdef SE_CLASS1

        if (size_tot != g.size_tot)
            return false;

        for (size_t i = 0 ; i < N ; i++)
        {
            if (sz_s[i] != g.sz_s[i])
                return false;
        }

#endif
        return true;
    }

    inline bool operator!=(const grid_sm<N,T> & g)
    {
        return ! this->operator==(g);
    }

    inline size_t size_s(unsigned int i) const
    {
        return sz_s[i];
    }

    inline size_t size(unsigned int i) const
    {
        return sz[i];
    }

    inline const size_t (& getSize() const)[N]
    {
        return sz;
    }

    inline grid_key_dx_iterator_sub<N> getSubIterator(grid_key_dx<N> & start, grid_key_dx<N> & stop) const
    {
        return grid_key_dx_iterator_sub<N>(*this,start,stop);
    }

    inline void swap(grid_sm<N,T> & g)
    {
        Box<N,size_t> box_t = box;
        box = g.box;
        g.box = box_t;

        size_t tmp = size_tot;
        size_tot = g.size_tot;
        g.size_tot = tmp;

        for (size_t i = 0 ; i < N ; i++)
        {
            tmp = sz[i];
            sz[i] = g.sz[i];
            g.sz[i] = tmp;

            tmp = sz_s[i];
            sz_s[i] = g.sz_s[i];
            g.sz_s[i] = tmp;
        }
    }

    std::string toString() const
    {
        std::stringstream str;

        for (size_t i = 0 ; i < N ; i++)
            str << "sz[" << i << "]=" << size(i) << " ";

        return str.str();
    }

    template <unsigned int,typename> friend class grid_sm;
};




class grid_key_dx_r
{
    size_t dim;

public:

    size_t getDim()
    {
        return dim;
    }

    grid_key_dx_r(grid_key_dx_r & key)
    :dim(key.dim)
    {
        // Allocate the key
        k = new mem_id[dim];

        // Copy the key
        for(unsigned int i = 0 ; i < dim ; i++)
        {
            k[i] = key.k[i];
        }
    }

    grid_key_dx_r(size_t dim)
    :dim(dim)
    {
        // Allocate the key
        k = new mem_id[dim];
    }

    ~grid_key_dx_r()
    {
        delete [] k;
    }

    template<typename a, typename ...T>void set(a v, T...t)
    {
        k[dim-1] = v;
        invert_assign(t...);
    }

    mem_id get(size_t i)
    {
        return k[i];
    }

    void set_d(size_t i, mem_id id)
    {
        k[i] = id;
    }

    mem_id * k;

private:

    template<typename a, typename ...T>void invert_assign(a v,T...t)
    {
        k[sizeof...(T)] = v;
        invert_assign(t...);
    }

    template<typename a, typename ...T>void invert_assign(a v)
    {
        k[0] = v;
    }

    void invert_assign()
    {
    }

};


class Iterator_g_const
{
    size_t dim;

    size_t sz;

    // Actual grid_key position
    grid_key_dx_r gk;

public:

    size_t getDim()
    {
        return dim;
    }

    Iterator_g_const(size_t n, size_t sz)
    :dim(n),sz(sz),gk(n)
    {
        // fill gk with the first grid element that satisfied the constrain: 0,1,2,3... dim

        for (size_t i = 0 ; i < dim ; i++)
        {
            gk.set_d(i,dim-i-1);
        }
    }

    Iterator_g_const & operator++()
            {

        gk.set_d(0,gk.get(0)+1);


        unsigned int i = 0;
        for ( ; i < dim-1 ; i++)
        {
            size_t id = gk.get(i);
            if (id >= sz)
            {
                // ! overflow, increment the next index

                id = gk.get(i+1);
                if (id+i+2 >= sz)
                {
                    // there is no-way to produce a valid key
                    // there is not need to check the previous index
                    // overflow i+1

                    gk.set_d(i+1,sz);
                }
                else
                {

                    // reinitialize the previous index

                    for (unsigned int s = 0 ; s <= i+1 ; s++)
                    {
                        gk.set_d(i+1-s,id+1+s);
                    }
                }
            }
            else
            {
                break;
            }
        }

        return *this;
            }

    bool isNext()
    {
        // If dimensionless return immediately
        if (dim == 0)
            return false;

        if (gk.get(dim-1) < static_cast<mem_id>(sz-dim+1))
        {

            return true;
        }

        return false;
    }

    grid_key_dx_r & get()
    {
        return gk;
    }
};

#endif