gsRationalTHBSplineBasis.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsCore/gsRationalBasis.h>
#include <gsHSplines/gsRationalTHBSpline.h>
#include <gsCore/gsForwardDeclarations.h>
 
 
namespace gismo
{
 
template<short_t d, class T>
class gsRationalTHBSplineBasis : public gsRationalBasis<gsTHBSplineBasis<d,T>>
{
 
public: 
 
    typedef gsRationalBasis<gsTHBSplineBasis<d,T>> Base;
 
    typedef gsTHBSplineBasis<d,T> Src_t;
 
    typedef Src_t Basis_t;
 
    typedef T Scalar_t;
 
    typedef typename gsHTensorBasis<d,T>::tensorBasis tensorBasis;
 
    typedef gsRationalTHBSpline<d,T> GeometryType;
 
    typedef gsRationalTHBSplineBasis< (1<d?d-1:d),T> BoundaryBasisType;
 
    typedef memory::shared_ptr< gsRationalTHBSplineBasis > Ptr;
 
    typedef memory::unique_ptr< gsRationalTHBSplineBasis > uPtr;
    
    // typedef typename Base::iterator iterator;
    // typedef typename Base::const_iterator const_iterator;
 
public:
    // Constructors forwarded from the base class
    gsRationalTHBSplineBasis() : Base() { };
 
    gsRationalTHBSplineBasis(const Src_t& basis ) : Base(basis) { }
 
    gsRationalTHBSplineBasis( Src_t* basis, gsMatrix<T> w ) : Base(basis, give(w)) { }
 
    gsRationalTHBSplineBasis( gsConstantBasis<T>* basis, gsMatrix<T>)
    {GISMO_ERROR("!!"); }
 
    gsRationalTHBSplineBasis(const gsRationalTHBSplineBasis & o) : Base(o) { }
 
    GISMO_CLONE_FUNCTION(gsRationalTHBSplineBasis)
  
    GISMO_MAKE_GEOMETRY_NEW
 
    std::ostream &print(std::ostream &os) const
    {
        os << "Rational THB-Spline basis: dim=" << this->dim()<< ", size="<< this->size() << ".";
        // for ( unsigned i = 0; i!=d; ++i )
        //     os << "\n  Direction "<< i <<": "<< this->m_src->component(i).knots() <<" ";
        os << "\n";
        return os;
    }
 
    GISMO_UPTR_FUNCTION_DEF(BoundaryBasisType, boundaryBasis, boxSide const &)
    {
        typename Src_t::BoundaryBasisType::uPtr bb = m_src->boundaryBasis(n1);
        gsMatrix<index_t> ind = m_src->boundary(n1);
 
        gsMatrix<T> ww( ind.size(),1);
        for ( index_t i=0; i<ind.size(); ++i)
            ww(i,0) = m_weights( (ind)(i,0), 0);
 
        return new BoundaryBasisType(bb.release(), give(ww));// note: constructor consumes the pointer
    }
public:
    void refine_withCoefs(gsMatrix<T> & coefs, gsMatrix<T> const & boxes)
    {
        auto tmp = m_src->clone();
        coefs = m_weights.asDiagonal() * coefs;
        tmp->refine_withCoefs(coefs, boxes);
        m_src->refine_withCoefs(m_weights, boxes);
        coefs.array().colwise() /= m_weights.col(0).array();
    }
 
    // void size_cwise(gsVector<index_t,d> & result) const
    // {
    //     // call the function of the underlying basis
    //     m_src->size_cwise(result);
    // }
 
    // /// Returns the strides for all dimensions
    // void stride_cwise(gsVector<index_t,d> & result) const
    // {
    //     // call the function of the underlying basis
    //     m_src->stride_cwise(result);
    // }
 
    // void swapDirections(const unsigned i, const unsigned j)
    // {
    //     gsVector<index_t, d> sz;
    //     size_cwise(sz);
 
    //     // First swap the weights
    //     swapTensorDirection(i, j, sz, m_weights);
 
    //     // Then swap the basis components
    //     m_src->swapDirections(i, j);
    // }
 
    // void uniformRefine_withCoefs(gsMatrix<T>& coefs, int numKnots=1, int mul=1, int dir=-1)
    // {
    //     GISMO_ASSERT( coefs.rows() == this->size() && m_weights.rows() == this->size(),
    //                   "Invalid dimensions" );
 
