gsTHBSplineBasis.h

Go to the documentation of this file.

#pragma once

#include <gsHSplines/gsHTensorBasis.h>
#include <gsHSplines/gsTHBSpline.h>


namespace gismo
{  

template<short_t d, class T>
class gsTHBSplineBasis : public gsHTensorBasis<d,T>
{
public:
    typedef gsTHBSpline<d,T> GeometryType;
    
    typedef typename gsHTensorBasis<d,T>::CMatrix CMatrix;

    typedef typename gsHTensorBasis<d,T>::cmatIterator cmatIterator;
    
    typedef typename gsHTensorBasis<d,T>::tensorBasis tensorBasis;

    typedef typename gsHTensorBasis<d,T>::point point;

    typedef memory::shared_ptr< gsTHBSplineBasis > Ptr;

    typedef memory::unique_ptr< gsTHBSplineBasis > uPtr;

    typedef typename
    util::conditional<d==1, gsConstantBasis<T>, gsTHBSplineBasis<static_cast<short_t>(d-1),T> >::type BoundaryBasisType;

    using gsHTensorBasis<d, T>::flatTensorIndexOf;
    using gsHTensorBasis<d, T>::m_manualLevels;
    using gsHTensorBasis<d, T>::_diadicIndexToKnotIndex;

    // polygon lines in parameter domain
    // the stucture is [levels [ line [ segments [ x y z w ] ] ] ],
    // where x y z w describes segment from (x, y) to (z, w)
    typedef typename std::vector< std::vector< std::vector< std::vector<T> > > > Polylines;

    // axis aligned bounding boxes in parameter domain
    // the structure is [levels [ boxes [ low_x low_y upp_x upp_y] ] ]
    // where the box is define by lower left corner (low_x, low_y) and upper right 
    // corner (upp_x, upp_y)
    typedef typename std::vector< std::vector< std::vector<index_t> > > AxisAlignedBoundingBox;

    //typedef typename std::vector< gsAabb<T> > AxisAlignedBoundingBox;

    // trimming curves in parameter domain
    // the stucture is [level [connected componenet [ line [ segments [ x y z w ] ] ] ] ],
    // where x y z w describes segment from (x, y) to (z, w), and first line 
    // indicates outer loop and next line describes holes
    typedef typename std::vector< std::vector< std::vector< std::vector< std::vector<T> > > > > TrimmingCurves;

public:

    gsTHBSplineBasis()
    { 
        representBasis(); 
    }

    gsTHBSplineBasis(gsTensorBSplineBasis<d,T> const&  tbasis, 
                     const std::vector<index_t> & boxes) 
    : gsHTensorBasis<d,T>(tbasis, boxes)
    { representBasis(); }

    gsTHBSplineBasis(gsTensorBSplineBasis<d,T> const&  tbasis, 
                     gsMatrix<T> const & boxes)
    : gsHTensorBasis<d,T>(tbasis, boxes) 
    {  representBasis(); }

    gsTHBSplineBasis( gsTensorBSplineBasis<d,T> const&  tbasis, 
                      gsMatrix<T> const & boxes, 
                      const std::vector<index_t> & levels)
    : gsHTensorBasis<d,T>(tbasis, boxes, levels)
    {  representBasis(); }

    gsTHBSplineBasis(gsBasis<T> const&  tbasis, bool manualLevels=false)
        : gsHTensorBasis<d,T>(tbasis, manualLevels)
    {  representBasis(); }

#ifdef __DOXYGEN__
    typename BoundaryBasisType::uPtr boundaryBasis(boxSide const & s);
#endif
    GISMO_UPTR_FUNCTION_DEF(BoundaryBasisType, boundaryBasis, boxSide const &)
    {
        return basisSlice(n1.direction(),n1.parameter());
    }

public:

    // Look at gsBasis.h for the documentation of this function
    gsMatrix<index_t> boundaryOffset(boxSide const & s, index_t offset ) const;

    BoundaryBasisType * basisSlice(index_t dir_fixed,T par ) const;

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

    index_t numActiveMax(const gsMatrix<T> & u, gsMatrix<index_t> & offset) const;

    // returns all actives at \a u from level \a lvl only
    // Does not clear result, only appends data
    void activeAtLevel_into(index_t lvl, const gsMatrix<T>& u,
                            std::vector<index_t> & result) const;

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

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

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

