gsHDomain.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>
#include <gsHSplines/gsHDomainLeafIter.h>
#include <gsCore/gsBoundary.h>

namespace gismo 
{

template<typename T, unsigned d>
class gsHDomainIterator;

template<class T>
class gsSegment;

// template <short_t d, class Z>
// class gsKdNode;

template <class T>
class gsVSegment;

template<short_t d, class Z = index_t>
class gsHDomain
{
public:
    typedef gsKdNode<d, Z> node;

    typedef typename node::point point; // it's a gsVector<index_t,d>

    typedef typename node::kdBox box; // it's a gsAABB<d,unsigned>

    typedef gsHDomainLeafIter<node,false> literator;

    typedef gsHDomainLeafIter<node,true> const_literator;

    template <class U, unsigned dd> friend class gsHDomainIterator;

private:

    struct numLeaves_visitor;
    struct numNodes_visitor;
    struct levelUp_visitor;

private:

    node * m_root;

    point m_upperIndex;
#   define Eigen gsEigen
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
#   undef Eigen

    unsigned m_indexLevel;

    unsigned m_maxInsLevel;

    unsigned m_maxPath;

public:

    gsHDomain() : m_indexLevel(0)
    {
        m_root = nullptr;
        m_maxInsLevel = 0;
        m_maxPath = 0;
    }
    
    gsHDomain(point const & upp)
    {
        m_root = nullptr;
        m_maxInsLevel = 0;
        m_maxPath = 0;
        init(upp);
    }

    gsHDomain( const gsHDomain & o) :
        m_upperIndex(o.m_upperIndex),
        m_indexLevel(o.m_indexLevel),
        m_maxInsLevel(o.m_maxInsLevel),
        m_maxPath(o.m_maxPath)
    {
        m_root = new node(*o.m_root);
    }

    gsHDomain& operator=( const gsHDomain & o)
    {
        if ( this == &o )
            return *this;
        
        m_root = new node(*o.m_root);

        m_upperIndex  = o.m_upperIndex;
        m_indexLevel  = o.m_indexLevel;
        m_maxInsLevel = o.m_maxInsLevel;
        m_maxPath    = o.m_maxPath;

        return *this;
    }

#if EIGEN_HAS_RVALUE_REFERENCES
    gsHDomain(gsHDomain&& o) :
    m_root(o.m_root),
    m_upperIndex(std::move(o.m_upperIndex)),
    m_indexLevel(o.m_indexLevel),
    m_maxInsLevel(o.m_maxInsLevel),
    m_maxPath(o.m_maxPath)
    {
        o.m_root = nullptr;
    }

    gsHDomain & operator=(gsHDomain&& o)
    {
        delete m_root; m_root = o.m_root; o.m_root = nullptr;
        m_upperIndex  = std::move(o.m_upperIndex);
        m_indexLevel  = o.m_indexLevel;
        m_maxInsLevel = o.m_maxInsLevel;
        m_maxPath     = o.m_maxPath;
        return *this;
    }
#endif

    void init(point const & upp, unsigned index_level)
    {
        m_indexLevel = index_level;
        m_maxInsLevel = 0;

        if ( m_root )
            delete m_root;

        for (short_t i=0; i<d; ++i)
            m_upperIndex[i] = (upp[i]<< m_indexLevel);

        m_root = new node(m_upperIndex);
        m_maxPath = 1;
    }

    void init(point const & upp)
    {
        // Idea: find index_level so that for each i, upp[i] * (2 ^ index_level) does not overflow in Z.
        // See issue #672 for more details.

        // backwards compatibility
        Z oldMax = 13;

        Z numMax = std::numeric_limits<Z>::max();

        std::vector<Z> logUpps(d);
        for(short_t i=0; i<d; i++)
        {
            // prevent division by zero
            if(upp[i] == 1)
                logUpps[i] = oldMax;
            else
            {
                // floor of log_2 (numMax / upp[i]):
                logUpps[i] = math::floor( (math::log(numMax) - math::log(upp[i])) / math::log(2) );
            }
        }

        // If the computed number would be too big we take 13 as we used to.
        init(upp, std::min( *std::min_element(logUpps.begin(), logUpps.end()), oldMax) );
    }

