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,
                                  std::vector<gsMatrix<T> >& result) const;
 
    GISMO_DEPRECATED
    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, short_t dir=-1);
 
    virtual void uniformRefine_withCoefs(gsMatrix<T>& coefs, int numKnots = 1, int mul = 1, short_t const 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 gsMatrix<T> & u) const;
    static 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