gsHTensorBasis.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsBasis.h>

#include <gsHSplines/gsHDomain.h>
#include <gsHSplines/gsHDomainIterator.h>
#include <gsHSplines/gsHDomainBoundaryIterator.h>

#include <gsCore/gsBoundary.h>

#include <gsNurbs/gsTensorBSplineBasis.h>
#include <gsNurbs/gsBSplineBasis.h> // for gsBasis::component(short_t)

#include <gsUtils/gsSortedVector.h>

namespace gismo
{

struct lvl_coef
{
    int pos; // flat index at grid of level \a lvl
    int unsigned lvl; // level
    real_t coef; // value of the coefficient (lvl,pos)
};

template<short_t d, class T>
class GISMO_DEFAULT_VIS gsHTensorBasis: public gsBasis<T>
{
public:
    typedef memory::shared_ptr< gsHTensorBasis > Ptr;

    typedef memory::unique_ptr< gsHTensorBasis > uPtr;

    typedef gsHTensorBasis<d,T> Self_t;

    typedef T Scalar_t;

    typedef gsHDomain<d> hdomain_type;

    typedef typename hdomain_type::point point;

    typedef typename hdomain_type::box   box;

    typedef std::vector< box > boxHistory;

    typedef gsSortedVector< index_t > CMatrix; // charMatrix_

    typedef typename CMatrix::const_iterator cmatIterator;

    typedef gsKnotVector<T> KnotVectorType;

    typedef typename gsBSplineTraits<static_cast<short_t>(d-1),T>::Basis BoundaryBasisType;

    typedef typename gsBSplineTraits<d,T>::Basis tensorBasis;

    static const short_t Dim = d;

    friend class gsHDomainIterator<T,d>;
    friend class gsHDomainBoundaryIterator<T,d>;
public:

    gsHTensorBasis()
    :
    m_manualLevels(false)
    {
        initialize_class(tensorBasis());
        update_structure();
    }

    gsHTensorBasis( gsBasis<T> const&  tbasis, bool manualLevels)
    :
    m_manualLevels(manualLevels)
    {
        if (!m_manualLevels) 
        {
            initialize_class(tbasis);
            // Build the characteristic matrices
        }
        else
        {
            // Degrees
            m_deg.resize(d);
            for( index_t i = 0; i < d; i++)
                m_deg[i] = tbasis.degree(i);

            // Construct the initial basis
            m_bases.reserve(3);
            if ( const tensorBasis * tb2 = dynamic_cast<const tensorBasis*>(&tbasis) )
            {
                m_bases.push_back(tb2->clone().release());
                // For the tbasis, count the unique knot values
                std::vector<std::vector<index_t>> lvlIndices(d);
                std::vector<index_t> dirIndices;
                for (short_t dim=0; dim!=d; dim++)
                {
                    dirIndices.resize(tb2->knots(dim).uSize());
                    std::iota(dirIndices.begin(),dirIndices.end(),0);
                    lvlIndices[dim] = dirIndices;
                }
                m_uIndices.push_back(lvlIndices);
            }
            else
            {
                GISMO_ERROR("Cannot construct a Hierarchical basis from "<< tbasis );
            }

            // Initialize the binary tree
            point upp;
            for ( index_t i = 0; i!=d; ++i )
                upp[i] = m_bases[0]->knots(i).numElements();

            m_tree.init(upp);
        }
        update_structure();
    }

    gsHTensorBasis( gsTensorBSplineBasis<d,T> const&  tbasis,
                    const std::vector<index_t> & boxes)
    :
    m_manualLevels(false)
    {
        initialize_class(tbasis);
        point i1;
        point i2;
        //i1.resize(d);
        //i2.resize(d);
        // Set all functions to active
        GISMO_ASSERT( (boxes.size()%(2*d+1))==0,
                      "The points did not define boxes properly. The basis was created without any domain structure.");

        for( size_t i = 0; i < (boxes.size()/(2*d+1)); i++)
        {
            for( short_t j = 0; j < d; j++)
            {
                i1[j] = boxes[(2*d+1)*i+j+1];
                i2[j] = boxes[(2*d+1)*i+j+d+1];
            }
            insert_box(i1,i2,boxes[i*(2*d+1)]);
        }

        // Build the characteristic matrices (note: call is non-vritual)
        update_structure();
    }

