gsBSplineBasis.h

Go to the documentation of this file.

#pragma once

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

#include <gsTensor/gsTensorBasis.h>
#include <gsTensor/gsTensorDomainIterator.h>
#include <gsTensor/gsTensorDomainBoundaryIterator.h>

#include <gsNurbs/gsKnotVector.h>

namespace gismo
{

template<short_t d, class T>
struct gsBSplineTraits
{
    typedef gsKnotVector<T> KnotVectorType;

    typedef gsTensorBSplineBasis<d,T> Basis;
    typedef gsTensorNurbsBasis<d,T>   RatBasis;
    typedef gsTensorBSpline<d,T>      Geometry;
    typedef gsTensorNurbs<d,T>        RatGeometry;
};
template<class T>
struct gsBSplineTraits<1,T>
{
    typedef gsKnotVector<T> KnotVectorType;
    typedef gsBSplineBasis<T>         Basis;
    typedef gsNurbsBasis<T>           RatBasis;
    typedef gsBSpline<T>              Geometry;
    typedef gsNurbs<T>                RatGeometry;
};
template<class T>
struct gsBSplineTraits<0,T>
{
    typedef gsKnotVector<T> KnotVectorType;
    typedef gsConstantBasis<T>                       Basis;
    typedef gsConstantBasis<T>                       RatBasis;
    typedef gsConstantFunction<T>                    Geometry;
    typedef gsConstantFunction<T>                    RatGeometry;
};

template<class T>
class gsTensorBSplineBasis<1,T> : public gsTensorBasis<1,T>
{
public:
    typedef gsKnotVector<T> KnotVectorType;

    typedef gsTensorBasis<1,T> Base;

    typedef gsBSplineBasis<T> Self_t;

    typedef gsTensorBSplineBasis TensorSelf_t;
    typedef gsTensorBSplineBasis<1,T> TensorSelf_tt;

    typedef T Scalar_t;

    typedef typename gsBSplineTraits<1,T>::Geometry GeometryType;

    typedef typename gsBSplineTraits<0,T>::Basis BoundaryBasisType;

    static const short_t Dim = 1;

    typedef memory::shared_ptr< Self_t > Ptr;

    typedef memory::unique_ptr< Self_t > uPtr;
  
public:

    // Look at gsBasis class for a description
    // Note: Specializing pointer type at return
    //GISMO_UPTR_FUNCTION_PURE(TensorSelf_tt, clone)
    private: virtual gsTensorBSplineBasis * clone_impl() const = 0;
    public: uPtr clone() const { return uPtr(dynamic_cast<Self_t*>(clone_impl())); }
    
    // gsTensorBSplineBasis( const Base & o)
    // { 
    //     const gsTensorBSplineBasis * a;
    //     if ( ( a = dynamic_cast<const gsTensorBSplineBasis *>( &o )) )
    //     {
    //         m_p        = a->degree() ;
    //         m_knots    = KnotVectorType( a->knots() );
    //         m_periodic = a->m_periodic;
    //     }
    //     else
    //         GISMO_ERROR("Cannot convert "<<o<<" to gsTensorBSplineBasis\n");
    // }

    static Self_t * New(std::vector<gsBasis<T>*> & bb );

    static Self_t * New(std::vector<Self_t*> & bb )
    { 
        return new Self_t(*bb.front());
    }

    static uPtr make(std::vector<gsBasis<T>*> & bb )
    {
        return uPtr(New(bb));
    }

    static uPtr make(std::vector<Self_t*> & bb )
    { 
        return uPtr(new Self_t(*bb.front()));
    }

    static uPtr make( const KnotVectorType & KV )
    {
        return uPtr(new Self_t(KV));
    }
    
    // Note: these casts can be dangerous
    // operator Self_t &() { return dynamic_cast<Self_t&>(*this);}
    // operator const Self_t &() const { return dynamic_cast<const Self_t&>(*this);}

public:

    void swap(gsTensorBSplineBasis& other)
    {
        std::swap(m_p, other.m_p);
        std::swap(m_periodic, other.m_periodic);
        m_knots.swap(other.m_knots);
    }

/* Virtual member functions required by the base class */

    // Look at gsBasis class for a description
    short_t domainDim() const { return Dim; }

    // Unhide/forward gsTensorBasis<1,T>::size(k), since the overload
    // with size() automatically hides it in this class
    using Base::size;

    // Look at gsBasis class for a description
    index_t size() const { return m_knots.size() - m_p - 1 - m_periodic; }

