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_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, gsMatrix<T> const & boxes)
    {
        this->refine_withTransfer(transfer,this->_boxToKnots(boxes));
    }
 
    void refine_withCoefs(gsMatrix<T> & coefs,const std::vector< std::vector<T> >& refineKnots);
    void refine_withCoefs(gsMatrix<T> & coefs, gsMatrix<T> const & boxes)
    {
        this->refine_withCoefs(coefs,this->_boxToKnots(boxes));
    }
 
    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:
 
    std::vector<std::vector<T>> _boxToKnots(gsMatrix<T> const & boxes);
 
    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
// *****************************************************************