    void evalAllDers_into(const gsMatrix<T> & u, int n,
                          std::vector<gsMatrix<T> > & result,
                          bool sameElement = false) const;

    // look at eval_into
    void fastEval_into(const gsMatrix<T>& u,
                       gsMatrix<T>& result) const
    {
        gsMatrix<index_t> indices;
        this->active_into(u, indices);

        result.setZero(indices.rows(), u.cols());

        const unsigned maxLvl = this->m_tree.getMaxInsLevel() + 1;
        std::vector< gsMatrix<T> > tmpResults(maxLvl, gsMatrix<T>());
        std::vector< gsMatrix<index_t> > tmpActive(maxLvl, gsMatrix<index_t>());
        gsVector<int> processed(maxLvl);

        for (int pt = 0; pt != u.cols(); pt++)
        {
            processed.setZero();

            for (int ind = 0; ind != indices.rows(); ind++)
            {
                unsigned index = indices(ind, pt);
                if (ind != 0 && index == 0)
                    break;

                unsigned lvl = getPresLevelOfBasisFun(index);

                if (processed(lvl) == 0)
                {
                    this->m_bases[lvl]->eval_into(u.col(pt), tmpResults[lvl]);
                    this->m_bases[lvl]->active_into(u.col(pt), tmpActive[lvl]);
                    processed(lvl) = 1;
                }

                if (m_is_truncated[index] == -1)
                {
                    index_t flatTenIndx = this->flatTensorIndexOf(index, lvl);
                    int localIndex = -1;
                    for (int row = 0; row != tmpActive[lvl].rows(); ++row)
                    {
                        if (tmpActive[lvl](row, 0) == flatTenIndx)
                        {
                            localIndex = row;
                            break;
                        }
                    }
                    result(ind, pt) = tmpResults[lvl](localIndex, 0);
                }
                else // basis function is truncated
                {
                    const gsMatrix<T>& basis = tmpResults[lvl];
                    const gsMatrix<index_t>& active = tmpActive[lvl];


                    const gsSparseVector<T>& coefs = getCoefs(index);
                    T tmp = coefs(active(0, 0)) * basis(0, 0);
                    for (int i = 1; i < active.rows(); i++)
                    {
                        tmp += coefs(active(i, 0)) * basis(i, 0);
                    }
                    result(ind, pt) = tmp;
                }
            }
        }
    }

    // look at deriv_into
    void fastDeriv_into(const gsMatrix<T>& u,
                        gsMatrix<T>& result) const
    {
        gsMatrix<index_t> indices;
        this->active_into(u, indices);

        result.setZero(indices.rows() * d, u.cols());

        const unsigned maxLvl = this->m_tree.getMaxInsLevel() + 1;
        std::vector< gsMatrix<T> > tmpDeriv( maxLvl, gsMatrix<T>());
        std::vector< gsMatrix<index_t> > tmpActive(maxLvl, gsMatrix<index_t>());
        gsVector<int> processed(maxLvl);

        for (index_t pt = 0; pt != u.cols(); pt++)
        {
            processed.setZero();

            for (int ind = 0; ind != indices.rows(); ind++)
            {
                index_t index = indices(ind, pt);
                if (ind != 0 && index == 0)
                    break;

                unsigned lvl = getPresLevelOfBasisFun(index);

                if (processed(lvl) == 0)
                {
                    this->m_bases[lvl]->deriv_into(u.col(pt), tmpDeriv[lvl]);
                    this->m_bases[lvl]->active_into(u.col(pt), tmpActive[lvl]);
                    processed(lvl) = 1;
                }

                if (m_is_truncated[index] == -1)
                {
                    index_t flatTenIndx = this->flatTensorIndexOf(index, lvl);
                    int localIndex = -1;
                    for (int row = 0; row != tmpActive[lvl].rows(); ++row)
                    {
                        if (tmpActive[lvl](row, 0) == flatTenIndx)
                        {
                            localIndex = row;
                            break;
                        }
                    }

                    result.block(ind * d, pt, d, 1) =
                            tmpDeriv[lvl].block(localIndex * d, 0, d, 1);
                }
                else // basis function is truncated
                {
                    const gsMatrix<T>& basis = tmpDeriv[lvl];
                    const gsMatrix<index_t>& active = tmpActive[lvl];
                    const gsSparseVector<T>& coefs = getCoefs(index);

                    for (unsigned dim = 0; dim != d; dim++) // for all deric
                    {
                        for (index_t i = 0; i != active.rows(); ++i)
                        {
                            result(ind * d + dim, pt) +=
                                    coefs(active(i, 0)) * basis(i * d + dim, 0);
                        }
                    }
                }
            }
        }
    }

