gsHDomainBoundaryIterator.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 gsHDomainBoundaryIterator: 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:

    gsHDomainBoundaryIterator(const gsHTensorBasis<d, T> & hbs, 
                              const boxSide & s )
        : 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>());

        // Set to one quadrature point by default
        m_quadrature.setNodes( gsVector<index_t>::Ones(d) );

        // Get the side information
        par = s.parameter();
        dir = s.direction();

        initLeaf(hbs.tree());
    }

    // ---> 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();

        return this->m_isGood;
    }

    // ---> Documentation in gsDomainIterator.h
    bool next(index_t increment)
    {
        for (index_t i = 0; i < increment; i++)
            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();

        return this->m_isGood;
    }
    
    void reset()
    {
        const gsHTensorBasis<d, T>* hbs =  dynamic_cast<const gsHTensorBasis<d, T> *>(m_basis);
        initLeaf(hbs->tree());
    }

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

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

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

private:

    gsHDomainBoundaryIterator();

    void initLeaf(const hDomain & tree_domain)
    {
        // Get the first leaf
        m_leaf = tree_domain.beginLeafIterator();

        for (; m_leaf.good(); m_leaf.next() )
        {
            // Check if this leaf is on our side
            if ( leafOnBoundary() )
            {
                updateLeaf();
                return;
            }
        }
        GISMO_ERROR("No leaves.\n");
    }


    bool nextLeaf()
    {
        for (m_leaf.next(); m_leaf.good(); m_leaf.next() )
        {
            // Check if this leaf is on our side
            if ( leafOnBoundary() )
            {
                updateLeaf();
                return true;
            }
        }
        return false;
    }

    bool leafOnBoundary() const
    {
        if ( par )
        {
            // AM: a little ugly for now, to be improved
            size_t diadicSize;
            const gsHTensorBasis<d,T> * hbasis = dynamic_cast<const gsHTensorBasis<d,T> * >(m_basis);
            if (basis().manualLevels() )
            {
                gsKnotVector<T> kv = hbasis->tensorLevel(m_leaf.level()).knots(dir);
                index_t start = 0;
                index_t end  = kv.uSize()-1;
                hbasis->_knotIndexToDiadicIndex(m_leaf.level(),dir,start);
                hbasis->_knotIndexToDiadicIndex(m_leaf.level(),dir,end);
                diadicSize = end - start;
            }
            else
                diadicSize = static_cast<const gsHTensorBasis<d,T>*>(m_basis)->tensorLevel(m_leaf.level()).knots(dir).uSize() - 1;

            return static_cast<size_t>(m_leaf.upperCorner().at(dir) ) == diadicSize;// todo: more efficient
        }
        else
        {
            return m_leaf.lowerCorner().at(dir) == 0;
        }
    }

    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();

            m_breaks[dim].clear();
            if ( dim == dir )
            {
                if ( par )
                {
                    m_breaks[dim].push_back( kv(end-1) );
                    m_breaks[dim].push_back( kv(end  ) );
                }
                else
                {
                    m_breaks[dim].push_back( kv(start)   );
                    m_breaks[dim].push_back( kv(start+1) );
                }
            }
            else
            {
                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;
        }

        // We are at a new element, so update cell data
        updateElement();
    }

    void updateElement()
    {
        // Update cell data
        for (unsigned i = 0; i < dir ; ++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]);
        }
        m_lower[dir] = 
        m_upper[dir] =
        center [dir] = (par ? *(m_curElement[dir]+1) : *m_curElement[dir] );
        for (unsigned i = dir+1; 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

private:

    // Boundary parameters
    unsigned dir; // direction normal to the boundary
    bool par;     // parameter value

    // 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;

    // Quadrature rule
    gsGaussRule<T> m_quadrature;

};

} // end namespace gismo