gsTensorBSpline.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsGeometry.h>
#include <gsNurbs/gsTensorBSplineBasis.h>

namespace gismo
{

template <short_t d, class T>
void constructCoefsForSlice(index_t dir_fixed, const index_t index,
                            const gsMatrix<T> & fullCoefs,
                            const gsVector<index_t, d> & sizes,
                            gsMatrix<T>& result);

template<short_t d, class T>
class gsTensorBSpline GISMO_FINAL : public gsGeoTraits<d,T>::GeometryBase
{
public: 
    typedef gsKnotVector<T> KnotVectorType;

    typedef T Scalar_t;

    typedef typename gsBSplineTraits<d,T>::Basis Basis;

    typedef typename gsGeoTraits<d,T>::GeometryBase Base;

    typedef gsBSplineBasis<T>  Family_t;

    typedef memory::shared_ptr< gsTensorBSpline > Ptr;

    typedef memory::unique_ptr< gsTensorBSpline > uPtr;

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

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

public:

    gsTensorBSpline() { }

    using gsGeometry<T>::swap;

#if !EIGEN_HAS_RVALUE_REFERENCES
    gsTensorBSpline & operator=(gsTensorBSpline other)
    { this->swap(other); return *this;}
#endif

    // Construct B-Spline by basis functions and coefficient matrix
    //gsTensorBSpline( const gsConstantBasis<T> & basis, const gsMatrix<T> & coefs )
    //{ }

    gsTensorBSpline( const Basis & basis, gsMatrix<T> coefs )
    : Base( basis, give(coefs) )
    { }
        
    template<typename U>
    gsTensorBSpline( KnotVectorType KV1, gsKnotVector<U> KV2, 
                     gsMatrix<T> tcoefs,
                     typename util::enable_if<d==2,U>::type * = NULL)
    {
        GISMO_ASSERT(d==2, "Wrong dimension: tried to make a "
                     << d<<"D tensor B-spline using 2 knot-vectors.");

        std::vector<Family_t*> cbases;
        cbases.push_back(new gsBSplineBasis<T>(give(KV1)) );
        cbases.push_back(new gsBSplineBasis<T>(give(KV2)) );
        Basis * tbasis = Basis::New(cbases); //d==2
        
        GISMO_ASSERT(tbasis->size()==tcoefs.rows(),
                     "Coefficient matrix for the tensor B-spline does not have "
                     "the expected number of control points (rows)." );
        
        this->m_basis = tbasis;
        this->m_coefs = give(tcoefs);
    }
    
    gsTensorBSpline(gsMatrix<T> const & corner,
                    KnotVectorType KV1,
                    KnotVectorType KV2);

    gsTensorBSpline( KnotVectorType KV1,
                     KnotVectorType KV2,
                     KnotVectorType KV3,
                     gsMatrix<T> tcoefs )
    {
        GISMO_ASSERT(d==3, "Wrong dimension: tried to make a "
                     << d<<"D tensor B-spline using 3 knot-vectors.");

        std::vector<Family_t*> cbases;
        cbases.push_back(new gsBSplineBasis<T>(give(KV1)) );
        cbases.push_back(new gsBSplineBasis<T>(give(KV2)) );
        cbases.push_back(new gsBSplineBasis<T>(give(KV3)) );
        Basis * tbasis = Basis::New(cbases); //d==3
        
        GISMO_ASSERT(tbasis->size()==tcoefs.rows(),
                     "Coefficient matrix for the tensor B-spline does not have "
                     "the expected number of control points (rows)." );
        
        this->m_basis = tbasis;
        this->m_coefs = give(tcoefs);
    }

    gsTensorBSpline( KnotVectorType KV1,
                     KnotVectorType KV2,
                     KnotVectorType KV3,
                     KnotVectorType KV4,
                     gsMatrix<T> tcoefs )
    {
        GISMO_ASSERT(d==4, "Wrong dimension: tried to make a "
                     << d<<"D tensor B-spline using 4 knot-vectors.");

        std::vector<Family_t*> cbases;
        cbases.reserve(4);
        cbases.push_back(new gsBSplineBasis<T>(give(KV1)) );
        cbases.push_back(new gsBSplineBasis<T>(give(KV2)) );
        cbases.push_back(new gsBSplineBasis<T>(give(KV3)) );
        cbases.push_back(new gsBSplineBasis<T>(give(KV4)) );
        Basis * tbasis = Basis::New(cbases); //d==4
        
        GISMO_ASSERT(tbasis->size()==tcoefs.rows(),
                     "Coefficient matrix for the tensor B-spline does not have "
                     "the expected number of control points (rows)." );
        
        this->m_basis = tbasis;
        this->m_coefs = give(tcoefs);
    }
    
    GISMO_CLONE_FUNCTION(gsTensorBSpline)

    GISMO_BASIS_ACCESSORS

public:

    // Look at gsGeometry class for a description
    void degreeElevate(short_t const i = 1, short_t const dir = -1);

    void insertKnot( T knot, int dir, int i = 1);

    KnotVectorType & knots(const int i) { return this->basis().knots(i); }

    const KnotVectorType & knots(const int i) const { return this->basis().knots(i); }

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

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

    /*** Additional members for tensor B-Splines ***/

    short_t degree(const unsigned & i) const
    { return this->basis().component(i).degree(); }

    void reverse(unsigned k);

    void swapDirections(const unsigned i, const unsigned j);

    void toggleOrientation();

    bool isPatchCorner(gsMatrix<T> const &v, T tol = 1e-3) const;

    void findCorner(const gsMatrix<T>   & v,
                    gsVector<index_t,d> & curr,
                    T tol = 1e-3);

    void setOriginCorner(gsMatrix<T> const &v);

    void setFurthestCorner(gsMatrix<T> const &v);
    
    inline void setPeriodic( int dir )
    {
        this->m_coefs = this->basis().perCoefs( this->m_coefs, dir );
        this->basis().setPeriodic( dir );
    }
    
    typename gsGeometry<T>::uPtr localRep(const gsMatrix<T> & u) const;

public:

    void slice(index_t dir_fixed,T par,BoundaryGeometryType & result) const;

    std::vector<gsGeometry<T>* > uniformSplit(index_t dir = -1) const;

    std::vector<gsGeometry<T>* > splitAtMult(index_t minMult = 1, index_t dir = -1) const;

    void splitAt( index_t dir,T xi, gsTensorBSpline<d,T>& left,  gsTensorBSpline<d,T>& right) const;

    typename gsGeometry<T>::uPtr iface(const boundaryInterface & bi,
                                       const gsGeometry<T> & other) const;

protected:
    // TODO Check function
    // check function: check the coefficient number, degree, knot vector ...

    using Base::m_basis;
    using Base::m_coefs;

}; // class gsTensorBSpline

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsTensorBSpline2(pybind11::module &m);
  void pybind11_init_gsTensorBSpline3(pybind11::module &m);
  void pybind11_init_gsTensorBSpline4(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo

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