gsKnotVector.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsForwardDeclarations.h>
#include <gsCore/gsDomain.h>

#include <gsNurbs/gsKnotIterator.h>

namespace gismo
{

template<typename T>
class gsKnotVector : public gsDomain<T>
{
public: // typedefs

    // IMPORTANT: mult_t must be a signed type, otherwise the
    // knotIterator business would not work.
    typedef index_t mult_t;

    // container types
    typedef std::vector<T>      knotContainer;
    typedef std::vector<mult_t> multContainer;

    // TODO: remove the non-const versions and think about the names.
    // E.g., reverse_uiterator clashes with urbegin().

    //typedef const T * iterator;
    typedef typename std::vector<T>::const_iterator iterator;
    typedef std::reverse_iterator<iterator> reverse_iterator;

    typedef     internal::gsUKnotIterator<T> uiterator;
    typedef std::reverse_iterator<uiterator> reverse_uiterator;

    typedef           internal::gsKnotIterator<T>  smart_iterator;
    typedef std::reverse_iterator<smart_iterator>  reverse_smart_iterator;

    //compatibility iterator typedefs
    typedef uiterator        const_uiterator;
    typedef iterator         const_iterator;
    typedef reverse_iterator const_reverse_iterator;
    
public: // iterator ends

    iterator               begin()   const;
    iterator               end()     const;    
    reverse_iterator       rbegin()  const;    
    reverse_iterator       rend()    const;

    iterator               beginAt(mult_t upos)   const;

    iterator               endAt(mult_t upos)     const;    
    
    uiterator              ubegin()  const;
    uiterator              uend()    const;
    reverse_uiterator      urbegin() const;
    reverse_uiterator      urend()   const;

    smart_iterator         sbegin()  const;
    smart_iterator         send()    const;    
    reverse_smart_iterator rsbegin() const;
    reverse_smart_iterator rsend()   const;

public: // constructors

    gsKnotVector() : m_deg(-1)
    { }

    explicit gsKnotVector(knotContainer knots, short_t degree = -1 );

    void swap( gsKnotVector& other );

    gsAsConstMatrix<T> asMatrix() const {return m_repKnots;}
    
public:
    gsKnotVector* clone() const;

public: // inserting and removing

    void insert(T knot, mult_t mult = 1 );

    template<typename iterType>
    void insert( iterType ibeg, iterType iend )
    {
        //size_t numKnots = size() + (iend - ibeg);
        // GISMO_ENSURE( numKnots < std::numeric_limits<mult_t>::max(),
        //               "Too many knots." );
        knotContainer temp(m_repKnots.size()+(iend-ibeg) );
        std::merge( begin(), end(), ibeg, iend, temp.begin() );

        m_repKnots.swap(temp);

        rebuildMultSum(); // Possibly sub-optimal but the code is clean.

        GISMO_ASSERT( check(), "Unsorted knots or invalid multiplicities.");
    }

    //TODO: insert with uniqIter and multiplicity.

    void remove( uiterator uit, mult_t mult = 1 );

    void remove( const T knot, mult_t mult = 1 );

    //TODO: remove with two iterators.
    
public: // multiplicities

    mult_t multiplicity( T u ) const;

    mult_t multFirst() const { return m_multSum.front(); }

    mult_t multLast() const { return m_multSum.back() - m_multSum.end()[-2]; }

    mult_t maxInteriorMultiplicity() const;

    mult_t minInteriorMultiplicity() const;

    mult_t multiplicityIndex( mult_t i ) const;

public: // modifiers

    void affineTransformTo( T newBeg, T newEnd );

public: // queries

    inline size_t size() const { return m_repKnots.size(); }

    inline size_t uSize() const { return m_multSum.size(); }

    const T& operator[](const mult_t i) const
    {
        GISMO_ASSERT( static_cast<size_t>(i) < m_repKnots.size(),
                      "Index "<<i<<" not in the knot vector.");
        return m_repKnots[i];
    }

    const T& operator()(const mult_t i) const
    {
        GISMO_ASSERT(i+numLeftGhosts() >=0  && static_cast<size_t>(i+numLeftGhosts()) < uSize(),
                              "Unique index "<<i<<" not in the knot vector.");
        return *( this->ubegin()+(numLeftGhosts()+i) );
    }

    inline T  uValue(const size_t & i) const
    { return this->operator()(i); }

    inline size_t numElements() const { return (domainUEnd() - domainUBegin()); }

public: // getters

    knotContainer unique () const
    {
        return knotContainer(this->ubegin(),this->uend());
    }

    multContainer multiplicities() const;

    operator const knotContainer& () const {return m_repKnots;}

    const T * data() const {return m_repKnots.data(); }

    const mult_t * multSumData() const {return m_multSum.data(); }

public: // Findspan and value query