    gsHTensorBasis( gsTensorBSplineBasis<d,T> const&  tbasis,
                    gsMatrix<T> const & boxes)
    :
    m_manualLevels(false)
    {
        //assert(boxes.rows() == 2);    //can accept only 2D coordinates- Delete during the generalization of the lac to nD
        GISMO_ASSERT(boxes.rows() == d, "Points in boxes need to be of dimension d.");
        GISMO_ASSERT(boxes.cols()%2 == 0, "Each box needs two corners but you don't provied gsHTensorBasis constructor with them.");
        initialize_class(tbasis);

        point k1;
        point k2;

        for(index_t i = 0; i < boxes.cols()/2; i++)
        {
            for(short_t j = 0; j < d; j++)
            {
                k1[j] = this->m_bases.back()->knots(j).uFind(boxes(j,2*i)).uIndex();
                k2[j] = this->m_bases.back()->knots(j).uFind(boxes(j,2*i+1)).uIndex()+1;
            }
            int level = m_tree.query3(k1,k2,m_bases.size()-1);
            for(short_t j = 0; j < d; j++)
            {
                k1[j] = this->m_bases[level+1]->knots(j).uFind(boxes(j,2*i)).uIndex();
                k2[j] = this->m_bases[level+1]->knots(j).uFind(boxes(j,2*i+1)).uIndex()+1;
            }

            insert_box(k1,k2,level+1);

            // Build the characteristic matrices (note: call is non-vritual)
            update_structure();

        }
    }

    gsHTensorBasis( gsTensorBSplineBasis<d,T> const& tbasis,
                    gsMatrix<T> const & boxes,
                    const std::vector<index_t> & levels)
    :
    m_manualLevels(false)
    {
        GISMO_ASSERT(boxes.rows() == d, "Points in boxes need to be of dimension d.");
        GISMO_ASSERT(boxes.cols()%2 == 0, "Each box needs two corners but you don't provied gsHTensorBasis constructor with them.");
        GISMO_ASSERT(unsigned (boxes.cols()/2) <= levels.size(), "We don't have enough levels for the boxes.");

        initialize_class(tbasis);

        gsVector<index_t,d> k1;
        gsVector<index_t,d> k2;

        const size_t mLevel = *std::max_element(levels.begin(), levels.end() );
        needLevel( mLevel );

        for(index_t i = 0; i < boxes.cols()/2; i++)
        {
            for(short_t j = 0; j < d; j++)
            {
                k1[j] = m_bases[levels[i]]->knots(j).uFind(boxes(j,2*i)).uIndex();
                k2[j] = m_bases[levels[i]]->knots(j).uFind(boxes(j,2*i+1)).uIndex()+1;
            }

            /* m_boxHistory.push_back( box(k1,k2,levels[i]) );  */
            this->m_tree.insertBox(k1,k2, levels[i]);

            // Build the characteristic matrices (note: call is non-vritual)
            update_structure();
        }
    }

    gsHTensorBasis( const gsHTensorBasis & o) : gsBasis<T>(o)
    {
        this->operator=(o);
    }

    gsHTensorBasis& operator=(const gsHTensorBasis & o)
    {
        if ( this != &o )
        {
            m_xmatrix_offset = o.m_xmatrix_offset;
            m_deg            = o.m_deg;
            m_tree           = o.m_tree;
            m_xmatrix        = o.m_xmatrix;
            m_manualLevels   = o.m_manualLevels; 
            m_uIndices       = o.m_uIndices;

            freeAll( m_bases );
            m_bases.resize( o.m_bases.size() );
            cloneAll(o.m_bases.begin(), o.m_bases.end(), m_bases.begin());
        }
        return *this;
    }

#if EIGEN_HAS_RVALUE_REFERENCES
    gsHTensorBasis(gsHTensorBasis&& other)
    {
        this->operator=(other);
    }

    gsHTensorBasis & operator=(gsHTensorBasis&& other)
    {
        m_deg     = std::move(other.m_deg);
        freeAll( m_bases );
        m_bases   = std::move(other.m_bases);
        m_xmatrix = std::move(other.m_xmatrix);
        m_tree    = std::move(other.m_tree);
        m_xmatrix_offset = std::move(other.m_xmatrix_offset);
        m_manualLevels   = std::move(other.m_manualLevels);
        m_uIndices   = std::move(other.m_uIndices);
        return *this;
    }
#endif

