BasisFuns Subroutine

private pure subroutine BasisFuns(me, u, N, span)

Computes the nonvanishing basis functions based on Algorithm A2.2 in the NURBS Book ¹.

Implementation Details of Algorithm A2.2

Let $U=\{u_{0}, \ldots ,u_{m}\}$ be a nondecreasing sequence of real numbers, i.e., $u_{i} \leq u_{i+1}$, $i=1, \ldots ,m-1$. The $u_{i}$ are called knots, and $U$ is the knot vector. The $i$th B-spline basis function of $p$-degree (order $p+1$), denoted by $N_{i,p}(u)$, is defined by Cox-de Boor recursion formula.

$$\begin{equation} \begin{gathered} N_{i, 0}(u)=\left\{\begin{array}{ll}{1} & {\text { if } u_{i} \leq u<u_{i+1}} \\ {0} & {\text { otherwise }}\end{array}\right. \\ N_{i, p}(u)=\dfrac{u-u_{i}}{u_{i+p}-u_{i}} N_{i, p-1}(u) +\dfrac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}} N_{i+1, p-1}(u). \end{gathered} \tag{2.5} \end{equation}$$

Note that

• $N_{i,0}(u)$ is a step function, equal to zero everywhere except on the half-open interval ${[u_{i}, u_{i+1})}$;
• for $p>0$, $N_{i,p}(u)$ is a linear combination of two ($p-1$)-degree basis functions
• computation of a set of basis functions requires specification of a knot vector, $U$, and the degree, $p$;
• Equation (2.5) can yield the quotient $0/0$; we define this quotient to be zero
• the half-open interval, ${[u_{i}, u_{i+1})}$, is called the $i$th knot span; it can have zero length, since knots need not be distinct

Assuming $u$ is in the $i$th span, $i=0,\ldots,n$, computation of the nonzero functions results in an inverted triangular scheme

Algorithm computes all of the non-vanishing basis functions in an efficient way, see example, and return them in the array N(0:p). There is no risk of division by zero as we can encounter in Cox-de Boor recursion formula.

Note: This algorithm applies to rational basis functions.

Ex2.3 Nonvanishing Basis Functions

There is a great deal of redundant computation inherent in Cox-de Boor recursion formula. For example, writing out the second-degree functions in general terms, we have

$$\begin{equation} {N_{i - 2,2}}(u) = \frac{{u - {u_{i - 2}}}}{{{u_i} - {u_{i - 2}}}}{N_{i - 2,1}}(u) - \frac{{{u_{i + 1}} - u}}{{{u_{i + 1}} - {u_{i - 1}}}}{N_{i - 1,1}}(u) \tag{2.14} \end{equation}$$

$$\begin{equation} {N_{i - 1,2}}(u) = \frac{{u - {u_{i - 1}}}}{{{u_{i + 1}} - {u_{i - 1}}}}{N_{i - 1,1}}(u) - \frac{{{u_{i + 2}} - u}}{{{u_{i + 2}} - {u_i}}}{N_{i,1}}(u) \tag{2.15} \end{equation}$$

$$\begin{equation} {N_{i,2}}(u) = \frac{{u - {u_i}}}{{{u_{i + 2}} - {u_i}}}{N_{i,1}}(u) - \frac{{{u_{i + 3}} - u}}{{{u_{i + 3}} - {u_{i + 1}}}}{N_{i + 1,1}}(u) \tag{2.16} \end{equation}$$

Note that

• the first term of Eq. (2.14) and the last term of Eq. (2.16) are not computed, since $N_{i-2,1}(u)$ $= N_{i+1,1}(u)$ $= 0$;
• the expression $\frac{{N_{i - 1,1}}(u)}{{{u_{i + 1}} - {u_{i - 1}}}}$ which appears in the second term of Eq. (2.14) appears in the first term of Eq. (2.15); a similar statement holds for the second term of Eq. (2.15) and the first term of Eq. (2.16).

We introduce the notation

$$\begin{equation} \mathop{\rm left}[j] = u - {u_{i + 1 - j}}\qquad \mathop{\rm right}[j] = {u_{i + j}} - u \notag \end{equation}$$

Eq. (2.15-2.16) are then

$$\begin{equation} {N_{i - 2,2}}(u) = \frac{{{\mathop{\rm left}} [ 3 ]}}{{{\mathop{\rm right}} [ 0 ] + {\mathop{\rm left}} [ 3 ]}}{N_{i - 2,1}}(u) + \frac{{{\mathop{\rm right}} [ 1 ]}}{{{\mathop{\rm right}} [ 1 ] + {\mathop{\rm left}} [ 2 ]}}{N_{i - 1,1}}(u) \notag \end{equation}$$

$$\begin{equation} {N_{i - 1,2}}(u) = \frac{{{\mathop{\rm left}} [ 2 ]}}{{{\mathop{\rm right}} [ 1 ] + {\mathop{\rm left}} [ 2 ]}}{N_{i - 1,1}}(u) + \frac{{{\mathop{\rm right}} [ 2 ]}}{{{\mathop{\rm right}} [ 2 ] + {\mathop{\rm left}} [ 1 ]}}{N_{i,1}}(u) \notag \end{equation}$$

$$\begin{equation} {N_{i,2}}(u) = \frac{{{\mathop{\rm left}} [ 1 ]}}{{{\mathop{\rm right}} [ 2 ] + {\mathop{\rm left}} [ 1 ]}}{N_{i,1}}(u) + \frac{{{\mathop{\rm right}} [ 3 ]}}{{{\mathop{\rm right}} [ 3 ] + {\mathop{\rm left}} [ 0 ]}}{N_{i + 1,1}}(u) \notag \end{equation}$$

Arguments

        pure subroutine BasisFuns(me,u,N,span)
            class(basis), intent(in)                :: me
            real(wp),       intent(in)              :: u            !< Given \( u \) value
            integer,        intent(in), optional    :: span         !< Knot span where \( u \) lies on
            real(wp),       intent(out)             :: N(0:me%p)    !< \( =  \{ N_{i-p,p}(u) \ldots N_{i,p}(u)\} \)
            !locals
            integer     :: i    !< Knot span index
            integer     :: j, r
            real(wp)    :: left(me%p), right(me%p), saved, temp!, w
            
            associate( p => me%p, Ui => me%kv ) !get rid of % signs
                
            ! set the knot span
            if(present(span))then
                i = span
            else
                i = me%FindSpan(u)
            end if

            N = 0.0_wp; right=0.0_wp; left=0.0_wp
            
            N(0) = 1.0_wp
            do j = 1, p
            left(j)  = u - Ui(i+1-j)
            right(j) = Ui(i+j) - u
            saved = 0.0_wp
            do r = 0, j-1
                temp = N(r) / (right(r+1) + left(j-r))
                N(r) = saved + right(r+1) * temp
                saved = left(j-r) * temp
            end do
            N(j) = saved
            end do

            ! if(me%rational)then
            !     w = sum(N(:)*me%w(i-p:i))
            !     do j = 0, me%p
            !         N(j) = N(j)*me%w(j) / w
            !     end do
            ! end if
                
            end associate
            
        end subroutine BasisFuns