gsBSplineSolver.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>
#include <gsCore/gsForwardDeclarations.h>

namespace gismo {


enum RootType
{
    odd,
    even,
    odd_interval,
    even_interval,
    invalid
};




template <typename T=real_t>
struct Root
{
    // requires cleanup
    RootType type;
    T parameter;
    T begParameter;
    T endParameter;
    gsVector<T> point;
    gsVector<T> begPoint;
    gsVector<T> endPoint;

    Root(bool is_odd, T par, const gsVector<T> &mpoint)
    : type(is_odd?odd:even), begParameter(0), endParameter(0)
        {
            parameter=par;
            point=mpoint;
        }
    Root(bool is_odd, T parB, T parE, const gsVector<T> &pointB, const gsVector<T> &pointE)
    : type(is_odd?odd_interval:even_interval), parameter(0)
        {
            begParameter=parB;
            endParameter=parE;
            begPoint=pointB;
            endPoint=pointE;
        }
    Root() : type(invalid){}
};


template <typename T>
unsigned findHyperPlaneIntersections (
        const gsBSpline<T>    &curve,
        const gsVector<T>     &normal,
        T                      reference,
        T                      tolerance,
        std::vector<Root<T> > &roots
        );


  
template<class T>
class gsBSplineSolver
{
    typedef gsKnotVector<T> knotVectorType;

public:
    gsBSplineSolver() : m_n(1),m_k(1), m_d(1), eps(1e-7) { }

    ~gsBSplineSolver() { } 
  
public:
    bool firstRoot(gsBSpline<T> const & bsp, int const & coord = 0,
                   T const & tr = 0, T const & tol = 1e-7, unsigned const & N=100)
    {
        initSolver(bsp,coord,tr,tol,N);
        return nextRoot();
    }

    bool nextRoot ();

    inline T value () { return x;}

    // Return a vector with all the roots
    void allRoots (gsBSpline<T> const & bsp, std::vector<T> & result,
                   int const & coord = 0, 
                   T const & tr = 0, T const & tol = 1e-7, unsigned const & N=100)
    {
        result.clear();
        initSolver(bsp,coord,tr,tol,N);

        while ( nextRoot() )
        {
            result.push_back( value() ) ;
        }
        //gsLog<< "gsBSplineSolver: found "<< result.size() <<" roots.\n  ";
    }

private:
    void initSolver(gsBSpline<T> const & bsp , int const & coord, 
                    T const & tr, T const & tol, unsigned const &N)
    {
        m_n  = bsp.coefsSize();
        m_d  = bsp.degree();
        m_k  = 1;
        eps  = tol;
        x    = 0.0;
        maxn = N + m_n;

        const knotVectorType & kv = bsp.knots();

        m_t.resize(m_n+N+m_d+1);
        m_c.resize(m_n+N);

        for ( unsigned i = 0; i < m_n; ++i )
        {
            m_c[i]= bsp.coef(i, coord) - tr;
            m_t[i]= kv[i];
        }

        for ( unsigned i = m_n; i < m_n + m_d +1 ; ++i )
            m_t[i]= kv[i];
    }

    int  insertKnot (int mu) ;

// Data members
private:

    // m_t: knot vector, m_c: coefficients
    std::vector<T> m_c, m_t;
    //gsVector<T> m_c, m_t;
    
    // m_n: coefficient size, maxn: maximum coeff. size // m_k span of root
    unsigned m_n, maxn, m_k ;

    int m_d;  // m_d: degree

    T eps; // tolerance

    T x; // root value
    
}; // class gsClass

} //namespace gismo

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