gsBoehm.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h> // just for index_t

namespace gismo
{


// =============================================================================

// Inserting knots in B-splines.

// =============================================================================


template<class T, class KnotVectorType, class Mat> 
void gsBoehm( 
    KnotVectorType & knots,
    Mat & coefs,
    T val,
    int r = 1,
    bool update_knots = true
    );
    

template<class T, class KnotVectorType, class Mat>
void gsBoehmSingle( 
    KnotVectorType & knots,
    Mat & coefs,
    T val,
    bool update_knots = true
    );



template<class T, class iter, class Mat>
void gsBoehmSingle( iter knot,   // Knot iterator
                    Mat & coefs, // Coefficients (p+1)
                    int p,       // degree
                    T val        // knot value to insert
    );


template<class KnotVectorType, class Mat, class ValIt>
void gsBoehmRefine( KnotVectorType & knots,
                    Mat & coefs, //all Coefficients
                    int p,       // degree
                    ValIt valBegin, ValIt valEnd, // knot values to insert
                    bool update_knots = true);



// =============================================================================

// Inserting knots in tensor-product B-splines.

// =============================================================================



//  Algorithm is based on gsBoehm and THE NURBS BOOK.
//
//  Some tests (for 2D and 3D) are written in gsTensorBoehm_test.
//
//  Possible improvement: function builds new matrix for coefficient (11), so it
//  uses at least [2 * memory(coefs) + epsilon] memory.
template<typename T, typename KnotVectorType, typename Mat>
void gsTensorBoehm(
        KnotVectorType& knots,
        Mat& coefs,
        T val,
        int direction,
        gsVector<unsigned> str,
        int r = 1,
        bool update_knots = true);



//
// Algorithm is based on gsBoehmRefine and THE NURBS BOOK
//
// Some tests (for 2D and 3D) are written in gsTensorBoehm_test
//
// Possible improvement: function builds new matrix for coefficients (!!), so it
// uses at least [2 * memory(coefs) + epsilon] memory.
template <typename KnotVectorType, typename Mat, typename ValIt>
void gsTensorBoehmRefine(
        KnotVectorType& knots,
        Mat& coefs,
        int direction,
        gsVector<unsigned> str,
        ValIt valBegin,
        ValIt valEnd,
        bool update_knots = true);


template <short_t d, typename KnotVectorType, typename Mat, typename ValIt>
void gsTensorBoehmRefineLocal(
        KnotVectorType& knots,
        const unsigned index,
        Mat& coefs,
        gsVector<index_t, d>& nmb_of_coefs,
        const gsVector<index_t, d>& act_size_of_coefs,
        const gsVector<index_t, d>& size_of_coefs,
        const unsigned direction,
        ValIt valBegin,
        ValIt valEnd,
        const bool update_knots);


template <short_t d, typename T, typename KnotVectorType, typename Mat>
void gsTensorInsertKnotDegreeTimes(
        const KnotVectorType& knots,
        Mat& coefs,
        const gsVector<index_t, d>& size_of_coefs,
        T val,
        const unsigned direction,
        gsVector<index_t, d>& start,
        gsVector<index_t, d>& end);



// -----------------------------------------------------------------------------
//
// Helper functions
//
// -----------------------------------------------------------------------------



namespace bspline
{

inline
int getIndex(const gsVector<unsigned>& stride,
                      const gsVector<int>& position)
{

    GISMO_ASSERT(stride.size() == position.size(),
                 "stide size is not equal to dimension of coefficients size");

    int ind = 0;
    unsigned n = stride.size();

    for (unsigned i = 0; i < n; ++i)
        ind += static_cast<int>(stride[i]) * position[i];

    return ind;
}


template<short_t d>
index_t getIndex(const gsVector<index_t, d>& stride,
                      const gsVector<index_t, d>& position)
{

    index_t ind = 0;

    for (short_t i = 0; i < d; ++i)
        ind += stride[i] * position[i];

    return ind;
}


template <typename T, typename KnotVectorType>
void computeAlpha(std::vector< std::vector<T> >& alpha,
                    const KnotVectorType& knots,
                    T value,
                    int r,
                    int k,
                    int p,
                    int s)
{

    for (int j = 1; j <= r; ++j)
    {
        int L = k - p + j;
        for (int i = 0; i <= p - j - s; ++i)
        {
            T knot_i = knots[L + i];
            alpha[j - 1][i] = (value - knot_i) / (knots[i + k + 1] - knot_i);
        }
    }

}


inline
void getLastIndex(const gsVector<unsigned>& stride,
                      const unsigned number_of_points,
                      gsVector<int>& last_point
                      )
{
    index_t stride_length = stride.size();
    for (index_t i = 0; i < stride_length; ++i)
    {
        if (i != stride_length - 1)
            // we need to subtrack 1, because we need actual last index
            last_point[i] = static_cast<int>(stride[i + 1] / stride[i]) - 1;
        else
            last_point[i] = static_cast<int>(number_of_points / stride[i]) - 1;
    }

}


template <typename T, typename KnotVectorType, typename ValIt,
          typename newKnotsType>
void computeTensorAlpha(std::vector< std::vector<T> >& alpha,
                        newKnotsType& nknots,
                        const KnotVectorType& knots,
                        ValIt valBegin,
                        ValIt valEnd,
                        bool sparse = false)
{
    const int a =  knots.iFind(*valBegin)     - knots.begin();
    const int b = (knots.iFind(*(valEnd - 1)) - knots.begin()) + 1;
    const int p = knots.degree(); // degree
    const int nik = std::distance(valBegin, valEnd); // number of inserted knots
    const int nk = knots.size();

    int i = b + p - 1;
    int k = b + p + nik - 1; // nik - 1 == r

    if (!sparse)
    {
        for (int j = 0; j <= a; j++)
            nknots[j] = knots[j];

        for (int j = b + p; j < nk; j++)
            nknots[j + nik] = knots[j];
    }


    for (int j = nik - 1; 0 <= j; --j)
    {
        const T newKnot = *(--valEnd); // X[j]

        // possible improvement: binary search instead of sequential
        while ((newKnot <= knots[i]) && (a < i))
        {
            nknots[k] = knots[i];
            k--;
            i--;
        }

        for (int ell = 1; ell <= p; ell++)
        {
            T alfa = nknots[k + ell] - newKnot;
            if (math::abs(alfa) != T(0.0))
                alfa /= nknots[k + ell] - knots[i - p + ell];

            alpha[j][ell - 1] = alfa;
        }

        nknots[k] = newKnot;
        k--;
    }
}


inline
void correctNewStride(gsVector<unsigned>& new_str,
                      const gsVector<unsigned>& str,
                      const int direction,
                      const int r)
{
    if (direction + 1 != str.size() )
    {
        new_str[direction + 1] += new_str[direction] * r;
    }
    for (index_t i = direction + 2; i < str.size(); ++i)
    {
        new_str[i] = new_str[i - 1] * (str[i] / str[i - 1]);
    }
}


template<short_t d>
void getLastIndexLocal(const gsVector<index_t, d>& size_of_coef,
                  gsVector<index_t, d>& last_point)
{
    for (short_t i = 0; i < d; ++i)
    {
        last_point[i] = size_of_coef[i] - 1;
    }
}


template <short_t d>
void buildCoeffsStrides(const gsVector<index_t, d>& size_of_coefs,
                         gsVector<index_t, d>& strides)
{
    strides(0) = 1;
    for (short_t dim = 1; dim < d; dim++)
        strides(dim) = size_of_coefs(dim - 1) * strides(dim - 1);
}



inline
unsigned numberOfIterations(const gsVector<index_t>& nmb_of_coefs,
                            const unsigned direction)
{
    unsigned nmb_of_iter = 1;
    for (index_t i = 0; i < nmb_of_coefs.size(); ++i)
    {
        if (static_cast<unsigned>(i) != direction)
            nmb_of_iter *= nmb_of_coefs[i];
    }

    return nmb_of_iter;
}

} //end namespace bspline

} //end namespace gismo

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