gsGridIterator.h

Go to the documentation of this file.

#pragma once

namespace gismo
{

enum gsGridIteratorMode
{
    CUBE   = 0, 
    BDR    = 1, 
    VERTEX = 2, 
    CWISE  = 3  
};
    
// note: default arguments are found in gsForwardDeclarations.h
template<typename Z, int mode, short_t d, bool> class gsGridIterator
{using Z::GISMO_ERROR_gsGridIterator_has_invalid_template_arguments;};

template<class Z, int mode, short_t d>
class gsGridIterator<Z,mode,d,true>
{
public:
    typedef gsVector<Z,d> point;

public:

    gsGridIterator()
    {
        GISMO_STATIC_ASSERT(std::numeric_limits<Z>::is_integer,"The template parameter needs to be an integer type.");
        GISMO_STATIC_ASSERT(mode > -1 && mode < 3, "The mode of gsGridIterator needs to be 0, 1 or 2.");
    }

    gsGridIterator(point const & a, point const & b, bool open = true)
    { reset(a, b, open); }

    explicit gsGridIterator(point const & b, bool open = true)
    { reset(point::Zero(b.size()), b, open); }

    explicit gsGridIterator(gsMatrix<Z,d,2> const & ab, bool open = true)
    { reset(ab.col(0), ab.col(1), open); }

    inline void reset(point const & a, point const & b, bool open = true)
    {
        GISMO_ASSERT(a.rows() == b.rows(), "Invalid endpoint dimensions");
        m_low = m_cur = a;
        if (open) m_upp = b.array() - 1; else m_upp = b;
        m_valid = ( (m_low.array() <= m_upp.array()).all() ? true : false );
    }

    void reset() { reset(m_low,m_upp, false); }

    // See http://eigen.tuxfamily.org/dox-devel/group__TopicStructHavingEigenMembers.html
    //EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF( (sizeof(point)%16)==0 );

public:

    operator bool() const {return m_valid;}
    
    const point & operator*() const {return m_cur;}

    const point * operator->() const {return &m_cur;}
    
    inline gsGridIterator & operator++()
    {
        // Multiple implementations according to mode
        switch (mode)
        {
        case 0: // ----------- Iteration over [m_low, m_upp]
        case 2: // iteration over vertices of the cube [m_low, m_upp]
            for (index_t i = 0; i != m_cur.size(); ++i)
                if ( m_cur[i] != m_upp[i] )
                {
                    if (0==mode) ++m_cur[i]; else m_cur[i] = m_upp[i];
                    return *this;
                }
                else
                    m_cur[i] = m_low[i];
            m_valid = 0;//done
            return *this;

        case 1: // ----------- Iteration over boundary of [m_low, m_upp]
        default:
            for (index_t i = 0; i != m_cur.size(); ++i)
            {        
                if ( m_cur[i] != m_upp[i] )
                {
                    if ( m_cur[i] == m_low[i] && (i+1!=m_cur.size() || 1==m_cur.size()) )
                    {
                        index_t c = i+1;
                        for (index_t j = c; j!=m_cur.size(); ++j)
                            if ( (m_cur[j] == m_low[j]) || 
                                 (m_cur[j] == m_upp[j]) )
                                ++c;
                        
                        if ( c==1 )
                            m_cur[i] = m_upp[i];
                        else
                            ++m_cur[i];
                    }
                    else
                        m_cur[i]++;

                    m_cur.head(i) = m_low.head(i);
                    return *this;
                }
            }
            /*fall through to default*/
            m_valid = 0;//done
            return *this;
        }
    }

    inline bool isFloor(int i) const { return m_cur[i] == m_low[i];}

    inline bool isCeil (int i) const
    { return mode ? m_cur[i] == m_upp[i] : m_cur[i] + 1 == m_upp[i];}

    bool isBoundary() const
    {
        return (m_cur.array() == m_low.array()).any() ||
            (m_cur.array() == m_upp.array()).any() ;
    }

