gsLagrangeBasis.h

#pragma once
#include <cmath>

//#include <gismo.h>
#include <gsCore/gsForwardDeclarations.h>
//#include "gsLagrangeBasis.hpp"
#include <gsUtils/gsCombinatorics.h>
#include <gsCore/gsBoundary.h>
#include <gsCore/gsBasis.h>

namespace gismo
{

// forward declaration
template< class T>  class gsLagrangeBasis;
template <class T>  class gsLagrangePoly;

/*template<short_t d, class T>
  struct gsTraits<gsLagrangeBasis<T>, d>
  {
  typedef gsLagrangeBasis<T>   TensorBasisType;
  typedef gsLagrange<T>        TensorGeometryType;
  typedef gsLagrange<T>        TensorBoundaryType;
  };

template<class T>
struct gsTraits<gsLagrangeBasis<T>,1>
{
typedef gsLagrangeBasis<T>   TensorBasisType;
typedef gsLagrange<T>        TensorGeometryType;
typedef gsLagrange<T>        TensorBoundaryType;
};*/

template<class T>
class gsLagrangeBasis : public gsBasis<T>
{
public:
    typedef gsBasis<T> Base;

    typedef T Scalar_t;

    typedef gsLagrangePoly<T> GeometryType;

    typedef gsLagrangeBasis<T> BoundaryBasisType;

    static const short_t Dim = 1;

    typedef memory::shared_ptr< gsLagrangeBasis > Ptr;

    typedef memory::unique_ptr< gsLagrangeBasis > uPtr;

//    static Ptr Make ( const gsVector<T> & vec, const T & start, const T & end)
//    { return Ptr( new gsLagrangeBasis(vec,start,end) ); }

public:

    gsLagrangeBasis()  : Base() { }

    gsLagrangeBasis ( const std::vector<T> & pars, const T & start, const T & end ) :
    m_p(pars.size()-1), m_breaks(pars), m_start(start), m_end(end)
    {
        check();
    }

    gsLagrangeBasis ( const T & start, const T & end, int amount_of_inner_breaks ) :
    m_p(amount_of_inner_breaks+1), m_start(start), m_end(end)
    {
        m_breaks.push_back(start);
        for(int i =0;i<amount_of_inner_breaks;i++)
        {
            m_breaks.push_back(
                start+( (i+1) * (end-start) / (amount_of_inner_breaks+1) ) );
        }
        m_breaks.push_back(end);
        check();
    }

    ~gsLagrangeBasis() { } //destructor

public:

    // Virtual member functions required by the base class

    short_t domainDim() const { return Dim; }

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

    void anchor_into(index_t i, gsMatrix<T>& result) const
    {
        result.resize(1,1);
        result(0,0)=m_breaks[i];
    }

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

    void numActive_into(const gsMatrix<T> & u, gsVector<index_t>& result) const
    {
        result.resize(1,u.cols());
        for(short_t i = 0; i<u.cols();++i)
        {
            result(0,i)=m_p+1;
        }
    }

    gsMatrix<index_t> boundary( ) const
    {
        gsMatrix<index_t> * res = new gsMatrix<index_t>(m_breaks.size(),1);
        for(short_t i = 0; i<=m_p;++i)
        {
            (*res)(i,0)=i;
        }
        return *res;
    }

    gsMatrix<index_t> boundary(boundary::side const & s ) const { return boundary(); }

    GISMO_UPTR_FUNCTION_DEC(gsBasis<T>, boundaryBasis, boxSide const &)

    
    gsMatrix<T> support() const
    {
        gsMatrix<T> * res = new gsMatrix<T>(1,2);
        *res << m_start , m_end ;
        return *res ;
    }


    gsMatrix<T> support(const index_t & i) const
    {
        return support();
    }

    gsMatrix<T> supportInterval(index_t dir) const;

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

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

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

    void derivSingle_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const
    {
        evalDerSingle_into(i,u,1,result);
    }

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

    void deriv2Single_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const
    {
        evalDerSingle_into(i,u,2,result);
    }

    void evalAllDers_into(const gsMatrix<T> & u, int n, gsMatrix<T>& result) const;

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

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

    GISMO_CLONE_FUNCTION(gsLagrangeBasis)

    GISMO_MAKE_GEOMETRY_NEW
    
    std::ostream &print(std::ostream &os) const
    {
        os << "Lagrange Basis: deg=" << this->degree()
           << ", size="<< this->size()
           << ", parameters= [";
        for ( unsigned i = 0; i<m_breaks.size(); ++i )
            os << m_breaks[i] << " ";
        os << "], interval= [" << m_start << "," << m_end << "].\n";
        return os;
    }


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

    // Virtual member that may be implemented or not by the derived class

    short_t size() const { return m_p+1; }

    gsLagrangeBasis<T>& component(index_t i) const
    {
        if ( i == 0 )
            return const_cast<gsLagrangeBasis&>(*this);
        else
            throw std::runtime_error("gsLagrangeBasis has only one component");
    }

    gsDomain<T> * domain() const
    { GISMO_NO_IMPLEMENTATION }

    short_t maxDegree() const
    {
        return degree();
    }

    short_t minDegree() const
    {
        return degree();
    }

    short_t degree() const
    {
        return size()-1;
    }

    gsGeometry<T> * interpolateParameters(gsMatrix<T> const& pts,
                                          gsMatrix<T> const& vals) const
    { GISMO_NO_IMPLEMENTATION }

    //MEMBER FUNCTIONS (SHOULD BE PRIVATE)

    void getTransformationLagrangeBezier(gsMatrix<T> & result) const
    {
        gsMatrix<T> L_to_M;
        _getTransformationLagrangeMonomial(L_to_M);
        gsMatrix<T> M_to_B;
        _getTransformationMonomialBezier(M_to_B);
        result=M_to_B*L_to_M;
    }

    bool check() const
    {
        bool consistent = true;
        // Checks if the breaks vector is strictly monotonic increasing.
        for(unsigned i=0;i<m_breaks.size()-1;i++)
            if(m_breaks[i]>=m_breaks[i+1])
            {
                GISMO_ERROR("gsLagrangeBasis Error: not strictly monoton increasing");
                consistent=false;
            }

        // Checks if the interval has zero length
        if(m_start==m_end)
        {
            GISMO_ERROR("gsLagrangeBasis Error: interval = 0");
            consistent=false;
        }
        return consistent;
    }

    void reparameterizeToZeroOne();

    // getters:

    short_t get_m_p() const
    {
        return m_p;
    }

    const std::vector<T> * get_m_breaks() const
    {
        return & m_breaks;
    }

    T get_m_start() const
    {
        return m_start;
    }

    T get_m_end() const
    {
        return m_end;
    }

private:

    void _getTransformationLagrangeMonomial(gsMatrix<T> & result) const;

    void _getTransformationMonomialBezier(gsMatrix<T> & result) const;

    bool _nextPoint(std::vector<int> & vec, int end) const;

    T _getFactor(int i) const
    {
        T prod = 1;
        for(short_t j = 0; j<=m_p;j++)
            if(i != j)
                prod*=(m_breaks[i]-m_breaks[j]);
        return prod;
    }

    // Data members
private:

    // Degree
    short_t m_p;
    // param vector
    std::vector<T> m_breaks;
    // Interval
    T m_start;
    T m_end;




}; // class gsLagrangeBasis

} // namespace gismo



//include the hpp file that contains the file definitions
#include "gsLagrangeBasis.hpp"

//#ifndef GISMO_HEADERS_ONLY
//#include GISMO_HPP_HEADER(gsLagrangeBasis.hpp)
//#endif