    virtual ~gsHTensorBasis()
    {
        freeAll( m_bases );
    }

    bool manualLevels() const { return m_manualLevels; }

    index_t numLevels() const { return m_bases.size(); }

    void addLevel( const gsTensorBSplineBasis<d, T>& next_basis);
    
    void only_insert_box(point const & k1, point const & k2, int lvl);

protected:

    // TO DO: remove these members after they are not used anymore
    std::vector<int> m_deg;

protected:

    mutable std::vector<tensorBasis*> m_bases;

    std::vector< CMatrix > m_xmatrix;

    hdomain_type m_tree;

    // // Stores the coordinates of all inserted boxes
    // (for debugging purposes)
    // boxHistory m_boxHistory;

    bool m_manualLevels;

    std::vector<std::vector<std::vector<index_t>>> m_uIndices;



public:
    // Needed since m_tree is 16B aligned
#   define Eigen gsEigen
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
#   undef Eigen
    protected:

    std::vector<index_t> m_xmatrix_offset;

    //------------------------------------

public:

    const std::vector< CMatrix >& getXmatrix() const
    {
        return m_xmatrix;
    }

    const std::vector<tensorBasis*>& getBases() const
    {
        return m_bases;
    }

    virtual short_t dim() const
    { return d; }

    int numBreaks(int lvl, int k) const
    {
        return m_tree.numBreaks(lvl,k);
    }

    int numKnots(int lvl, int k) const
    {
        needLevel(lvl);

        return m_bases[lvl]->knots(k).size();
    }

    T knot(int lvl, int k, int i) const
    {
        needLevel(lvl);

        return m_bases[lvl]->component(k).knot(i);
        //return m_bases[lvl]->knot(k,i);
    }

    virtual void anchors_into(gsMatrix<T> & result) const
    {
        result.resize( d, this->size()) ;
        unsigned k(0);

        gsVector<index_t, d> ind;
        for(size_t i = 0; i < m_xmatrix.size(); i++)
        {
            for( CMatrix::const_iterator it =
                     m_xmatrix[i].begin(); it != m_xmatrix[i].end(); it++)
            {
                ind = m_bases[i]->tensorIndex(*it);
                for ( short_t r = 0; r!=d; ++r )
                    result(r,k) = m_bases[i]->knots(r).greville( ind[r] );
                k++;
            }
        }
    }

    virtual void anchor_into(index_t i, gsMatrix<T> & result) const
    {
        index_t lvl = levelOf(i);
        index_t ind = flatTensorIndexOf(i);

        m_bases[lvl]->anchor_into(ind,result);
    }

    virtual void connectivity(const gsMatrix<T> & nodes, gsMesh<T> & mesh) const;
    void connectivity(const gsMatrix<T> & nodes, int level, gsMesh<T> & mesh) const;

    // Prints the characteristic matrices (ie. the indices of all basis
    // functions in the basis)
    void printCharMatrix(std::ostream &os = gsInfo) const
    {
        os<<"Characteristic matrix:\n";
        for(size_t i = 0; i!= m_xmatrix.size(); i++)
        {
            if ( m_xmatrix[i].size() )
            {
                os<<"- level="<<i<<
                    ", size="<<m_xmatrix[i].size() << ":\n";
                os << "("<< m_bases[i]->tensorIndex(*m_xmatrix[i].begin()).transpose() <<")";
                for( CMatrix::const_iterator  it =
                         m_xmatrix[i].begin()+1; it != m_xmatrix[i].end(); it++)
                {
                    os << ", ("<< m_bases[i]->tensorIndex(*it).transpose() <<")";
                }
                os <<"\n";
            }
            else
            {
                os<<"- level="<<i<<" is empty.\n";
            }
        }
    }

    void printSpaces(std::ostream &os = gsInfo) const
    {
        os<<"Spline-space hierarchy:\n";
        for(size_t i = 0; i!= m_xmatrix.size(); i++)
        {
            if ( m_xmatrix[i].size() )
            {
                os<<"- level="<<i<<
                    ", size="<<m_xmatrix[i].size() << ":\n";
                os << "Space: "<< * m_bases[i] <<")\n";
                if (m_manualLevels)
                {
                    os << "Indices:\n";
                    for (size_t dim=0; dim!=d; dim++)
                        os << "Dir "<<dim<<": "<<gsAsConstVector<index_t>(m_uIndices[i][dim]).transpose()<<"\n";
                }
            }
            else
            {
                os<<"- level="<<i<<" is empty.\n";
            }
        }
    }