    index_t numPoints() const
    {
        switch (mode)
        {
        case 0: return ((m_upp-m_low).array() + 1).prod();
        case 1: return ((m_upp-m_low).array() + 1).prod() - ((m_upp-m_low).array()-1).prod();
        case 2: return ( 1 << m_low.rows() );
        default: GISMO_ERROR("gsGridIterator::numPoints");
        }
    }

    point numPointsCwise() const { return (m_upp-m_low).array() + 1;}

    const point & lower() const {return m_low;}

    const point & upper() const {return m_upp;}

    point strides() const
    {
        point res(m_low.rows());
        res[0] = 1;
        for (index_t i = 1; i != res.size(); ++i )
            res[i] = res[i-1] * ( m_upp[i-1] - m_low[i-1] + 1 );
        return res;
    }

private:
    
    point m_low, m_upp; 
    point m_cur;        
    bool  m_valid;      
};


template<class T, int mode, short_t d>
class gsGridIterator<T,mode,d,false>
{
public:

    typedef gsVector<T,d> point;

    typedef gsGridIterator<index_t, mode, d> integer_iterator;

    typedef typename integer_iterator::point point_index;
public:

    gsGridIterator(point const & a, 
                   point const & b, 
                   point_index const & np)
    : m_iter(np, 1)
    {
        reset(a, b);
        GISMO_STATIC_ASSERT(mode > -1 && mode < 3, "The mode of gsGridIterator needs to be 0, 1 or 2.");
    }

    gsGridIterator(gsMatrix<T,d,2> const & ab, 
                   point_index const & np)
    : m_iter(np, 1)
    {
        reset(ab.col(0), ab.col(1));
    }

    gsGridIterator(gsMatrix<T,d,2> const & ab, unsigned numPoints)
    : m_low(ab.col(0)), m_upp(ab.col(1))
    {
        // deduce the number of points per direction for an approximately uniform grid
        const gsVector<double,d> L = (ab.col(1) - ab.col(0)).template cast<double>();
        const double h = std::pow(L.prod()/numPoints, 1.0 / m_low.rows());
        const point_index npts = (L/h).unaryExpr((double(*)(double))std::ceil).template cast<index_t>();
        m_iter = integer_iterator(npts, 1);
        reset(ab.col(0), ab.col(1));
    }

    void reset() { m_cur = m_low; m_iter.reset(); }

    void reset(point const & a, 
               point const & b)
    {
        m_cur = m_low = a;
        m_upp = b;
        m_step = (b-a).array() / (m_iter.numPointsCwise().array() - 1)
            .matrix().cwiseMax(1).template cast<T>().array() ;
        m_iter.reset();
    }

    gsMatrix<T> toMatrix() const
    {
        gsMatrix<T> res(m_cur.rows(), numPoints());
        integer_iterator iter = m_iter;
        iter.reset();

        for(index_t c = 0; iter; ++iter, ++c)
            update(*iter,res.col(c).data());
        
        return res;
    }
    
    // See http://eigen.tuxfamily.org/dox-devel/group__TopicStructHavingEigenMembers.html
    //EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF( (sizeof(point)%16)==0 );

public:

    operator bool() const {return m_iter;}
        
    const gsMatrix<T> & operator*() const {return m_cur;}

    const gsMatrix<T> * operator->() const {return &m_cur;}

    inline gsGridIterator & operator++()
    {
        if ( ++m_iter ) update(*m_iter, m_cur.data());
        return *this;        
    }

    inline bool isFloor(int i) const { return m_iter.isFloor(i);}

    inline bool isCeil (int i) const { return m_iter.isCeil(i);}

    inline bool isBoundary() const { return m_iter.isBoundary();}

    index_t numPoints() const {return m_iter.numPoints();}
    
    point_index numPointsCwise() const {return m_iter.numPointsCwise();}

    const point & lower() const {return m_low;}

    const point & upper() const {return m_upp;}

    const point_index & tensorIndex() const { return *m_iter;}

    index_t index(const index_t i) const { return m_iter->at(i);}

    const point & step() const {return m_step;}

    const integer_iterator & index_iterator() const { return m_iter;}

    inline const gsMatrix<T> operator [] (const point_index & ti) const
    {
        gsMatrix<T> res(m_low.rows(),1);
        update(ti, res.data());
        return res;
    }

private:
    