    // look at deriv2_into
    void fastDeriv2_into(const gsMatrix<T>& u,
                         gsMatrix<T>& result) const
    {
        gsMatrix<index_t> indices;
        this->active_into(u, indices);
        const unsigned numDers = (d * (d + 1)) / 2;

        result.setZero(indices.rows() * numDers, u.cols());

        const unsigned maxLvl = this->m_tree.getMaxInsLevel() + 1;
        std::vector< gsMatrix<T> > tmpDeriv2(maxLvl, gsMatrix<T>());
        std::vector< gsMatrix<index_t> > tmpActive(maxLvl, gsMatrix<index_t>());
        gsVector<int> processed(maxLvl);

        for (int pt = 0; pt != u.cols(); pt++)
        {
            processed.setZero();

            for (int ind = 0; ind != indices.rows(); ind++)
            {

                unsigned index = indices(ind, pt);
                if (ind != 0 && index == 0)
                    break;

                unsigned lvl = getPresLevelOfBasisFun(index);

                if (processed(lvl) == 0)
                {
                    this->m_bases[lvl]->deriv2_into(u.col(pt), tmpDeriv2[lvl]);
                    this->m_bases[lvl]->active_into(u.col(pt), tmpActive[lvl]);
                    processed(lvl) = 1;
                }

                if (m_is_truncated[index] == -1)
                {
                    index_t flatTenIndx = this->flatTensorIndexOf(index, lvl);
                    int localIndx = -1;
                    for (int row = 0; row != tmpActive[lvl].rows(); ++row)
                    {
                        if (tmpActive[lvl](row, 0) == flatTenIndx)
                        {
                            localIndx = row;
                            break;
                        }
                    }

                    result.block(ind * numDers, pt, numDers, 1) =
                            tmpDeriv2[lvl].block(localIndx * numDers, 0, numDers, 1);
                }
                else // basis function is truncated
                {
                    const gsMatrix<T>& basis = tmpDeriv2[lvl];
                    const gsMatrix<index_t>& active = tmpActive[lvl];
                    const gsSparseVector<T>& coefs = getCoefs(index);

                    for (unsigned der = 0; der != numDers; der++)
                    {
                        for (index_t i = 0; i != active.rows(); ++i)
                        {
                            result(ind * numDers + der, pt) +=
                                    coefs(active(i, 0)) *
                                    basis(i * numDers + der, 0);
                        }
                    }
                }
            }
        }
    }

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

    // Because of overriding one of the "eval_into" functions, all
    // functions in the base class with this name are hidden from the
    // derived class: Compiler does not search the base class as soon
    // as the function name is found in the derived class.
    // Therefore we need to add a "using" declaration to "reveal" the
    // function again. THis brings into scope all definitions of
    // eval_into from the base class.
    using gsBasis<T>::eval_into;

    unsigned numTruncated() const
    { return m_presentation.size(); }

    bool isTruncated(unsigned i) const
    {
        return (this->m_is_truncated[i] != -1);
    }

    typename std::map<index_t, gsSparseVector<T> >::const_iterator truncatedBegin() const
    { return m_presentation.begin(); }

    typename std::map<index_t, gsSparseVector<T> >::const_iterator truncatedEnd() const
    { return m_presentation.end(); }

    const gsSparseVector<T>& getCoefs(unsigned i) const
    {
        if (this->m_is_truncated[i] == -1)
        {
            GISMO_ERROR("This basis function has no sparse representation. "
                        "It is not truncated.");
        }
        else
        {
            return this->m_presentation.find(i)->second;
        }
    }


    void evalSingle_into(index_t i,
                         const gsMatrix<T>& u,
                         gsMatrix<T>& result) const;

private:

    index_t getPresLevelOfBasisFun(const index_t index) const
    {
        return (m_is_truncated[index] == -1 ?
                this->levelOf(index) : m_is_truncated[index] );
    }

    void representBasis(); // rename: precompute coeffs


    void _representBasisFunction(const unsigned j,
                                const unsigned pres_level,
                                const gsVector<index_t, d>& finest_low,
                                const gsVector<index_t, d>& finest_high);