    uiterator uFind( const T u ) const;

    iterator iFind( const T u ) const;

    uiterator uUpperBound( const T u ) const;

    uiterator uLowerBound( const T u ) const;
    
public: // miscellaneous

    std::ostream &print(std::ostream &os = gsInfo ) const;

    bool check() const;

    typedef typename knotContainer::iterator nonConstIterator    ;
    typedef typename multContainer::iterator nonConstMultIterator;

    iterator domainBegin() const
    {
        GISMO_ASSERT( size() > static_cast<size_t>(2*m_deg+1), "Not enough knots.");
        return begin() + m_deg;
    }

    iterator domainEnd() const
    {
        GISMO_ASSERT( size() > static_cast<size_t>(2*m_deg+1), "Not enough knots.");
        return end() - (m_deg + 1);
    }

    smart_iterator domainSBegin() const
    { return sbegin() + m_deg; }

    smart_iterator domainSEnd() const
    { return send() - (m_deg + 1); }

    uiterator domainUBegin() const
    {
        return domainSBegin().uIterator();
        // equivalent:
        //return ubegin() + domainSBegin().uIndex();
    }

    uiterator domainUEnd() const
    {
        return domainSEnd().uIterator();
        // equivalent:
        //return ubegin() + domainSEnd().uIndex();
    }

    inline bool inDomain(const T u) const
    {
        return u >= *domainBegin() && u <= *domainEnd();
        // equivalent:
        // return u >= m_repKnots[m_deg] && u <= m_repKnots.end()[-m_deg-1];        
    }

    void erase(const mult_t first, const mult_t last);

    void trimLeft (const mult_t numKnots);
    
    void trimRight(const mult_t numKnots);

    index_t numLeftGhosts() const
    {
        smart_iterator it(*this,0,0);
        it += math::min( (size_t)m_deg, size() );
        return std::distance( uiterator(*this,0,0), it.uIterator() );
    }

    index_t numRightGhosts() const
    {
        return std::distance(domainUEnd(), uend()) - 1;
    }

public:

    static bool isConsistent(const knotContainer & repKnots,
                             const multContainer & multSums);

    inline bool operator== (const gsKnotVector<T>& other) const
    {
        return (m_repKnots == other.m_repKnots) &&
               (m_multSum  == other.m_multSum ) &&
               (m_deg == other.m_deg);
    }

private: // internals

    // Creates a vector of multiplicity sums from a scratch and m_repKnots.
    void rebuildMultSum();

private: // members

    // Knots including repetitions.
    knotContainer m_repKnots;

    // m_multSum[i] = cardinality of { knots <= unique knot [i] }.
    multContainer m_multSum;




//=======================================================================//
//+++++++ Here follow the functions I would like to get rid of. +++++++++//
//=======================================================================//



public: // Deprecated functions required by gsKnotVector.

    explicit gsKnotVector(short_t degree)
    {
        m_deg = degree;
    }

    gsKnotVector( T first,
                  T last,
                  unsigned interior,
                  mult_t mult_ends=1,
                  mult_t mult_interior=1,
                  short_t degree = -1 );

    template<typename iterType>
    gsKnotVector(short_t deg, const iterType begOfKnots, const iterType endOfKnots)
    {
        insert(begOfKnots,endOfKnots);
        m_deg = (deg == - 1 ? deduceDegree() : deg);
    }

public:
    // TODO If stays, make it private.
    
    void initUniform( T first,
                      T last,
                      unsigned interior,
                      unsigned mult_ends,
                      unsigned mult_interior,
                      short_t degree=-1);

    void initUniform( unsigned numKnots,
                      unsigned mult_ends,
                      unsigned mult_interior = 1,
                      short_t degree = - 1)
    {
        initUniform(0.0, 1.0, numKnots - 2, mult_ends, mult_interior, degree );
    }


    void initGraded(unsigned numKnots, int degree,
                    T grading = 0.5, unsigned mult_interior = 1)
    {
        initGraded( 0.0, 1.0, numKnots - 2, degree, grading, mult_interior);
    }

    void initGraded(T u0, T u1, unsigned interior, int degree,
                    T grading, unsigned mult_interior = 1)
    {
        GISMO_ASSERT(u0<u1,"Knot vector must be an interval.");

        m_deg = degree;
        m_repKnots.reserve( 2*(m_deg+1) + interior*mult_interior );
        m_multSum .reserve(interior+2);

        const T h = (u1-u0) / (T)(interior+1);

        m_repKnots.insert(m_repKnots.begin(), m_deg+1, u0);
        m_multSum .push_back(m_deg+1);

        for ( unsigned i=1; i<=interior; i++ )
        {
            m_repKnots.insert(m_repKnots.end(), mult_interior,
                              math::pow(T(i)*h, 1.0/grading) );
            m_multSum .push_back( mult_interior + m_multSum.back() );
        }
        m_repKnots.insert(m_repKnots.end(), m_deg+1, u1);
        m_multSum .push_back( m_deg+1 + m_multSum.back() );
    }

