gsTensorBSplineBasis.h

Go to the documentation of this file.

#pragma once

#include <gsTensor/gsTensorBasis.h>
#include <gsNurbs/gsBSplineBasis.h>
#include <gsNurbs/gsTensorBSpline.h>

namespace gismo
{

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

    typedef gsTensorBasis<d,T> Base;
    
    typedef gsBSplineBasis<T>  Family_t;

    typedef gsTensorBSplineBasis Self_t;

    typedef gsBSplineBasis<T> CoordinateBasis;
    typedef CoordinateBasis                  Basis_t;

    typedef T Scalar_t;

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

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

    typedef typename Base::iterator        iterator;
    typedef typename Base::const_iterator  const_iterator;

    typedef memory::shared_ptr< Self_t > Ptr;

    typedef memory::unique_ptr< Self_t > uPtr;

public:

    gsTensorBSplineBasis() : Base()
    {
        for(short_t i = 0; i!=d; ++i)
            this->m_bases[i] = new Basis_t();
    }

    void swap(gsTensorBSplineBasis & other)
    {
        this->Base::swap(static_cast<Base&>(other));
        std::swap(m_isPeriodic, other.m_isPeriodic);
    }
    
#if !EIGEN_HAS_RVALUE_REFERENCES
    gsTensorBSplineBasis & operator=(gsTensorBSplineBasis other)
    { gsTensorBSplineBasis::swap(other); return *this;}
// #else // defined implicitly
//     gsTensorBSplineBasis(gsTensorBSplineBasis && other) : Base(give(other)) { }
//     gsTensorBSplineBasis(const gsTensorBSplineBasis &  other) : Base(other) { }
//     gsTensorBSplineBasis & operator=(gsTensorBSplineBasis&& other)
//     { return (gsTensorBSplineBasis &)Base::operator=(give(other)); }
//     gsTensorBSplineBasis & operator=(const gsTensorBSplineBasis& other)
//     { return (gsTensorBSplineBasis &)Base::operator=(other); }
#endif
    
    template<typename U>
    gsTensorBSplineBasis( KnotVectorType KV1, gsKnotVector<U> KV2,
                          typename util::enable_if<d==2,U>::type * = NULL )
    : Base( new Basis_t(give(KV1)), new Basis_t(give(KV2)) )
    { m_isPeriodic = -1; }

    gsTensorBSplineBasis( KnotVectorType KV1, 
                          KnotVectorType KV2, 
                          KnotVectorType KV3 )
    : Base( new Basis_t(give(KV1)), new Basis_t(give(KV2)), new Basis_t(give(KV3)) )
    { m_isPeriodic = -1; }

    gsTensorBSplineBasis( KnotVectorType KV1, 
                          KnotVectorType KV2, 
                          KnotVectorType KV3,
                          KnotVectorType KV4)
    : Base( new Basis_t(give(KV1)), new Basis_t(give(KV2)),
            new Basis_t(give(KV3)), new Basis_t(give(KV4)) )
    { m_isPeriodic = -1; }

    explicit gsTensorBSplineBasis(std::vector<KnotVectorType> KV)
    { 
        GISMO_ENSURE(d == KV.size(), "Invalid number of knot-vectors given." );
        for(short_t i = 0; i!=d; ++i)
            this->m_bases[i] = new Basis_t( give(KV[i]) );
        m_isPeriodic = -1; 
    }

    gsTensorBSplineBasis( Basis_t* x,  Basis_t*  y) : Base(x,y) 
    { 
        GISMO_ENSURE(d==2,"Invalid constructor." );
        setIsPeriodic();
    }
    
    gsTensorBSplineBasis( Basis_t* x,  Basis_t* y, Basis_t* z ) : Base(x,y,z) 
    { 
        GISMO_ENSURE(d==3,"Invalid constructor." );
        setIsPeriodic();
        
    }
    
    gsTensorBSplineBasis( Basis_t* x,  Basis_t* y, Basis_t* z, Basis_t* w ) : Base(x,y,z,w) 
    { 
        GISMO_ENSURE(d==4,"Invalid constructor." );
        setIsPeriodic();
    }
    
    explicit gsTensorBSplineBasis(std::vector< gsBasis<T>*> & bb );
    
    explicit gsTensorBSplineBasis(std::vector< Basis_t*> & bb );

#ifdef __DOXYGEN__
    typename gsBSplineTraits<static_cast<short_t>(d-1),T>::Basis::uPtr boundaryBasis(boxSide const & s);
#endif
    GISMO_UPTR_FUNCTION_DEF(BoundaryBasisType, boundaryBasis, boxSide const &)
    {
        std::vector<gsBasis<T>*> rr;
        this->getComponentsForSide(n1,rr);
        return BoundaryBasisType::New(rr);
    }

