gsBasis.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsFunctionSet.h>
#include <gsCore/gsBoundary.h>//for boxSide

#define GISMO_MAKE_GEOMETRY_NEW    \
virtual memory::unique_ptr<gsGeometry<T> > makeGeometry( gsMatrix<T>coefs ) const      \
    { return memory::unique_ptr<gsGeometry<T> >(new GeometryType(*this, give(coefs))); }

namespace gismo
{

template<class T>
class gsBasis : public gsFunctionSet<T>
{

public:

    typedef gsFunctionSet<T> Base;

    typedef memory::shared_ptr< gsBasis > Ptr;

    typedef memory::unique_ptr< gsBasis > uPtr;

    typedef T Scalar_t;

    static const bool IsRational = false;

    typedef memory::unique_ptr< gsDomainIterator<T> > domainIter;

    virtual ~gsBasis();

public:

    /*
      Member functions with non-virtual implementations
      (override the _into versions in derived classes).
    */

    const gsBasis<T> & piece(const index_t k) const
    {
        GISMO_ENSURE(0==k, "Single basis is defined on single subdomain, received: "<<k );
        return *this;
    }

    //
    gsBasisFun<T> function(index_t i) const;


    gsMatrix<T> evalSingle(index_t i, const gsMatrix<T> & u) const
    {
        gsMatrix<T> result;
        this->evalSingle_into(i, u, result);
        return result;
    }

    gsMatrix<T> derivSingle(index_t i, const gsMatrix<T> & u) const
    {
        gsMatrix<T> result;
        this->derivSingle_into(i, u, result);
        return result;
    }

    gsMatrix<T> deriv2Single(index_t i, const gsMatrix<T> & u) const
    {
        gsMatrix<T> result;
        this->deriv2Single_into(i, u, result);
        return result;
    }

    gsVector<index_t> numActive(const gsMatrix<T> & u) const
    {
        gsVector<index_t> result;
        this->numActive_into(u, result);
        return result;
    }





    gsMatrix<T> evalFunc(const gsMatrix<T> & u, const gsMatrix<T> & coefs) const
    {
        gsMatrix<T> result;
        this->evalFunc_into(u, coefs, result);
        return result;
    }

    virtual void evalFunc_into(const gsMatrix<T> & u,
                               const gsMatrix<T> & coefs,
                               gsMatrix<T>& result) const;


    gsMatrix<T> derivFunc(const gsMatrix<T> & u, const gsMatrix<T> & coefs) const
    {
        gsMatrix<T> result;
        this->derivFunc_into(u, coefs, result);
        return result;
    }

    virtual void derivFunc_into(const gsMatrix<T> & u,
                                const gsMatrix<T> & coefs,
                                gsMatrix<T>& result ) const;

    virtual void jacobianFunc_into(const gsMatrix<T> & u,
                                   const gsMatrix<T> & coefs,
                                   gsMatrix<T>& result ) const;


    //  second derivatives as follows (eg. for a surface (x,y) --> (f_1, f_2, f_3):
    //     *   ( dxxf_1 dyyf_1 dxyf_1 dxxf_2 dyyf_2 dxy2f_2 dxxf_3 dyyf_3 dxyf_3 )^T
    gsMatrix<T> deriv2Func(const gsMatrix<T> & u, const gsMatrix<T> & coefs) const
    {
        gsMatrix<T> result;
        this->deriv2Func_into(u, coefs, result);
        return result;
    }

    virtual void deriv2Func_into(const gsMatrix<T> & u, const gsMatrix<T> & coefs, gsMatrix<T>& result ) const;

    virtual void evalAllDersFunc_into(const gsMatrix<T> & u,
                                      const gsMatrix<T> & coefs,
                                      const unsigned n,
                                      std::vector< gsMatrix<T> >& result,
                                      bool sameElement = false) const;