    ~gsHDomain() { delete m_root; }

    gsHDomain * clone() const;

public:

    void computeFinestIndex(gsVector<Z, d> const & index,
                            unsigned lvl,
                            gsVector<Z, d> & result) const;

    void computeLevelIndex(gsVector<Z, d> const & index,
                           unsigned lvl,
                           gsVector<Z, d> & result) const;

    void local2globalIndex(gsVector<Z, d> const & index,
                           unsigned lvl,
                           gsVector<Z, d> & result) const;

    void global2localIndex(gsVector<Z, d> const & index,
                           unsigned lvl,
                           gsVector<Z, d> & result) const;

    const point & upperCorner() const
    {
        return m_upperIndex;
    }

    const point upperCornerIndex() const
    {
        point ind = m_upperIndex;
        for (short_t i=0; i<d; ++i)
            ind[i] = (ind[i] >> m_indexLevel);
        return ind;
    }

    void insertBox (point const & lower, point const & upper,
                    node * _node, int lvl);

    void insertBox (point const & lower, point const & upper, int lvl)
    { insertBox(lower, upper, m_root, lvl); }

    void clearBox (point const & lower, point const & upper,
                    int lvl);

    void sinkBox (point const & lower, point const & upper, int lvl);

    void internalIndex (point const & point_idx, int lvl, point & internal_idx)
    {
        for ( short_t i = 0; i!=d; ++i )
            internal_idx[i] = point_idx[i] << (m_indexLevel-lvl);
    }


    bool query1 (point const & lower, point const & upper,
                 int level, node * _node ) const;

    bool query1 (point const & lower, point const & upper,
                 int level) const;

    bool query2 (point const & lower, point const & upper,
                 int level, node *_node ) const;

    bool query2 (point const & lower, point const & upper,
                 int level) const;

    int query3(point const & k1, point const & k2, 
               int level, node *_node ) const;

    int query3(point const & lower, point const & upper,
               int level) const;

    int query4(point const & lower, point const & upper,
               int level, node  *_node) const;

    int query4(point const & lower, point const & upper,
               int level) const;

    std::pair<point,point> queryLevelCell(point const & lower,
                                          point const & upper,
                                          int level) const;

    int levelOf(point const & p, int level) const
    { return pointSearch(p,level,m_root)->level;}

    void incrementLevel();
    
    void multiplyByTwo();

    void divideByTwo();

    void decrementLevel();

    literator beginLeafIterator()
    {
        return literator(m_root, m_indexLevel);
    }

    const_literator beginLeafIterator() const
    {
        return const_literator(m_root, m_indexLevel);
    }

    void makeCompressed();
    
    int size() const;

    int numBreaks(int lvl, int k) const
    {
        return m_upperIndex[k] >> (m_indexLevel - lvl);
    }

    int leafSize() const;

    std::pair<int,int> minMaxPath() const;

    void printLeaves() const;

    // getBoxes-functions might get removed at some point of time.
    // Use iterators instead whenever possible.
    void getBoxes(gsMatrix<Z>& b1,
                  gsMatrix<Z>& b2,
                  gsVector<Z>& level) const;

    // getBoxes-functions might get removed at some point of time.
    // Use iterators instead whenever possible.
    void getBoxesOnSide(boundary::side s,
                        gsMatrix<Z>& b1,
                        gsMatrix<Z>& b2,
                        gsVector<Z>& level) const;



    // getBoxes-functions might get removed at some point of time.
    // Use iterators instead whenever possible.
    void getBoxesInLevelIndex(gsMatrix<Z>& b1,
                              gsMatrix<Z>& b2,
                              gsVector<index_t>& level) const;

    std::vector<std::vector<std::vector<std::vector<Z>>>> getPolylines() const;

    std::vector<std::vector<std::vector<Z>>> getPolylinesSingleLevel(std::vector<gsVSegment<Z>>& seg) const;


    inline unsigned getIndexLevel() const
    {
        return m_indexLevel;
    }

    inline unsigned getMaxInsLevel() const
    {
        return m_maxInsLevel;
    }

    void computeMaxInsLevel();

private:
    
    static bool haveOverlap(box const & box1, box const & box2);