    //     gsSparseMatrix<T, RowMajor> transfer;
    //     if (dir==-1)
    //     {
    //         m_src->uniformRefine_withTransfer(transfer, numKnots, mul);
 
    //         coefs     = transfer * ( m_weights.asDiagonal() * coefs);
    //         m_weights = transfer * m_weights;
    //         // Alternative way
    //         // gsBasis<T> * tmp = m_src->clone();
    //         // tmp->uniformRefine_withCoefs(coefs, numKnots);
    //         // delete tmp;
    //         // m_src->uniformRefine_withCoefs(m_weights, numKnots);
 
    //         // back to affine coefs
    //         coefs.array().colwise() /= m_weights.col(0).array();
    //         // equiv:
    //         // for (int i = 0; i < coefs.rows(); ++i)
    //         //    coefs.row(i) /= m_weights.at(i);
    //     }
    //     else
    //     {
    //         GISMO_ASSERT( dir >= 0 && static_cast<unsigned>(dir) < d,
    //                       "Invalid basis component "<< dir <<" requested for degree elevation" );
 
    //         gsVector<index_t,d> sz;
    //         m_src->size_cwise(sz);
    //         m_src->component(dir).uniformRefine_withTransfer( transfer, numKnots, mul );
 
    //         const index_t coefs_cols = coefs.cols();
    //         const index_t weights_cols = m_weights.cols();
 
    //         coefs = m_weights.asDiagonal() * coefs; //<<<<-----this goes wrong!!
    //         swapTensorDirection(0, dir, sz, coefs);
    //         coefs.resize( sz[0], coefs_cols * sz.template tail<static_cast<short_t>(d-1)>().prod() );
    //         coefs     = transfer * coefs;
 
    //         swapTensorDirection(0, dir, sz, m_weights);
    //         m_weights.resize( sz[0], weights_cols * sz.template tail<static_cast<short_t>(d-1)>().prod() );
    //         m_weights = transfer * m_weights;
 
    //         sz[0] = coefs.rows();
 
    //         coefs.resize( sz.prod(), coefs_cols );
    //         m_weights.resize( sz.prod(), weights_cols );
    //         swapTensorDirection(0, dir, sz, coefs);
    //         swapTensorDirection(0, dir, sz, m_weights);
 
    //         coefs.array().colwise() /= m_weights.col(0).array();
    //     }
    // }
 
#ifdef __DOXYGEN__
    typename BoundaryBasisType::uPtr boundaryBasis(boxSide const & s);
#endif
 
    // GISMO_UPTR_FUNCTION_DEF(BoundaryBasisType, boundaryBasis, boxSide const &)
    // {
    //     typename Src_t::BoundaryBasisType::uPtr bb = m_src->boundaryBasis(n1);
    //     gsMatrix<index_t> ind = m_src->boundary(n1);
        
    //     gsMatrix<T> ww( ind.size(),1);
    //     for ( index_t i=0; i<ind.size(); ++i)
    //         ww(i,0) = m_weights( (ind)(i,0), 0);
        
    //     return new BoundaryBasisType(bb.release(), give(ww));// note: constructor consumes the pointer
    // }
 
    // void matchWith(const boundaryInterface & bi, const gsBasis<T> & other,
    //                gsMatrix<index_t> & bndThis, gsMatrix<index_t> & bndOther) const
    // {
    //     this->matchWith(bi,other,bndThis,bndOther,0);
    // }
 
    // // see gsBasis for documentation
    // void matchWith(const boundaryInterface & bi, const gsBasis<T> & other,
    //                gsMatrix<index_t> & bndThis, gsMatrix<index_t> & bndOther, index_t offset) const
    // {
    //     if ( const gsRationalTHBSplineBasis<d,T> * _other = dynamic_cast<const gsRationalTHBSplineBasis<d,T> *>(&other) )
    //         m_src->matchWith(bi,_other->source(),bndThis,bndOther,offset);
    //     else if ( const gsTensorBasis<d,T> * __other = dynamic_cast<const gsTensorBasis<d,T> *>(&other) )
    //         m_src->matchWith(bi,*__other,bndThis,bndOther,offset);
    //     else
    //         gsWarn<<"Cannot match with "<<other<<"\n";
    // }
 
 
protected:
    using Base::m_src;
    using Base::m_weights;
 
};
 
 
} // namespace gismo