    // Look at gsBasis class for a description
    size_t numElements(boxSide const & s = 0) const { return (0==s?m_knots.numElements():1); }

    // Look at gsBasis class for a description
    size_t elementIndex(const gsVector<T> & u ) const;

    // Same as gsBasis::elementIndex but argument is a value instead of a vector
    size_t elementIndex(T u ) const;

    // Look at gsBasis class for a description
    gsMatrix<T> elementInSupportOf(index_t j) const;
    
    void elementSupport_into(const index_t i, gsMatrix<index_t,1,2>& result) const
    {
        gsMatrix<index_t> tmp_vec;
        m_knots.supportIndex_into(i, tmp_vec);
        result = tmp_vec.cwiseMax(0).cwiseMin(m_knots.numElements());
    }

    template <int _Rows>
    gsMatrix<T> elementDom(const gsMatrix<index_t,_Rows,2> & box) const
    {
        gsMatrix<T> rvo(1,2);
        rvo.at(0) = m_knots.uValue(box.at(0));
        rvo.at(1) = m_knots.uValue(box.at(1));
        return rvo;
    }

    // Look at gsBasis class for a description
    const TensorSelf_t & component(short_t i) const = 0;

    // Look at gsBasis class for a description
    TensorSelf_t & component(short_t i) = 0;

    void anchors_into(gsMatrix<T> & result) const 
    { 
        m_knots.greville_into(result); 
    }

    void anchor_into(index_t i, gsMatrix<T> & result) const
    { 
        result.resize(1,1);
        result(0,0) = m_knots.greville(i);
    }

    // Look at gsBasis class for a description
    void connectivity(const gsMatrix<T> & nodes, 
                      gsMesh<T> & mesh) const;

    // Look at gsBasis class for a description
    void active_into(const gsMatrix<T> & u, gsMatrix<index_t>& result) const;

    // Look at gsBasis class for a description
    bool isActive(const index_t i, const gsVector<T> & u) const;

    // Look at gsBasis class for a description
    gsMatrix<index_t> allBoundary( ) const ;

    // Look at gsBasis class for a description
    gsMatrix<index_t> boundaryOffset(boxSide const & s,index_t offset) const;

#ifdef __DOXYGEN__
    typename BoundaryBasisType::uPtr boundaryBasis(boxSide const & s);
#endif
    GISMO_UPTR_FUNCTION_DEC(gsConstantBasis<T>, boundaryBasis, boxSide const &)

    // Look at gsBasis class for a description
    gsMatrix<T> support() const ;

    // Look at gsBasis class for a description
    gsMatrix<T> support(const index_t & i ) const ;

    index_t twin(index_t i) const ;

    // Look at gsBasis class for a description
    virtual void eval_into(const gsMatrix<T> & u, gsMatrix<T>& result) const;

    // Look at gsBasis class for a description
    virtual void evalSingle_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result) const;

    // Look at gsBasis class for a description
    virtual void evalFunc_into(const gsMatrix<T> & u, const gsMatrix<T> & coefs, gsMatrix<T>& result) const;

    // Look at gsBasis class for a description
    void deriv_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const ;

    // Look at gsBasis class for a description
    void derivSingle_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const ;

    // Look at gsBasis class for a description
    void deriv_into(const gsMatrix<T> & u, const gsMatrix<T> & coefs, gsMatrix<T>& result ) const ;

    // Look at gsBasis class for a description
    void deriv2_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const ;

    // Look at gsBasis class for a description
    void deriv2Single_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const ;

    // Look at gsBasis class for a description
    void deriv2_into(const gsMatrix<T> & u, const gsMatrix<T> & coefs, gsMatrix<T>& result ) const ;

    // Look at gsBasis class for a description
    gsMatrix<T> laplacian(const gsMatrix<T> & u ) const ;

    // Look at gsBasis class for a description
    typename gsBasis<T>::uPtr tensorize(const gsBasis<T> & other) const;
    
    bool check() const
    { 
        if ( m_periodic > 0 )
        {
            // Periodicity check wrt knot values
            return ( 
                m_knots.degree()      == m_p &&
                (int)m_knots.size()   >  2*m_p+1
                );           
        }
        else
        {
            return ( 
                m_knots.degree()      == m_p &&
                (int)m_knots.size()   >  2*m_p+1
                ); 
        }
    }
  
    std::ostream &print(std::ostream &os) const;

    std::string detail() const;