    static bool isContained(box const & box1, box const & box2);

    static std::pair<point,point> select_part(point const & k1, point const & k2,
                                              point const & k3, point const & k4);
    
    static void bisectBox(box const & original, 
                          int k, Z coord,
                          box & leftBox,
                          box & rightBox );

    static void setLevel(node *_node, int lvl);

    static bool isDegenerate(box const & someBox);

    void setIndexLevel(int) const
    {
        GISMO_NO_IMPLEMENTATION
    }

    // getBoxes-functions might get removed at some point of time.
    // Use iterators instead whenever possible.
    void getBoxes_vec(std::vector<std::vector<Z>>& boxes) const;

    void connect_Boxes(std::vector<std::vector<Z> > &boxes) const;
    void connect_Boxes2d(std::vector<std::vector<Z> > &boxes) const;

    void connect_Boxes_2(std::vector<std::vector<Z> > &boxes) const;

    void getRidOfOverlaps( std::list< std::list< gsVSegment<Z> > >& vert_seg_lists ) const;

    void sweeplineConnectAndMerge( std::vector< std::vector< std::vector<Z> > >& result,
                                   std::list< std::list< gsVSegment<Z> > >& vert_seg_lists ) const;
    

private:

    template<typename visitor>
    typename visitor::return_type
    boxSearch(point const & k1, point const & k2, 
              int level, node  *_node) const;

    template<typename visitor>
    typename visitor::return_type
    leafSearch() const;

    template<typename visitor>
    typename visitor::return_type
    nodeSearch() const;

    node * pointSearch(const point & p, int level, node  *_node) const;

    struct maxLevel_visitor
    {
        typedef int return_type;
        static return_type init() {return 0;}
        
        static void visitLeaf(gsKdNode<d, Z> * leafNode, return_type &i)
        {
            if (leafNode->level>i) i=leafNode->level;
        }
    };
    
    struct levelUp_visitor
    {
        typedef int return_type;
        static return_type init() {return 0;}
        
        static void visitLeaf(gsKdNode<d, Z> * leafNode, return_type &)
        {
            leafNode->level++;
        }
    };

    struct levelDown_visitor
    {
        typedef int return_type;
        static return_type init() {return 0;}
        
        static void visitLeaf(gsKdNode<d, Z> * leafNode, return_type &)
        {
            leafNode->level--;
        }
    };

    struct numLeaves_visitor
    {
        typedef int return_type;
        static return_type init() {return 0;}
        
        static void visitLeaf(gsKdNode<d, Z> * , return_type & i)
        {
            i++;
        }
    };

    struct numNodes_visitor
    {
        typedef int return_type;
        static return_type init() {return 0;}
        
        static void visitNode(gsKdNode<d, Z> * , return_type & i)
        {
            i++;
        }
    };

    struct liftCoordsOneLevel_visitor
    {
        typedef int return_type;
        static return_type init() {return 0;}
        
        static void visitNode(gsKdNode<d, Z> * leafNode, return_type &)
        {
            leafNode->multiplyByTwo();
        }
    };

    struct reduceCoordsOneLevel_visitor
    {
        typedef int return_type;
        static return_type init() {return 0;}

        static void visitNode(gsKdNode<d, Z> * leafNode, return_type &)
        {
            leafNode->divideByTwo();
        }
    };

    struct printLeaves_visitor
    {
        typedef int return_type;
        static return_type init() {return 0;}
        
        static void visitLeaf(gsKdNode<d, Z> * leafNode, return_type &)
        {
            gsInfo << *leafNode;
        }
    };

    // TODO: stop traverse as soon as it is found for the first time..
    struct get_cell_visitor
    {
        typedef std::pair<point,point> return_type;

        // initialize result
        static return_type init()
        {
            return std::make_pair(point::Zero(),point::Zero());
            //return return_type();//!does not properly initialize the points
        }

        static void visitLeaf(gismo::gsKdNode<d, Z> * leafNode , int level, return_type & res)
        {
            if ( leafNode->level == level )
            {
                res.first  = leafNode->lowCorner();
                res.second = leafNode->uppCorner();
            }
        }
    };

};


}// end namespace gismo


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