    void _saveNewBasisFunPresentation(const gsMatrix<T>& coefs,
                                      const gsVector<index_t, d>& act_size_of_coefs,
                                      const unsigned j,
                                      const unsigned pres_level,
                                      const gsVector<index_t, d>& finest_low);



    unsigned _basisFunIndexOnLevel(const gsVector<index_t, d>& index,
                                   const unsigned level,
                                   const gsVector<index_t, d>& fin_low,
                                   const unsigned new_level);



    void _truncate(gsMatrix<T>& coefs,
                   const gsVector<index_t, d>& act_size_of_coefs,
                   const gsVector<index_t, d>& size_of_coefs,
                   const unsigned level,
                   const gsVector<index_t, d>& bspl_vec_ti,
                   const unsigned bspl_vec_ti_level,
                   const gsVector<index_t, d>& finest_low);


    unsigned _updateSizeOfCoefs(const unsigned clevel,
                                const unsigned flevel,
                                const gsVector<index_t, d>& finest_low,
                                const gsVector<index_t, d>& finest_high,
                                gsVector<index_t, d>& size_of_coefs);

public:

  short_t domainDim() const { return d; }

    GISMO_CLONE_FUNCTION(gsTHBSplineBasis)

  
  std::ostream &print(std::ostream &os) const
  {
      os << "Truncated ";
      gsHTensorBasis<d,T>::printBasic(os);
      //this->printCharMatrix(os);
      return os;
  }

    GISMO_MAKE_GEOMETRY_NEW


    void getBsplinePatchGlobal(gsVector<index_t> b1,
                          gsVector<index_t> b2,
                          unsigned level,
                          const gsMatrix<T>& geom_coef,
                          gsMatrix<T>& cp, gsKnotVector<T>& k1,
                               gsKnotVector<T>& k2) const;

  void getBsplinePatches(const gsMatrix<T>& geom_coef, gsMatrix<T>& cp, gsMatrix<index_t>& b1, gsMatrix<index_t>& b2, gsVector<index_t>& level, gsMatrix<index_t>& nvertices) const;

  gsMultiPatch<T> getBsplinePatchesToMultiPatch(const gsMatrix<T>& geom_coef) const;

    void getBsplinePatches_trimming(const gsMatrix<T>& geom_coef, gsMatrix<T>& cp, gsMatrix<index_t>& b1, gsMatrix<index_t>& b2, gsVector<index_t>& level, gsMatrix<index_t>& nvertices,
                           std::vector<std::vector<std::vector< std::vector<T> > > >& trim_curves) const;

    gsMultiPatch<T> getBsplinePatchesToMultiPatch_trimming(
            const gsMatrix<T>& geom_coef,
            std::vector<std::vector<std::vector< std::vector<T> > > >& trim_curves) const;

    /*
       * @brief Return the connected components of domain levels in knot vector indices (the boundary of a CC and the holes in the corresponding component)
       * @param[out] level levels of the boxes (level[i]: level of the i-th box,)
       * @param[out] connectedComponents the connected components in format vector<connected_component<polylines<segments<index_t> > > > where the first polyline is the outer curve and the rest are holes
    */
    void getConnectedComponents(std::vector<std::vector<std::vector< std::vector<index_t> > > >& connectedComponents, gsVector<index_t>& level) const;


   void transferbyLvl (std::vector<gsSparseMatrix<T> >& result);

    void decomposeDomain(typename gsTHBSplineBasis::AxisAlignedBoundingBox& boundaryAABB,
                         typename gsTHBSplineBasis::TrimmingCurves& trimCurves) const;

    gsTensorBSpline<d,T> getBSplinePatch(const std::vector<index_t>& boundingBox,
                                         const unsigned level,
                                         const gsMatrix<T>& geomCoefs) const;

private:
    void update_structure() 
    {
        gsHTensorBasis<d,T>::update_structure(); 
        representBasis();
    }

    // todo: rename as: representAtLevel
    void globalRefinement(const gsMatrix<T> & thbCoefs, int level, 
                          gsMatrix<T> & lvlCoefs) const;

    gsSparseMatrix<T> coarsening(const std::vector<gsSortedVector<index_t> >& old,
                           const std::vector<gsSortedVector<index_t> >& n,
                           const gsSparseMatrix<T,RowMajor> & transfer) const;