    void printBases(std::ostream &os = gsInfo) const
    {
        os<<"Spline-space hierarchy:\n";
        for(unsigned i = 0; i< m_bases.size(); i++)
        {
            os<<"- level="<<i<<
                ", size="<<m_bases[i]->size() << ":\n";
            os << "Space: "<< * m_bases[i] <<")\n";
            if (m_manualLevels)
            {
                os << "Indices:\n";
                for (size_t dim=0; dim!=d; dim++)
                    os << "Dir "<<dim<<": "<<gsAsConstVector<index_t>(m_uIndices[i][dim]).transpose()<<"\n";
            }
        }
    }

    void printBasic(std::ostream &os = gsInfo) const
    {
        os << "basis of dimension " <<this->dim()<<
            ",\nlevels="<< this->m_tree.getMaxInsLevel()+1 <<", size="
           << this->size()<<", tree_nodes="<< this->m_tree.size()
            //   << ", leaf_nodes="<< this->m_tree.leafSize();
            //const std::pair<int,int> paths  = this->m_tree.minMaxPath();
            //os << ", path lengths=("<<paths.first<<", "<<paths.second
           << ").\n";
        const gsMatrix<T> supp  = this->support();
        os << "Domain: ["<< supp.col(0).transpose()<< "]..["<<
            supp.col(1).transpose()<< "].\n";
        os <<"Size per level: ";
        for(size_t i = 0; i!= m_xmatrix.size(); i++)
            os << this->m_xmatrix[i].size()<< " ";
        os<<"\n";
    }

    // Look at gsBasis.h for the documentation of this function
    virtual void active_into(const gsMatrix<T> & u, gsMatrix<index_t>& result) const;


    // Look at gsBasis.h for the documentation of this function
    gsMatrix<index_t> allBoundary( ) const;

    // Look at gsBasis.h for the documentation of this function
    virtual gsMatrix<index_t> boundaryOffset(boxSide const & s, index_t offset ) const;

    virtual gsMatrix<index_t> boundaryOffsetLevel(boxSide const & s, index_t offset , index_t level) const;

    virtual index_t levelAtCorner(boxCorner const & c) const;

    virtual index_t functionAtCorner(boxCorner const & c) const;
    virtual index_t functionAtCorner(boxCorner const & c, index_t level) const;//..


    // Look at gsBasis.h for the documentation of this function
    //void evalAllDers_into(const gsMatrix<T> & u, int n,
    //                      std::vector<gsMatrix<T> >& result, bool sameElement) const;

    const gsHDomain<d> & tree() const { return m_tree; }

    gsHDomain<d> &       tree()       { return m_tree; }

    void makeCompressed();


    // GISMO_UPTR_FUNCTION_DEC(gsHTensorBasis<d,T>, boundaryBasis, boxSide const &)

    gsMatrix<T> support() const;

    gsMatrix<T> support(const index_t & i) const;

    void elementSupport_into(const index_t i, gsMatrix<index_t, d, 2>& result) const
    {
        index_t lvl = levelOf(i);
        m_bases[lvl]->elementSupport_into(m_xmatrix[lvl][ i - m_xmatrix_offset[lvl] ],
                                          result);
    }

    gsMatrix<T> elementInSupportOf(index_t j) const
    {
        index_t lvl = levelOf(j);
        gsMatrix<index_t,d,2> sup;
        m_bases[lvl]->elementSupport_into(m_xmatrix[lvl][j-m_xmatrix_offset[lvl]], sup);
        std::pair<point,point> box =  m_tree.queryLevelCell(sup.col(0),sup.col(1),lvl);
        for ( short_t i = 0; i!=d; ++i) //get intersection
        {
            box.first[i]  = ( sup(i,0) >= box.first[i]  ? sup(i,0) : box.first[i] );
            box.second[i] = ( sup(i,1) <= box.second[i] ? sup(i,1) : box.second[i]);
        }
        sup.col(0) = (box.first+box.second)/2;
        sup.col(1) = sup.col(0).array() + 1;
        return m_bases[lvl]->elementDom(sup);
    }