    // Look at gsBasis class for a description
    virtual void evalDerSingle_into(index_t i, const gsMatrix<T> & u,
                                    int n, gsMatrix<T>& result) const;

    // Look at gsBasis class for a description
    virtual void evalAllDers_into(const gsMatrix<T> & u, int n,
                                  std::vector<gsMatrix<T> >& result,
                                  bool sameElement = false) const;

    // Look at gsBasis class for a description
    virtual void evalAllDersSingle_into(index_t i, const gsMatrix<T> & u,
                                        int n, gsMatrix<T>& result) const;

    // Look at gsBasis class for a description
    short_t degree(short_t i) const
    {
        GISMO_UNUSED(i);
        GISMO_ASSERT(i==0,"Asked for degree(i) in 1D basis.");
        return m_p; 
    }

    short_t degree() const {return m_p;}

    // Look at gsBasis class for a description
    short_t maxDegree()   const { return m_p; }

    // Look at gsBasis class for a description
    short_t minDegree()   const { return m_p; }

    // Look at gsBasis class for a description
    short_t totalDegree() const { return m_p; }
 
    inline unsigned order() const { return m_p+1; }

    inline bool inDomain(T const & pp) const 
    { return ( (pp >= *(m_knots.begin()+m_p)) &&  (pp <= *(m_knots.end()-m_p-1) ) ); }

    T domainStart() const { return *(m_knots.begin()+m_p); }

    T domainEnd() const { return *(m_knots.end()-m_p-1); }

    T _activeLength() const { return domainEnd() - domainStart(); }

    inline index_t firstActive(T u) const {
        return ( inDomain(u) ? (m_knots.iFind(u)-m_knots.begin()) - m_p : 0 );
    }

    // Number of active functions at any point of the domain
    inline index_t numActive() const { return m_p + 1; }

    // Look at gsBasis class for a description
    gsDomain<T> * domain() const { return const_cast<KnotVectorType *>(&m_knots); }

    const KnotVectorType & knots (int i  = 0) const 
    { 
        GISMO_ENSURE(i==0, "Invalid knots requested");
        return m_knots;
    }
    KnotVectorType & knots (int i  = 0)
    {
        GISMO_ENSURE(i==0, "Invalid knots requested");
        return m_knots;
    }

    T knot(index_t const & i) const { return m_knots[i];}

    void insertKnot(T knot, int mult=1)
    { m_knots.insert( knot, mult); }

    void removeKnot(T knot, int mult=1)
    { m_knots.remove( knot, mult); }

    // compatibility with tensor-bsplines
    void insertKnots(const std::vector< std::vector<T> >& refineKnots)
    {
        GISMO_ASSERT( refineKnots.size() == 1, "refineKnots vector has wrong size" );
        this->knots().insert(refineKnots.front());
    }

    // Look at gsBasis class for a description
    void refineElements(std::vector<index_t> const & elements)
    { m_knots.refineSpans(elements); }

    // Look at gsBasis class for a description
    void uniformRefine(int numKnots = 1, int mul=1, int dir=-1)
    { GISMO_UNUSED(dir); m_knots.uniformRefine(numKnots,mul); }

    // Look at gsBasis class for a description
    void uniformRefine_withCoefs(gsMatrix<T>& coefs, int numKnots = 1, int mul=1, int dir=-1);