    gsSparseMatrix<T> coarsening_direct( const std::vector<gsSortedVector<index_t> >& old,
                                   const std::vector<gsSortedVector<index_t> >& n, 
                                   const std::vector<gsSparseMatrix<T,RowMajor> >& transfer) const;

    gsSparseMatrix<T> coarsening_direct2( const std::vector<gsSortedVector<index_t> >& old,
                                   const std::vector<gsSortedVector<index_t> >& n,
                                   const std::vector<gsSparseMatrix<T,RowMajor> >& transfer) const;
    

    // ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    // Utility functions for decomposeDomain
    // ......................................................................
    
    bool isFirstBoxCompletelyInsideSecond(const std::vector<index_t>& firstBox,
                                          const std::vector<index_t>& secondBox) const
    {
        return secondBox[0] < firstBox[0] && secondBox[1] < firstBox[1] &&
               firstBox[2] < secondBox[2] && firstBox[3] < secondBox[3];   
    }
    

    bool areBoxesTheSame(const std::vector<index_t>& firstBox,
                         const std::vector<index_t>& secondBox) const
    {
        return firstBox[0] == secondBox[0] && firstBox[1] == secondBox[1] &&
               firstBox[2] == secondBox[2] && firstBox[3] == secondBox[3];
    }
    
    void breakCycles(typename gsTHBSplineBasis::AxisAlignedBoundingBox& aabb,
                     typename gsTHBSplineBasis::Polylines& polylines) const;


    index_t identifyCycle(const std::vector< std::vector< T> >& polyline,
                          std::pair<T, T>& pt) const;

    
    void breakPolylineIntoTwoParts(const std::vector< std::vector<T> >& polyline,
                                   const index_t segment, 
                                   const std::pair<T, T>& pt,
                                   std::vector< std::vector<T> >& part1,
                                   std::vector< std::vector<T> >& part2) const;


    void findNewAABB(const std::vector< std::vector<T> >& polyline,
                     std::vector<index_t>& aabb) const;


private:

    template<short_t dd>
    typename util::enable_if<dd==2,void>::type
    getBsplinePatchGlobal_impl(gsVector<index_t> b1,
                               gsVector<index_t> b2,
                               unsigned level,
                               const gsMatrix<T>& geom_coef,
                               gsMatrix<T>& cp, gsKnotVector<T>& k1,
                               gsKnotVector<T>& k2) const;

    template<short_t dd>
    typename util::enable_if<dd!=2,void>::type
    getBsplinePatchGlobal_impl(gsVector<index_t> ,
                               gsVector<index_t> ,
                               unsigned ,
                               const gsMatrix<T>&,
                               gsMatrix<T>&, gsKnotVector<T>&,
                               gsKnotVector<T>&) const { GISMO_NO_IMPLEMENTATION }

    template<short_t dd>
    typename util::enable_if<dd==2,gsTensorBSpline<d,T> >::type
    getBSplinePatch_impl(const std::vector<index_t>& boundingBox,
                         const unsigned level,
                         const gsMatrix<T>& geomCoefs) const;

    template<short_t dd>
    typename util::enable_if<dd!=2,gsTensorBSpline<d,T> >::type
    getBSplinePatch_impl(const std::vector<index_t>&,
                    const unsigned,
                    const gsMatrix<T>&) const { GISMO_NO_IMPLEMENTATION }

    
private:

    // m_is_truncated(j)
    // if -1   : j-th basis function is not truncated,
    // if  \em k>0 : j-th basis function is truncated and its representation level is \em k
    gsVector<int> m_is_truncated;


    // m_presentation[j] is presentation of the j-th basis function in terms of
    // B-Splines at level m_is_truncated[j]
    //
    // if m_is_truncated[j] is equal to -1, then there is no entry
    // m_presentation[j]
    std::map<index_t, gsSparseVector<T> > m_presentation;

    using gsHTensorBasis<d,T>::m_bases;
    using gsHTensorBasis<d,T>::m_xmatrix;
    using gsHTensorBasis<d,T>::m_xmatrix_offset;
    using gsHTensorBasis<d,T>::m_deg;
};
#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsTHBSplineBasis2(pybind11::module &m);
  void pybind11_init_gsTHBSplineBasis3(pybind11::module &m);
  void pybind11_init_gsTHBSplineBasis4(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo

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