    gsMatrix<T> greville() const
    {
        gsMatrix<T> gr( 1,this->size() - m_deg - 1 );
        this->greville_into(gr);
        return gr;
    }

    // TODO Write properly.
     
/*
    void makeCompressed(const T tol = 1e-7)
    {
        // TO DO: Check and improve, use tmp storage
        nonConstIterator     k = m_repKnots.begin();
        nonConstMultIterator m = m_multSum  .begin();
        T prev  = *k;
        while( k+1 != end() )
        {
            if ( math::abs( *k - *(k+1) ) <= tol )
            {
                *m = *(m+1) ;
                m_repKnots.erase(k+1);//all app.
                m_multSum.erase(m+1);
            }
            else
            {
                k++;
                m++;
            }
        }
        m_repKnots.resize( k - m_repKnots.begin() );
        m_multSum .resize( m - m_multSum .begin() );
    }
*/



    bool includes(const gsKnotVector<T> & other) const;

    void difference(const gsKnotVector<T> & other,
                   std::vector<T>& result) const;

    void symDifference(const gsKnotVector<T> & other,
                       std::vector<T>& result) const;

    gsKnotVector knotUnion(const gsKnotVector & b) const;

    gsKnotVector knotIntersection(const gsKnotVector & b) const;

    void reverse();

    void insert( const knotContainer &knots, int mult = 1 )
    {
        for( int i = 0; i < mult; ++i)
            insert( knots.begin(), knots.end() ); // inefficient
    }

    void set_degree(short_t p)
    {
        m_deg = p;
    }

    void supportIndex_into(const mult_t &i, gsMatrix<index_t>& result) const;
    // TODO: If the function stays, make a unit test from what is in unifiedKnotVector.cpp.

    void uniformRefine( mult_t numKnots = 1, mult_t mult = 1);

    template<typename iterType>
    void refineSpans( iterType begin, iterType end, mult_t knotsPerSpan )
    {
        // We sort the input to be on the safe side.
        multContainer input;
        input.insert(input.begin(), begin, end); // Pitfall: using std::copy requires reserve or somesuch beforehand.
        std::sort(input.begin(), input.end());

        knotContainer newKnots;
        T segmentsPerSpan = knotsPerSpan + 1;
        T newKnot;
        T spanBegin;
        T spanEnd;
        for( typename multContainer::const_iterator it= input.begin();
             it != input.end();
             ++it )
            for( mult_t k = 1; k <= knotsPerSpan; ++k )
            {
                spanBegin =*(this->ubegin()+*it);
                spanEnd  =*(this->ubegin()+*it+1);
                newKnot = ( (segmentsPerSpan-(T)(k)) * spanBegin + (T)(k) * spanEnd ) / segmentsPerSpan;
                newKnots.push_back( newKnot );
            }

        insert(newKnots.begin(), newKnots.end());
    }

    void addConstant( T amount );

    void addConstant( T start, T amount );

public: // things required by gsKnotVector

    gsKnotVector( const knotContainer& uKnots,
                         int degree,
                         int regularity);

    void initClamped(int degree, unsigned numKnots = 2, unsigned mult_interior = 1);

    void initClamped(T u0, T u1, int degree, unsigned interior = 0,
                     unsigned mult_interior = 1)
    {
        return initUniform(u0, u1, interior, degree + 1, mult_interior, degree );
    }

    int degree () const;


    void greville_into(gsMatrix<T> & result) const;

    void centers_into(gsMatrix<T> & result) const;

    int deduceDegree() const;

    T greville(int i) const;

    T first () const
    {
        GISMO_ASSERT(this->size()>=1, "I need at least one knot.");
        return m_repKnots.front();
    }
    
    T last  () const
    {
        GISMO_ASSERT(this->size()>=1, "I need at least one knot.");
        return m_repKnots.back();
    }

    void transform(T c, T d)
    {
        affineTransformTo(c,d);
    }

    void refineSpans( const std::vector<unsigned> & spanIndices, mult_t knotsPerSpan = 1)
    {
        multContainer transformedIndices;
        transformedIndices.reserve(spanIndices.size());
        for( std::vector<unsigned>::const_iterator it = spanIndices.begin();
             it != spanIndices.end();
             ++it )
            transformedIndices.push_back(static_cast<mult_t>(*it));

        return refineSpans( transformedIndices, knotsPerSpan );
    }

    void refineSpans( const multContainer & spanIndices, mult_t knotsPerSpan = 1);