    // Look at gsBasis class for a description
    void uniformRefine_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, int numKnots = 1, int mul=1);
    
    // Look at gsBasis class for a description
    void uniformCoarsen(int numKnots = 1)
    { m_knots.coarsen(numKnots); }

    // Look at gsBasis class for a description
    void uniformCoarsen_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, int numKnots = 1);

    void refine_withCoefs(gsMatrix<T>& coefs, const std::vector<T>& knots);

    // compatibility with tensor-bsplines
    void refine_withCoefs(gsMatrix<T> & coefs, 
                          const std::vector< std::vector<T> >& refineKnots)
    {
        refine_withCoefs(coefs,refineKnots.front() );
    }

    void refine_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, const std::vector<T>& knots);

    void refine_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, 
                             const std::vector<std::vector<T> >& knots)
    {
        GISMO_ASSERT(knots.size()==1, "the knots you want to insert do not have the right dimension");
        refine_withTransfer(transfer,knots.front());
        //GISMO_NO_IMPLEMENTATION
    }

    void refine_k(const TensorSelf_t & other, int const & i = 1)
    { 
        GISMO_ASSERT( m_p >= other.m_p, "Degree of other knot-vector should be lower.");
        //for (typename std::vector<T>::iterator it = 
        //         m_knots.begin(); it != m_knots.end(); ++it)
        //    GISMO_ASSERT( has(*it), "Knot "<< *it<<" is not in the knot vector.");

        // grab unique knots
        const std::vector<T> knots = other.m_knots.unique();

        // Increase the degree without adjusting any knot
        m_p += i; 
        m_knots.set_degree(m_p);
        // Adjust (reduce) smoothness to satisfy initial constraint knots
        m_knots.insert(knots,i);
    }

    void refine_p(short_t const & i = 1)
    { degreeElevate(i); }

    void refine_h(short_t const & i = 1)
    { uniformRefine(i); }
  
    void degreeElevate(short_t const & i = 1, short_t const dir = -1)
    {
        GISMO_UNUSED(dir);
        GISMO_ASSERT( dir == -1 || dir == 0, "Invalid direction");
        m_p+=i; m_knots.degreeElevate(i); 
    }

    // Look at gsBasis for documentation
    void degreeReduce (short_t const & i = 1, short_t const dir = -1)
    {
        GISMO_UNUSED(dir);
        GISMO_ASSERT( dir == -1 || dir == 0, "Invalid direction");
        GISMO_ASSERT( i<=m_p, "Cannot reduce degree to negative");
        m_p-=i; m_knots.degreeReduce(i);
        //m_periodic =
    }

    // Look at gsBasis for documentation
    void degreeIncrease(short_t const & i = 1, short_t const dir = -1)
    {
        GISMO_UNUSED(dir);
        GISMO_ASSERT( dir == -1 || dir == 0, "Invalid direction");
        m_p+=i; m_knots.degreeIncrease(i);
    }

    // Look at gsBasis for documentation
    void degreeDecrease(short_t const & i = 1, short_t const dir = -1)
    {
        GISMO_UNUSED(dir);
        GISMO_ASSERT( dir == -1 || dir == 0, "Invalid direction");
        m_p-=i; m_knots.degreeDecrease(i);
    }

    void elevateContinuity(int const & i = 1) 
    { 
        GISMO_ASSERT( i>=0 && ( m_knots.size()>static_cast<size_t>(2*(m_p+1)) || i<=m_p ), 
                      "Cannot achieve continuity less than C^{-1} at interior knots.");
        m_knots.reduceMultiplicity(i);
    }

    void reduceContinuity(int const & i = 1) 
    { 
        GISMO_ASSERT( i>=0 && ( m_knots.size()>static_cast<size_t>(2*(m_p+1)) || i<=m_p ), 
                      "Cannot achieve continuity less than C^{-1} at interior knots.");
        // TODO check: max interior mult + i <= m_p+1
        m_knots.increaseMultiplicity(i);
    }

    bool isPeriodic() const
    {
        return (m_periodic > 0);
    }

    // Look at gsBasis class for a description
    index_t functionAtCorner(boxCorner const & c) const;

    int numCrossingFunctions () const
    {
        return m_periodic;
    }

    bool isClamped() const
    {
        if( m_knots[0] != m_knots[m_p])
            return false;

        else if( m_knots[m_knots.size() - m_p -1] != m_knots[m_knots.size()-1])
            return false;

        else
            return true;
    }

    void setPeriodic(bool flag = true)
    {
        if ( flag )
            _convertToPeriodic();
        else
            m_periodic = 0;
    }

    // Compatible with tensor B-spline basis
    void setPeriodic(int dir)
    {
        GISMO_UNUSED(dir);
        GISMO_ASSERT(dir==0, "Invalid direction");
            _convertToPeriodic();
    }

    int borderKnotMult() const;

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

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

    void enforceOuterKnotsPeriodic();

    // Look at gsBasis class for a description
    void reverse() { m_knots.reverse(); }

    void matchWith(const boundaryInterface &, const gsBasis<T> &,
                   gsMatrix<index_t> &, gsMatrix<index_t> &, index_t) const { }

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

protected:

    void _convertToPeriodic();

    void _stretchEndKnots();

