gsHDomainIterator.h

Go to the documentation of this file.

#pragma once

#include <gsHSplines/gsHDomain.h>
#include <gsHSplines/gsKdNode.h>
#include <gsNurbs/gsTensorBSplineBasis.h>

#include <gsCore/gsDomainIterator.h>

namespace gismo
{

// Documentation in gsDomainIterator
template<typename T, unsigned d>
class gsHDomainIterator: public gsDomainIterator<T>
{
public:

    typedef gsKdNode<d, index_t> node;

    typedef typename node::point point;

    typedef typename std::vector<T>::const_iterator  uiter;

    typedef gsHDomain<d,index_t> hDomain;

    typedef typename hDomain::const_literator leafIterator;

public:

    gsHDomainIterator(const gsHTensorBasis<d, T> & hbs)
    : gsDomainIterator<T>(hbs)
    {
        // Initialize mesh data
        m_meshStart.resize(d);
        m_meshEnd  .resize(d);

        // Initialize cell data
        m_curElement.resize(d);
        m_lower     .resize(d);
        m_upper     .resize(d);

        // Allocate breaks
        m_breaks = std::vector<std::vector<T> >(d, std::vector<T>());

        m_leaf = hbs.tree().beginLeafIterator();
        updateLeaf();
        updateElement();
    }

    // ---> Documentation in gsDomainIterator.h
    bool next()
    {
        this->m_isGood = nextLexicographic(m_curElement, m_meshStart, m_meshEnd);

        if (this->m_isGood) // new element in m_leaf
            updateElement();
        else // went through all elements in m_leaf
        {
            this->m_isGood = nextLeaf();
            if (this->m_isGood)
                updateElement();
        }

        ++m_id; //increment id
        return this->m_isGood;
    }

    // ---> Documentation in gsDomainIterator.h
    bool next(index_t increment)
    {
        for (index_t i = 0; i != increment && this->m_isGood; ++i)
        {
            this->m_isGood = nextLexicographic(m_curElement, m_meshStart, m_meshEnd);
            if (!this->m_isGood)
                this->m_isGood = nextLeaf();
        }

        if (this->m_isGood)
            updateElement();

        m_id += increment; //increment id
        return this->m_isGood;
    }

    void reset()
    {
        m_id = 0;
        const gsHTensorBasis<d, T>* hbs =  dynamic_cast<const gsHTensorBasis<d, T> *>(m_basis);
        m_leaf = hbs->tree().beginLeafIterator();
        updateLeaf();
        updateElement();
    }

    const gsVector<T>& lowerCorner() const { return m_lower; }

    const gsVector<T>& upperCorner() const { return m_upper; }

    int getLevel() const
    {
        return m_leaf.level();
    }

    // Returns the element multi-index at the current level
    // If you need the element at the level above, divide this all indices by 2
    gsVector<index_t> elementMultiIndex() const
    {
        gsVector<index_t> res(d);
        for (index_t i = 0; i!=d; ++i)
        {
            res[i] =  std::distance(m_breaks[i].begin(), m_curElement[i]);
        }
        return res;
    }

private:

    gsHDomainIterator();

    bool nextLeaf()
    {
        this->m_isGood = m_leaf.next();

        if ( m_leaf.good() )
            updateLeaf();

        return this->m_isGood;
    }

    void updateLeaf()
    {
        const point & lower = m_leaf.lowerCorner();
        const point & upper = m_leaf.upperCorner();
        // gsDebug<<"leaf "<<  lower.transpose() <<", "
        //        << upper.transpose() <<"\n";

        const int level2 = m_leaf.level();

        // Update leaf box
        for (unsigned dim = 0; dim < d; ++dim)
        {
            index_t start = lower(dim);
            index_t end  = upper(dim) ;

            if (basis().manualLevels() )
            {
                static_cast<const gsHTensorBasis<d,T>*>(m_basis)->
                    _diadicIndexToKnotIndex(level2,dim,start);
                static_cast<const gsHTensorBasis<d,T>*>(m_basis)->
                    _diadicIndexToKnotIndex(level2,dim,end);
            }

            const gsKnotVector<T> & kv =
                static_cast<const gsHTensorBasis<d,T>*>(m_basis)
                ->tensorLevel(level2).component(dim).knots();

            // knotVals = kv.unique()

            m_breaks[dim].clear();
            for (index_t index = start; index <= end; ++index)
                m_breaks[dim].push_back( kv(index) );// unique index

            m_curElement(dim) =
            m_meshStart(dim)  = m_breaks[dim].begin();

            // for n breaks, we have n - 1 elements (spans)
            m_meshEnd(dim) =  m_breaks[dim].end() - 1;
        }
    }

    void updateElement()
    {
        // Update cell data
        for (unsigned i = 0; i < d ; ++i)
        {
            m_lower[i]  = *m_curElement[i];
            m_upper[i]  = *(m_curElement[i]+1);
            center[i] = (T)(0.5) * (m_lower[i] + m_upper[i]);
        }
    }

// =============================================================================
// members
// =============================================================================

    const gsHTensorBasis<d,T> & basis() const { return *static_cast<const gsHTensorBasis<d,T>*>(m_basis); }

public:

    using gsDomainIterator<T>::center;
    using gsDomainIterator<T>::m_basis;

#   define Eigen gsEigen
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
#   undef Eigen

protected:
    using gsDomainIterator<T>::m_id;

private:

    // The current leaf node of the tree
    leafIterator m_leaf;

    // Coordinates of the grid cell boundaries
    // \todo remove this member
    std::vector< std::vector<T> > m_breaks;

    // Extent of the tensor grid
    gsVector<uiter, d> m_meshStart, m_meshEnd;

    // Current element as pointers to it's supporting mesh-lines
    gsVector<uiter, d> m_curElement;

    // parameter coordinates of current grid cell
    gsVector<T> m_lower, m_upper;
};

} // end namespace gismo