    GISMO_UPTR_FUNCTION_PURE(gsHTensorBasis, clone)

    
    index_t size() const;

    int treeSize() const
    {
        return m_tree.size();
    }

    void numActive_into(const gsMatrix<T> & u, gsVector<index_t>& result) const;

    virtual gsBSplineBasis<T> & component(short_t i)
    {
        return m_bases[ this->maxLevel() ]->component(i);
    }

    virtual const gsBSplineBasis<T> & component(short_t i) const
    {
        return m_bases[ this->maxLevel() ]->component(i);
    }

    tensorBasis & tensorLevel(index_t i) const
    {
        needLevel( i );
        return *this->m_bases[i];
    }

    // Refine the basis uniformly by inserting \a numKnots new knots on each knot span
    virtual void uniformRefine(int numKnots = 1, int mul=1, int dir=-1);

    // Refine the basis uniformly and adjust the given matrix of coefficients accordingly
    //virtual void uniformRefine_withCoefs(gsMatrix<T>& coefs, int numKnots = 1, int mul = 1)

    // Refine the basis uniformly and produce a sparse matrix which
    // maps coarse coefficient vectors to refined ones
    //virtual void uniformRefine_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, int numKnots = 1, int mul = 1)

    // Refine the basis uniformly and adjust the given matrix of coefficients accordingly
    void uniformRefine_withCoefs(gsMatrix<T>& coefs, int numKnots = 1, int mul = 1, int dir = -1);

    // Refine the basis and adjust the given matrix of coefficients accordingly
    void refine_withCoefs(gsMatrix<T> & coefs, gsMatrix<T> const & boxes);
    // void unrefine_withCoefs(gsMatrix<T> & coefs, gsMatrix<T> const & boxes);

    void refineElements_withCoefs   (gsMatrix<T> & coefs,std::vector<index_t> const & boxes);
    void refineElements_withTransfer(std::vector<index_t> const & boxes, gsSparseMatrix<T> &transfer);
    void refineElements_withTransfer2(std::vector<index_t> const & boxes, gsSparseMatrix<T> &transfer);

    void refineElements_withCoefs2(gsMatrix<T> & coefs,std::vector<index_t> const & boxes);

    void unrefineElements_withCoefs   (gsMatrix<T> & coefs,std::vector<index_t> const & boxes);
    void unrefineElements_withTransfer(std::vector<index_t> const & boxes, gsSparseMatrix<T> &transfer);

    // Coarsens the basis uniformly by removing \a numKnots knots on each knot span
    virtual void uniformCoarsen(int numKnots = 1);

    // Coarsen the basis uniformly and adjust the given matrix of coefficients accordingly
    void uniformCoarsen_withCoefs(gsMatrix<T>& coefs, int numKnots = 1);

    // see gsBasis for documentation
    void matchWith(const boundaryInterface & bi, const gsBasis<T> & other,
                   gsMatrix<index_t> & bndThis, gsMatrix<index_t> & bndOther) const;

    // see gsBasis for documentation
    void matchWith(const boundaryInterface & bi, const gsBasis<T> & other,
                   gsMatrix<index_t> & bndThis, gsMatrix<index_t> & bndOther, index_t offset) const;

    short_t maxDegree() const
    {
        short_t td = m_bases[0]->degree(0);
        // take maximum of coordinate bases degrees
        for (short_t k=1; k!=d; ++k)
            td = math::max(td, m_bases[0]->degree(k));
        return td;
    }

    short_t minDegree() const
    {
        short_t td = m_bases[0]->degree(0);
        // take maximum of coordinate bases degrees
        for (short_t k=1; k!=d; ++k)
            td = math::min(td, m_bases[0]->degree(k));
        return td;
    }

    void reduceContinuity(int const & i = 1)
    {
        for (unsigned int lvl = 0; lvl <= maxLevel(); lvl++)
        {
            for (unsigned int dir = 0; dir < d; dir++)
            {
                // TODO check: max interior mult + i <= m_p+1

                // We iterate through unique knots, skipping the first and last knot
                // At level 0 we iterate through all unique knots,
                // At level >0 we iterate through all knots that are new, i.e. every other knot starting from 1
                for (typename gsKnotVector<T>::uiterator it =
                         m_bases[lvl]->knots(dir).ubegin() + 1;
                     it < m_bases[lvl]->knots(dir).uend() - 1;
                     it += (lvl == 0? 1 : 2))
                {
                    for(unsigned int j =lvl;j < m_bases.size();j++)
                        m_bases[j]->component(dir).insertKnot(*it,i);

                }
            }
        }
        update_structure();
    }