    // Update the point \a pt to the position \a ti
    inline void update(const point_index & ti, T * pt) const
    {
        for (index_t i = 0; i != m_low.size(); ++i)
            if ( ti[i] == m_iter.lower()[i] )
                *(pt++) = m_low[i]; // avoid numerical error at first val
            else if ( ti[i] == m_iter.upper()[i] )
                *(pt++) = m_upp[i]; // avoid numerical error at last val
            else
                *(pt++) = m_low[i] + static_cast<T>(ti[i]) * m_step[i];
    }

private:

    point  m_low, m_upp, m_step; 
    
    integer_iterator m_iter;     
    gsMatrix<T>      m_cur;      
};


template<class T, short_t d> // mode == 3 == CWISE
class gsGridIterator<T,CWISE,d,false> 
{
public:

    typedef gsGridIterator<index_t, 0, d> integer_iterator;

    typedef typename integer_iterator::point point_index;

    typedef gsVector<const T *,d, 2> CwiseData; //non-aligned array
    
public:

    template<class CwiseContainer>
    explicit gsGridIterator(const CwiseContainer & cwise)
    : m_cwise(cwise.size()), m_cur(m_cwise.size(),1)
    {
        point_index npts(m_cwise.size());
        for (index_t i = 0; i != npts.size(); ++i)
        {
            m_cwise[i] = cwise[i].data();
            npts[i]    = cwise[i].size() - 1;
            GISMO_ASSERT(cwise[i].cols()==1 || cwise[i].rows()==1, "Invalid input");
        }
        m_iter = integer_iterator(npts, 0);
        //if (1==cwise.front().cols())
        //    m_cur.derived().resize(1, npts.size());
        //else
        //m_cur.derived().resize(npts.size(), 1);
        update(*m_iter, m_cur.data());
    }

    template<class Matrix_t>
    explicit gsGridIterator(const std::vector<Matrix_t*> & cwise)
    : m_cwise(cwise.size()), m_cur(m_cwise.size(),1)
    {
        point_index npts(m_cwise.size());
        for (index_t i = 0; i != npts.size(); ++i)
        {
            m_cwise[i] = cwise[i]->data();
            npts[i]    = cwise[i]->size() - 1;
            GISMO_ASSERT(cwise[i]->cols()==1 || cwise[i]->rows()==1, "Invalid input");
        }
        m_iter = integer_iterator(npts, 0);
        update(*m_iter, m_cur.data());
    }

    void reset() { m_iter.reset(); update(*m_iter, m_cur.data());}

    //void restart() { m_iter.reset(); update(*m_iter, m_cur.data());}

public:

    operator bool() const {return m_iter;}
        
    const gsMatrix<T> & operator*() const {return m_cur;}

    const gsMatrix<T> * operator->() const {return &m_cur;}
    
    inline gsGridIterator & operator++()
    {
        if (++m_iter) update(*m_iter, m_cur.data());
        return *this;
    }

    inline bool isFloor(int i) const { return m_iter.isFloor(i);}

    inline bool isCeil (int i) const { return m_iter.isCeil(i);}

    inline bool isBoundary() const { return m_iter.isBoundary();}

    index_t numPoints() const {return m_iter.numPoints();}
    
    point_index numPointsCwise() const {return m_iter.numPointsCwise();}

    const point_index & tensorIndex() const { return *m_iter;}

    index_t index(const index_t i) const { return m_iter->at(i);}

    const integer_iterator & index_iterator() const { return m_iter;}

    inline const gsMatrix<T> operator [] (const point_index & ti) const
    {
        gsMatrix<T> res(m_cur.rows(),1);
        update(ti, res.data());
        return res;
    }

private:

    // Update the point \a pt to the position \a ti
    inline void update(const point_index & ti, T * pt)
    {
        for (index_t i = 0; i != m_cur.rows(); ++i)
            *(pt++) = m_cwise[i][ti[i]];
    }

private:

    CwiseData m_cwise; 
    
    integer_iterator m_iter; 
    gsMatrix<T>      m_cur;  
};


} // namespace gismo