    void degreeIncrease(int const & i = 1)
    {
        // update knots
        m_repKnots.reserve(size()+2*i);
        m_repKnots.insert(m_repKnots.begin(), i, m_repKnots.front() );
        m_repKnots.insert(m_repKnots.end()  , i, m_repKnots.back()  );

        // update multiplicity sum
        std::transform(m_multSum.begin(), m_multSum.end(), m_multSum.begin(),
                       GS_BIND2ND(std::plus<mult_t>(), i) );
        m_multSum.back() += i;

        m_deg += i;
    }

    void degreeDecrease(int const & i = 1, bool updateInterior = false)
    {
        remove( ubegin()  , i );
        remove( uend() - 1, i );
        m_deg -= i;

        if (updateInterior)
        {        
            if ( m_deg <= 0 )
                initUniform(first(), last(), 0, 1, 0, 0);
            else
                for (uiterator itr = ubegin()+1; itr != uend()-1; ++itr)
                    if ( itr.multiplicity() > m_deg )
                        remove( itr, itr.multiplicity() - m_deg );
        }
    }

    void increaseMultiplicity(const mult_t i = 1, bool boundary = false);

    void reduceMultiplicity(const mult_t i = 1, bool boundary = false);

    template<typename iterType>
    void append ( iterType begin, iterType end )
    {
        insert( begin, end );
    }

    //void _stretchEndKnots()

    const knotContainer& get() const // to be removed since we have implicit cast
    {
        return m_repKnots;
    }

public: // Deprecated functions required by gsCompactKnotVector.

    inline bool has(T knot) const
    {
        return std::binary_search( ubegin(), uend(), knot);
        // equivalent:
        // return 0 != multiplicity(knot);
    }

    unsigned u_multiplicityIndex(size_t const & i) const
    {
        return (ubegin() + i).multiplicity();
    }

    inline unsigned firstKnotIndex(const size_t & i) const
    {
        return (ubegin()+i).firstAppearance();
        // equivalent:
        // return 0 == i ? 0 : m_multSum[i-1];
    }
     
    inline unsigned lastKnotIndex(const size_t & i) const
    {
        return (ubegin()+i).lastAppearance();
        // equivalent:
        // return m_multSum[i] - 1;
    }
     
    inline unsigned knotsUntilSpan(const size_t & i) const
    {
        return (ubegin()+i).multSum();
    }

    bool operator != (const gsKnotVector<T>& other) const
    {
        return ! ((*this)==other);
    }    

    inline T at (const size_t & i) const
    {
        return m_repKnots.at(i);
    }

    bool isUniform(T tol = 1e-9) const
    {
        const T df = *(ubegin() + 1) - *ubegin();
        for( uiterator uit = ubegin() + 1; uit != uend(); ++uit )
            if( math::abs(*uit - (*uit-(T)(1)) - df) > tol )
                return false;
        return true;
    }
     
    bool isOpen() const
    {
        const int dp1 = m_deg + 1;
        return (multFirst() == dp1 &&
                multLast () == dp1 );
        // equivalent
        //return ( ubegin  .multiplicity() == dp1 &&
        //         (--uend).multiplicity() == dp1 );
        // equivalent
        //return m_multSum.front() == dp1 && 
        //    m_multSum.back() - m_multSum.end()[-2] == dp1;
    }

    // \sa gsDomain
    virtual knotContainer breaks() const
    {
        return knotContainer(domainUBegin(), domainUEnd()+1);
    }

public: // others

    void getUniformRefinementKnots(mult_t knotsPerSpan, knotContainer& result, 
                                   mult_t mult = 1) const;

    std::string detail() const;

    T maxIntervalLength() const;

    T minIntervalLength() const;

    void degreeElevate(const short_t & i = 1);

    void degreeReduce(const short_t & i);

    std::vector<T> coarsen(index_t knotRemove = 1, index_t knotSkip = 1, mult_t  mul = -1);

public: // members

    // TODO remove!
    short_t m_deg;
};

template<typename T>
std::ostream& operator << (std::ostream& out, const gsKnotVector<T> KV )
{
    KV.print(out);
    return out;
}


#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsKnotVector(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo


// *****************************************************************
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsKnotVector.hpp)
#else
#ifdef gsKnotVector_EXPORT
#include GISMO_HPP_HEADER(gsKnotVector.hpp)
#undef  EXTERN_CLASS_TEMPLATE
#define EXTERN_CLASS_TEMPLATE CLASS_TEMPLATE_INST
#endif
namespace gismo
{
EXTERN_CLASS_TEMPLATE gsKnotVector<real_t>;
//EXTERN_CLASS_TEMPLATE internal::gsXml<gsKnotVector<real_t> >;
}
#endif
// *****************************************************************