    virtual inline short_t degree(short_t i) const
    { return m_bases[0]->degree(i);}



    // S.K.
    index_t getLevelAtPoint(const  gsMatrix<T> & Pt ) const;

    // S.K.
    index_t getLevelAtIndex(const point & Pt ) const;

    // S.K.
    void getLevelUniqueSpanAtPoints(const  gsMatrix<T> & Pt,
                                    gsVector<index_t> & lvl,
                                    gsMatrix<index_t> & loIdx ) const;

    unsigned maxLevel() const
    {
        return m_tree.getMaxInsLevel();
    }

    inline index_t levelOf(index_t i) const
    {
        return std::upper_bound(m_xmatrix_offset.begin(),
                                m_xmatrix_offset.end(), i)
            - m_xmatrix_offset.begin() - 1;
    }

/*
  const boxHistory & get_inserted_boxes() const
  {
  return m_boxHistory;
  }
*/

    void degreeElevate(int const & i = 1, int const dir = -1);
    void degreeReduce(int const & i = 1, int const dir = -1);

    void degreeIncrease(int const & i= 1, int const dir = -1);
    void degreeDecrease(int const & i = 1, int const dir = -1);

    virtual void refine(gsMatrix<T> const & boxes, int refExt);
    virtual void unrefine(gsMatrix<T> const & boxes, int refExt);

    std::vector<index_t> asElements(gsMatrix<T> const & boxes, int refExt = 0) const;
    std::vector<index_t> asElementsUnrefine(gsMatrix<T> const & boxes, int refExt = 0) const;

    // std::vector<index_t> asElements(gsMatrix<T> const & boxes, int refExt = 0) const;

    virtual void refine(gsMatrix<T> const & boxes);
    // virtual void unrefine(gsMatrix<T> const & boxes);

    virtual void refineElements(std::vector<index_t> const & boxes);

    virtual void unrefineElements(std::vector<index_t> const & boxes);

    void refineSide(const boxSide side, index_t lvl);

    void refineBasisFunction(const index_t i);

    // Look at gsBasis.h for the documentation of this function
    //virtual void uniformRefine(int numKnots = 1);

    typename gsBasis<T>::domainIter makeDomainIterator() const
    {
        return typename gsBasis<T>::domainIter(new gsHDomainIterator<T, d>(*this));
    }

    typename gsBasis<T>::domainIter makeDomainIterator(const boxSide & s) const
    {
        return ( s == boundary::none ?
                 typename gsBasis<T>::domainIter(new gsHDomainIterator<T, d>(*this)) :
                 typename gsBasis<T>::domainIter(new gsHDomainBoundaryIterator<T, d>(*this,s) )
            );
    }

    inline
    index_t flatTensorIndexOf(const index_t i) const
    {
        return flatTensorIndexOf(i, this->levelOf(i) );
    }

    inline
    index_t flatTensorIndexOf(const index_t i, const index_t level) const
    {
        const index_t offset = this->m_xmatrix_offset[level];
        const index_t ind_in_level = this->m_xmatrix[level][i - offset];
        return ind_in_level;
    }

    std::vector< std::vector< std::vector<index_t > > > domainBoundariesParams( std::vector< std::vector< std::vector< std::vector< T > > > >& result) const;

    std::vector< std::vector< std::vector<index_t > > > domainBoundariesIndices( std::vector< std::vector< std::vector< std::vector<index_t > > > >& result) const;
    // TO DO: use gsHDomainLeafIterator for a better implementation
    size_t numElements(boxSide const & s = 0) const
    {
        typename gsBasis<T>::domainIter domIter =
            s == boundary::none ?
            typename gsBasis<T>::domainIter(new gsHDomainIterator<T, d>(*this)) :
            typename gsBasis<T>::domainIter(new gsHDomainBoundaryIterator<T, d>(*this,s) );

        size_t numEl = 0;
        for (; domIter->good(); domIter->next())
        {
            numEl++;
        }

        return numEl;
    }