    GISMO_CLONE_FUNCTION(gsTensorBSplineBasis)
    
    static Self_t * New(std::vector<gsBasis<T>*> & bb )
    { return new Self_t(bb); }

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

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

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

    KnotVectorType & knots (int i)
    { return Self_t::component(i).knots(); }

    const KnotVectorType & knots (int i) const 
    { return Self_t::component(i).knots(); }

    // knot \a k of direction \a i
    T knot(int i, int k) const 
    { return Self_t::component(i).knots()[k]; }


    const Basis_t & component(short_t dir) const
    {
        return static_cast<const Basis_t &>(Base::component(dir));
    }

    Basis_t & component(short_t dir)
    {
        return static_cast<Basis_t &>(Base::component(dir));
    }

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

    void active_cwise(const gsMatrix<T> & u, gsVector<index_t,d>& low,
                      gsVector<index_t,d>& upp ) const;

    std::ostream &print(std::ostream &os) const
    {
        os << "TensorBSplineBasis: dim=" << this->dim()<< ", size="<< this->size() <<".";
        if( m_isPeriodic != -1 )
            os << "Periodic in " << m_isPeriodic << "-th direction.\n";
        for ( short_t i = 0; i!=d; ++i )
            os << "\n  Direction "<< i <<": "<< Self_t::component(i).knots() <<" ";
        os << "\n";
        return os;
    }

    //
    // param other parent/reference mesh determining the smoothness at the inner knots.
    // param i number of k-refinement steps to perform

    void k_refine(Self_t & other, int const & i = 1)
    { 
        for (short_t j = 0; j < d; ++j)
            Self_t::component(j).refine_k(other.component(j), i);
    }

    void refine_p(int const & i = 1)
    {
        for (short_t j = 0; j < d; ++j)
            Self_t::component(j).refine_p(i);
    }
    
    void refine_h(int const & i = 1)
    {
        for (short_t j = 0; j < d; ++j)
            Self_t::component(j).refine_h(i);
    }

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


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

    void insertKnot(T knot, index_t dir, int mult=1)
    { this->knots(dir).insert( knot, mult); }

    void removeKnot(T knot, index_t dir, int mult=1)
    { this->knots(dir).remove( knot, mult); }

    void insertKnots(const std::vector< std::vector<T> >& refineKnots)
    {
        GISMO_ASSERT( refineKnots.size() == d, "refineKnots vector has wrong size" );
        for (short_t j = 0; j < d; ++j) // refine basis in each direction
            this->knots(j).insert(refineKnots[j]);
    }

    void refine( gsMatrix<T> const & boxes, int refExt = 0);

    GISMO_MAKE_GEOMETRY_NEW

    void reduceContinuity(int const & i = 1) 
    { 
        for (short_t j = 0; j < d; ++j)
            Self_t::component(j).reduceContinuity(i);
    }

    template <int _Rows>
    void elementSupport_into(const index_t i, gsMatrix<index_t, _Rows, 2> & result) const
    {
        //result.resize(d,2);
        gsMatrix<index_t> tmp_vec;
        const gsVector<index_t, d> ti = this->tensorIndex(i);

        for (short_t dim = 0; dim < d; ++dim)
        {
            const gsKnotVector<T> & kv = Self_t::component(dim).knots();
            kv.supportIndex_into(ti[dim], tmp_vec);
            result.row(dim).noalias() =
                tmp_vec.cwiseMax(0).cwiseMin(kv.numElements());
        }
    }

    gsMatrix<index_t, d, 2> elementSupport(const index_t & i) const
    {
        gsMatrix<index_t, d, 2> result(d, 2);
        elementSupport_into(i, result);
        return result;
    }