    static void linearCombination_into(const gsMatrix<T> & coefs,
                                       const gsMatrix<index_t> & actives,
                                       const gsMatrix<T> & values,
                                       gsMatrix<T> & result, bool sameElement = false);


    inline short_t dim() const {return this->domainDim();}

    /*
      Member functions that need to be redefined in the derived class.
    */

    virtual const gsMatrix<T> & weights() const 
    {
        static gsMatrix<T> dummy;
        return dummy; 
    }

    virtual gsMatrix<T> & weights()
    {
        static gsMatrix<T> dummy;
        return dummy; 
    }

    virtual bool isRational() const { return false;}

    gsMatrix<T> anchors() const
    {
        gsMatrix<T> result;
        this->anchors_into(result);
        return result;
    }

    gsMatrix<T> anchor(index_t i) const
    {
        gsMatrix<T> result;
        this->anchor_into(i, result);
        return result;
    }

    virtual void anchors_into(gsMatrix<T>& result) const;

    virtual void anchor_into(index_t i, gsMatrix<T>& result) const;

    virtual void connectivityAtAnchors(gsMesh<T> & mesh) const;

    virtual void connectivity(const gsMatrix<T> & nodes, gsMesh<T> & mesh) const;


    virtual void active_into(const gsMatrix<T> & u, gsMatrix<index_t>& result) const;

    virtual void numActive_into(const gsMatrix<T> & u, gsVector<index_t>& result) const;

    virtual bool isActive(const index_t i, const gsVector<T> & u) const;

    virtual void activeCoefs_into(const gsVector<T> & u, const gsMatrix<T> & coefs,
                                  gsMatrix<T>& result) const;


    virtual gsMatrix<index_t> allBoundary( ) const;

    virtual gsMatrix<index_t> boundaryOffset(boxSide const & s, index_t offset) const;

    gsMatrix<index_t> boundary(boxSide const & s) const
    { return this->boundaryOffset(s,0); }

    virtual index_t functionAtCorner(boxCorner const & c) const;

#ifdef __DOXYGEN__
    gsBasis<T>::uPtr boundaryBasis(boxSide const & s);
#endif
    GISMO_UPTR_FUNCTION_DEC(gsBasis<T>, boundaryBasis, boxSide const &)

    
    virtual uPtr componentBasis(boxComponent b) const;

    virtual uPtr componentBasis_withIndices(boxComponent b, gsMatrix<index_t>& indices, bool noBoundary = true) const;

    virtual gsMatrix<T> support() const;

    virtual gsMatrix<T> support(const index_t & i) const;

    gsMatrix<T> supportInterval(index_t dir) const;


    virtual void eval_into(const gsMatrix<T> & u, gsMatrix<T>& result) const;

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

    virtual void deriv_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const;

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

    virtual void deriv2_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const;

    // hessianSingle_into
    virtual void deriv2Single_into(index_t i,
                                   const gsMatrix<T> & u,
                                   gsMatrix<T>& result ) const;

    virtual void evalAllDersSingle_into(index_t i, const gsMatrix<T> & u,
                                        int n, gsMatrix<T>& result) const;

    virtual void evalDerSingle_into(index_t i, const
                                    gsMatrix<T> & u, int n,
                                    gsMatrix<T>& result) const;

    virtual gsMatrix<T> laplacian(const gsMatrix<T> & u ) const;


    GISMO_UPTR_FUNCTION_PURE(gsBasis, clone)

    virtual memory::unique_ptr<gsGeometry<T> > makeGeometry(gsMatrix<T> coefs) const = 0;

    virtual typename gsBasis::uPtr create() const;

    virtual typename gsBasis::uPtr tensorize(const gsBasis & other) const;

    virtual const gsBasis & source () const { return *this; }

    virtual gsBasis & source () { return *this; }

    virtual memory::unique_ptr<gsBasis<T> > makeNonRational() const { return clone(); }

    virtual domainIter makeDomainIterator() const;

    virtual domainIter makeDomainIterator(const boxSide & s) const;