    void flatTensorIndexesToHierachicalIndexes(gsSortedVector< int > & indexes,const int level) const;

    int flatTensorIndexToHierachicalIndex(index_t index,const int level) const;

    void activeBoundaryFunctionsOfLevel(const unsigned level,const boxSide & s,std::vector<bool>& actives) const;


    virtual void increaseMultiplicity(index_t lvl, int dir, T knotValue, int mult = 1);

    virtual void increaseMultiplicity(index_t lvl, int dir, const std::vector<T> & knotValue, int mult = 1);

protected:

    virtual void update_structure(); // to do: rename as updateCharMatrices

    void needLevel(int maxLevel) const;

    void createMoreLevels(int numLevels) const;

    void getBoxesAlongSlice( int dir, T par,std::vector<index_t>& boxes ) const;

protected:

    void _knotIndexToDiadicIndex(const index_t level, const index_t dir, index_t & knotIndex) const;
    void _knotIndexToDiadicIndex(const index_t level, gsVector<index_t,d> & diadicIndex) const;

    void _diadicIndexToKnotIndex(const index_t level, gsVector<index_t,d> & diadicIndex) const;
    void _diadicIndexToKnotIndex(const index_t level, const index_t dir, index_t & diadicIndex) const;

private:

    void insert_box(point const & k1, point const & k2, int lvl);

    void initialize_class(gsBasis<T> const&  tbasis);

    void functionOverlap(const point & boxLow, const point & boxUpp,
                         const int level, point & actLow, point & actUpp);

    // \brief Sets all functions of \a level to active or passive- one by one
    void set_activ1(int level);

    // \brief Computes the set of active basis functions in the basis
    void setActive();

    // \brief Computes the connectivity for a level on a mesh that has all vertices
    void addConnectivity(int level, gsMesh<T> & mesh) const;

    virtual gsSparseMatrix<T> coarsening(const std::vector<CMatrix>& old,
                                         const std::vector<CMatrix>& n,
                                         const gsSparseMatrix<T,RowMajor> & transfer) const = 0;

    virtual gsSparseMatrix<T> coarsening_direct(const std::vector<gsSortedVector<index_t> >& old,
                                                const std::vector<gsSortedVector<index_t> >& n,
                                                const std::vector<gsSparseMatrix<T,RowMajor> >& transfer) const = 0;

    virtual gsSparseMatrix<T> coarsening_direct2(const std::vector<gsSortedVector<index_t> >& old,
                                                 const std::vector<gsSortedVector<index_t> >& n,
                                                 const std::vector<gsSparseMatrix<T,RowMajor> >& transfer) const = 0;

    std::vector< std::vector< std::vector<index_t > > > domainBoundariesGeneric(std::vector< std::vector< std::vector< std::vector<index_t > > > >& indices,
                                                                                      std::vector< std::vector< std::vector< std::vector< T > > > >& params,
                                                                                      bool indicesFlag ) const;

public:
    void transfer (const std::vector<gsSortedVector<index_t> > &old, gsSparseMatrix<T>& result);

    void transfer2 (const std::vector<gsSortedVector<index_t> > &old, gsSparseMatrix<T>& result);

    void setActiveToLvl(int level, std::vector<CMatrix>& x_matrix_lvl) const;

public:
    bool testPartitionOfUnity(const index_t npts = 100, const T tol = 1e-12) const;

//    void local2globalIndex( gsVector<index_t,d> const & index,
//                            index_t lvl,
//                            gsVector<index_t,d> & result
//        ) const;

//    void global2localIndex( gsVector<index_t,d> const & index,
//                            index_t lvl,
//                            gsVector<index_t,d> & result
//        ) const;

}; // class gsHTensorBasis

// Next line disallows instantization of gsTensorBasis<0,T>
template<typename T> class gsHTensorBasis<0,T>
{using T::GISMO_ERROR_gsHTensorBasis_cannot_have_dimension_zero;};


#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsHTensorBasis2(pybind11::module &m);
  void pybind11_init_gsHTensorBasis3(pybind11::module &m);
  void pybind11_init_gsHTensorBasis4(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11


} // namespace gismo

// ************************************************
// ************************************************

#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsHTensorBasis.hpp)
#endif