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