    virtual std::ostream &print(std::ostream &os) const = 0;

    virtual std::string detail() const
    {
        // By default just uses print(..)
        std::ostringstream os;
        print(os);
        return os.str();
    }

    /*
      Member functions that may be implemented or not in the derived class
    */

    virtual size_t numElements(boxSide const & s = 0) const; // = none

    virtual size_t elementIndex(const gsVector<T> & u ) const;

    virtual gsMatrix<T> elementInSupportOf(index_t j) const;
    
    virtual const gsBasis<T> & component(short_t i) const;

    virtual gsBasis<T> & component(short_t i);

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

    virtual std::vector<index_t> asElements(gsMatrix<T> const & boxes, int refExt = 0) const;
    virtual std::vector<index_t> asElementsUnrefine(gsMatrix<T> const & boxes, int refExt = 0) const;

    virtual void refineElements(std::vector<index_t> const & boxes);
    virtual void unrefineElements(std::vector<index_t> const & boxes);

    virtual void refineElements_withCoefs(gsMatrix<T> & coefs,std::vector<index_t> const & boxes);
    virtual void unrefineElements_withCoefs(gsMatrix<T> & coefs,std::vector<index_t> const & boxes);

    virtual void uniformRefine(int numKnots = 1, int mul=1, int dir=-1);

    virtual void uniformRefine_withCoefs(gsMatrix<T>& coefs, int numKnots = 1, int mul = 1, int dir=-1);

    virtual void uniformRefine_withTransfer(gsSparseMatrix<T,RowMajor>& transfer,
                                            int numKnots = 1, int mul = 1);

    virtual void uniformCoarsen(int numKnots = 1);

    virtual void uniformCoarsen_withCoefs(gsMatrix<T>& coefs, int numKnots = 1);

    virtual void uniformCoarsen_withTransfer(gsSparseMatrix<T,RowMajor>& transfer,
                                             int numKnots = 1);

    virtual void degreeElevate(short_t const & i = 1, short_t const dir = -1);

    virtual void degreeReduce(short_t const & i = 1, short_t const dir = -1);

    virtual void degreeIncrease(short_t const & i = 1, short_t const dir = -1);

    virtual void degreeDecrease(short_t const & i = 1, short_t const dir = -1);

    void setDegree(short_t const& i);

    void setDegreePreservingMultiplicity(short_t const& i);

    virtual void elevateContinuity(int const & i = 1);

    virtual void reduceContinuity(int const & i = 1);

    virtual gsDomain<T> * domain() const;

    virtual short_t maxDegree() const;

    virtual short_t minDegree() const;

    virtual short_t totalDegree() const;

    virtual short_t degree(short_t i) const;

    memory::unique_ptr<gsGeometry<T> > interpolateData(gsMatrix<T> const& vals,
                                    gsMatrix<T> const& pts ) const;

    virtual memory::unique_ptr<gsGeometry<T> > interpolateAtAnchors(gsMatrix<T> const& vals) const;

    //gsGeometry<T> * projectL2(gsFunction<T> const & func) const;

    gsSparseMatrix<T> collocationMatrix(gsMatrix<T> const& u) const;

    std::vector<gsSparseMatrix<T> > collocationMatrixWithDeriv(const gsBasis<T> & b, const gsMatrix<T> & u);
    
    virtual void reverse();

    virtual void matchWith(const boundaryInterface & bi, const gsBasis<T> & other,
                           gsMatrix<index_t> & bndThis, gsMatrix<index_t> & bndOther, index_t offset = 0) const;


    virtual T getMinCellLength() const;

    virtual T getMaxCellLength() const;

protected:

    // inline void getLinearCombination(
    // const gsMatrix<T>         & scalars,
    // const gsMatrix<T> * const & coefs,
    // const gsMatrix<index_t>  & indices,
    // gsMatrix<T>&                result );

}; // class gsBasis

#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsBasis(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11

} // namespace gismo


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