public:

    gsMatrix<T> perCoefs( const gsMatrix<T>& coefs ) const
    {
        gsMatrix<T> per_coefs = coefs;
        per_coefs.bottomRows( m_periodic ) = coefs.topRows( m_periodic );
        return per_coefs;
    }

    gsMatrix<T> perCoefs(const gsMatrix<T>& coefs, int dir) const
    {
        GISMO_UNUSED(dir);
        GISMO_ASSERT(dir==0, "Error");
        return perCoefs(coefs);
    }

    void expandCoefs(gsMatrix<T> & coefs) const
    {
        const index_t sz = coefs.rows();
        coefs.conservativeResize(sz+m_periodic, gsEigen::NoChange);
        coefs.bottomRows( m_periodic ) = coefs.topRows( m_periodic );
    }

    void trimCoefs(gsMatrix<T> & coefs) const
    {
        const index_t sz = coefs.rows();
        coefs.conservativeResize(sz-m_periodic, gsEigen::NoChange);
    }

    int trueSize() const
    { return this->size() + m_periodic; }
 
// Data members
protected:

    short_t m_p;

    KnotVectorType m_knots;
    
    int m_periodic;

    /*/// @brief Multiplicity of the p+1st knot from the beginning and from the end.
      int m_bordKnotMulti;*/

}; // class gsTensorBSplineBasis<1>


//Using C++11 alias:
// template<class T, class KnotVectorType>
// using gsBSplineBasis = gsTensorBSplineBasis<1,T>

template<class T>
class gsBSplineBasis : public gsTensorBSplineBasis<1,T>
{
public:
    typedef memory::shared_ptr< gsBSplineBasis > Ptr;
    typedef memory::unique_ptr< gsBSplineBasis > uPtr;

    typedef gsKnotVector<T> KnotVectorType;
    typedef gsTensorBSplineBasis<1,T> Base;
    typedef gsBSplineBasis<T> Self_t;

    typedef typename gsBSplineTraits<1,T>::Geometry GeometryType;

    typedef typename gsBSplineTraits<0,T>::Basis BoundaryBasisType;

public:

    // Default empty constructor
    explicit gsBSplineBasis(const bool periodic = false )
    { 
        m_p = 0;
        m_knots.initClamped(0);
        m_periodic = 0;

        if( periodic )
            this->_convertToPeriodic();

        if( ! this->check() )
            gsWarn << "Warning: Inconsistent "<< *this<< "\n";
    }

    explicit gsBSplineBasis(KnotVectorType KV, const bool periodic = false)
    { 
        m_p        = KV.degree();
        m_knots.swap(KV);
        m_periodic = 0;

        if( periodic )
            this->_convertToPeriodic();

        if( ! this->check() )
            gsWarn << "Warning: Insconsistent "<< *this<< "\n";
    }

    explicit gsBSplineBasis(std::vector<KnotVectorType> KV)
    { 
        GISMO_ASSERT(1 == KV.size(), "Expecting a single knotvector." );

        m_p        = KV.front().degree();
        m_knots.swap(KV.front());
        m_periodic = 0;
    }

    gsBSplineBasis(const T u0, const T u1, const unsigned interior, 
                   const int degree, const unsigned mult_interior=1,
                   const bool periodic = false )
    { 
        m_p = degree;
        m_knots.initUniform(u0, u1, interior, m_p+1, mult_interior, m_p);
        m_periodic = 0;

        if( periodic )
            this->_convertToPeriodic();
        
        if( ! this->check() )
            gsWarn << "Warning: Insconsistent "<< *this<< "\n";
    }

/*
    gsBSplineBasis( const gsBSplineBasis & o)
    : Base(o)
    { }
*/

    // Look at gsBasis class for a description
    GISMO_CLONE_FUNCTION(gsBSplineBasis)

    // Look at gsBasis class for a description
    Self_t & component(short_t i);
    
    // Look at gsBasis class for a description
    const Self_t & component(short_t i) const;

    memory::unique_ptr<gsGeometry<T> > makeGeometry( gsMatrix<T> coefs ) const;
            
private:
    
    using Base::m_p;
    using Base::m_knots;
    using Base::m_periodic;
};

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsBSplineBasis(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo


// *****************************************************************
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsBSplineBasis.hpp)
#else
#ifdef gsBSplineBasis_EXPORT
#include GISMO_HPP_HEADER(gsBSplineBasis.hpp) // needed on thebrain?
#undef  EXTERN_CLASS_TEMPLATE
#define EXTERN_CLASS_TEMPLATE CLASS_TEMPLATE_INST
#endif
namespace gismo
{
template<class T> const short_t gsTensorBSplineBasis<1,T>::Dim; //-O3 (SLE11) fix
EXTERN_CLASS_TEMPLATE gsTensorBSplineBasis<1,real_t>;
EXTERN_CLASS_TEMPLATE gsBSplineBasis<real_t>;
}
#endif
// *****************************************************************