    template <int _Rows>
    void elementActive_into(const gsMatrix<index_t,_Rows,2> & box,
                             gsMatrix<index_t> & result) const
    {
        GISMO_ASSERT( box.rows() == static_cast<index_t>(d), "Invalid input box");
        gsMatrix<index_t,d,2> tmp;
        
        gsVector<index_t,d> str;
        this->stride_cwise(str);

        for (short_t dm = 0; dm != d; ++dm)
        {
            tmp(dm,0) = Self_t::component(dm).knots().lastKnotIndex (box(dm,0)) - this->degree(dm);
            tmp(dm,1) = Self_t::component(dm).knots().firstKnotIndex(box(dm,1)) - 1;
        }

        const gsVector<index_t,d> sz = tmp.col(1)- tmp.col(0) + gsVector<index_t,d>::Ones();

        result = gsVector<index_t>::LinSpaced(sz[0], tmp(0,0), tmp(0,1));
        for (short_t dm = 1; dm != d; ++dm)
            result = result.replicate(1,sz[dm]) + 
                gsVector<index_t>::Constant(result.rows(), str[dm] )
                * gsVector<index_t>::LinSpaced(sz[dm], tmp(dm,0), tmp(dm,1))
                .transpose();
    }

    template <int _Rows>
    gsMatrix<T> elementDom(const gsMatrix<index_t,_Rows,2> & box) const
    {
        GISMO_ASSERT( box.rows() == static_cast<index_t>(d), "Invalid input box");
        gsMatrix<T> rvo;
        rvo.resize(d,2);
        for (short_t dm = 0; dm != d; ++dm)
        {
            rvo(dm,0) = Self_t::component(dm).knots().uValue(box(dm,0));
            rvo(dm,1) = Self_t::component(dm).knots().uValue(box(dm,1));
        }
        return rvo;
    }
    
    inline bool isPeriodic() const { return m_isPeriodic != -1; }

    inline int periodicDirection() const { return m_isPeriodic; }

    inline void setPeriodic( const int dir )
    {
        Self_t::component(dir).setPeriodic();
        if( Self_t::component(dir).isPeriodic() ) // Only when succeeded when converting to periodic.
            m_isPeriodic = dir;
    }

    gsMatrix<T> perCoefs( const gsMatrix<T>& originalCoefs, short_t dir ) const
    {
        // Identify which coefficients to copy and where to copy them.
        std::vector<index_t> sourceSliceIndices;
        std::vector<index_t> targetSliceIndices;
        index_t numPeriodic = Self_t::component(dir).numCrossingFunctions();

        const index_t sz = this->size(dir) - numPeriodic;
        for( index_t i = 0; i < numPeriodic; i++ )
        {
            gsMatrix<index_t> currentSourceSlice = this->coefSlice(dir,i);
            gsMatrix<index_t> currentTargetSlice = this->coefSlice(dir, sz  + i );

            for( index_t j = 0; j != currentSourceSlice.size(); j++ )
            {
                sourceSliceIndices.push_back( static_cast<index_t>( currentSourceSlice(j) ) );
                targetSliceIndices.push_back( static_cast<index_t>( currentTargetSlice(j) ) );
            }
        }

        // Copy the chosen coefficients.
        gsMatrix<T> result = originalCoefs;
        for( size_t i = 0; i != sourceSliceIndices.size(); i++ )
        {
            //gsDebug << "source: " << sourceSliceIndices[i]  << "\n";
            //gsDebug << "target: " << targetSliceIndices[i]  << "\n";
            result.row( targetSliceIndices[ i ] ) = originalCoefs.row( sourceSliceIndices[ i ] );
        }

        return result;
    }

private:

    void setIsPeriodic()
    {
        m_isPeriodic = -1;
        for( short_t i = 0; i < this->dim(); i++ )
        {
            if( Self_t::component(i).isPeriodic() )
            {
                if( m_isPeriodic == -1 )
                    m_isPeriodic = i;
                else
                    gsWarn << "Cannot handle a basis that is periodic in more than one direction.\n";
            }
        }
    }

protected:

    short_t m_isPeriodic;

};

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsTensorBSplineBasis2(pybind11::module &m);
  void pybind11_init_gsTensorBSplineBasis3(pybind11::module &m);
  void pybind11_init_gsTensorBSplineBasis4(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo


// *****************************************************************
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsTensorBSplineBasis.hpp)
/*
#else
#ifdef gsTensorBSplineBasis_EXPORT
#include GISMO_HPP_HEADER(gsTensorBSplineBasis.hpp)
#undef  EXTERN_CLASS_TEMPLATE
#define EXTERN_CLASS_TEMPLATE CLASS_TEMPLATE_INST
#endif
namespace gismo
{
EXTERN_CLASS_TEMPLATE gsTensorBSplineBasis<2,real_t>;
EXTERN_CLASS_TEMPLATE gsTensorBSplineBasis<3,real_t>;
EXTERN_CLASS_TEMPLATE gsTensorBSplineBasis<4,real_t>;
}
*/
#endif
// *****************************************************************