DersBasisFuns – B-spline Basis Function Library

private pure subroutine DersBasisFuns(me, u, n, ders, span)

Computes all nonzero basis functions and their derivatives based on Algorithm A2.3 in the NURBS Book [^1].

Implementation Details of Algorithm A2.3

Return value, real :: ders(j,k) $( 0 \leq j \leq p$ and $0 \leq k \leq n )$ , is the $k$ th derivative of the function $N_{i-p+j,p}(\xi)$ , where $i$ is the knot span of $\xi$ , by using the following equation,

$\begin{equation} \begin{gathered} N_{i,p}^{(k)} = \dfrac{p!}{(p - k)!} \sum \limits_{j = 0}^k {a_{k,j}} {N_{i + j,p - k}} \\ \\ a_{0,0} = 1 \\ \\ a_{k,j} = \begin{cases} \dfrac{a_{k-1,0}}{{u_{i+p-k+1}} - {u_i}} & j = 0 \\ \\ \dfrac{a_{k-1,j}-a_{k-1,j-1}}{{u_{i+p-k+1}} - {u_i}} & j = 1,\ldots,k-1 \\ \\ \dfrac{a_{k-1,k-1}}{u_{i+p+1} - u_{i+k}} & j = k \end{cases} \end{gathered} \tag{2.10} \end{equation}$

Remarks;

$k$ should not exceed $p$ (all higher derivatives are zero);
the denominators involving knot differences can become zero; the quotient is defined to be zero in this case

To avoids memory usage as much as possible in the computiation, two local arrays "a" and "ndu" introduced;

a(0:1,0:p), stores the two most recently computed rows, $a_{k,j}$ and $a_{k-1,j}$ ;
$a_{k,j}$ depend on the $a_{k-1,j}$ but not the $a_{s(< k-1),j}$ and utilize the knot differences seen in the lower triangle part of the ndu(0:p,0:p) matrix by using ${{\mathop{\rm right}}}$ and ${{\mathop{\rm left}}}$ functions
- The knot differences $={{\mathop{\rm right}} [r+1]}$ $+ {{\mathop{\rm left}} [j-r]}$ $=\xi_{i+r+1}$ $- \xi_{i+r+1-j}$
ndu(0:p,0:p), stores the basis functions and knot differences on its upper and lower triangle, respectively;


$N_{i+j-r,r}(\xi) \\ (r,j)=(0,0)$	$N_{i-1,1}(\xi) \\ (r,j)=(1,0)$	$N_{i-2,2}(\xi) \\ (r,j)=(2,0)$
${\xi_{i+1}-\xi_{i}} \\ (r,j)=(0,1)$	$N_{i,1}(\xi) \\ (r,j)=(1,1)$	$N_{i-1,2}(\xi) \\ (r,j)=(2,1)$
${\xi_{i+1}-\xi_{i-1}} \\ (r,j)=(0,2)$	${\xi_{i+2}-\xi_{i}} \\ (r,j)=(1,2)$	$N_{i,2}(\xi) \\ (r,j)=(2,2)$
ndu array (See example)

Note: This algorithm applies to rational basis functions.

Ex2.4 Computing Derivatives of the Basis Functions

Let $p = 2$ , $\mathbf{\xi} = \{0,0,0,1,2,3,4,4,5,5,5\}$ , and $\xi = 2.5$ . Then $\xi \in {[\xi_{4}, \xi_{5})}$ , and the array becomes


$N_{4,0}(2.5)=1$	$N_{3,1}(2.5)=0.5$	$N_{2,2}(2.5)=0.125$
${\xi_{5}-\xi_{4}}=1$	$N_{4,1}(2.5)=0.5$	$N_{3,2}(2.5)=0.750$
${\xi_{5}-\xi_{3}}=2$	${\xi_{6}-\xi_{4}}=2$	$N_{4,2}(2.5)=0.125$

Now compute $N_{4,2}^{(1)}(2.5)$ and $N_{4,2}^{(2)}(2.5)$ ; with $i = 4$ in Eq. (2.10), we have

$\begin{equation} a_{1,0}=\frac{1}{\xi_{6}-\xi_{4}}=\frac{1}{2} \\ a_{1,1}=-\frac{1}{\xi_{7}-\xi_{5}}=-1 \\ a_{2,0}=\frac{a_{1,0}}{\xi_{5}-\xi_{4}}=\frac{1}{2}\\ a_{2,1}=\frac{a_{1,1}-a_{1,0}}{\xi_{6}-\xi_{5}}=\frac{-1-\frac{1}{2}}{4-3}=-\frac{3}{2}\\ a_{2,2}=-\frac{a_{1,1}}{\xi_{7}-\xi_{6}}=\frac{1}{4-4}=\frac{1}{0}\\ N_{4,2}^{(1)}=2{[a_{1,0} N_{4,1}(2.5)+a_{1,1} N_{5,1}(2.5)]}\\ N_{4,2}^{(2)}=2{[a_{2,0} N_{4,0}(2.5)+a_{2,1} N_{5,0}(2.5)+a_{2,2} N_{6,0}(2.5)]} \notag \end{equation}$

Now $a_{1,1}$ , $a_{2,1}$ and $a_{2,2}$ all use knot differences which are not in the array, but they are multiplied respectively by $N_{5,1}(2.5)$ , $N_{5,0}(2.5)$ , and $N_{6,0}(2.5)$ , which are also not in the array. These terms are defined to be zero, and we are left with

$\begin{equation} N_{4,2}^{(1)}=2a_{1,0} N_{4,1}(2.5)=0.5\\ N_{4,2}^{(2)}=2a_{2,0} N_{4,0}(2.5)=1 \notag \end{equation}$

To check these values, recall from Cox-de Boor formula that $N_{4,2}(\xi)=0.5((\xi)-2)^{2}$ on $\xi \in {[2,3)}$ . The computation of $N_{3,2}^{(1)}(2.5)$ , $N_{3,2}^{(2)}(2.5)$ , $N_{2,2}^{(1)}(2.5)$ and $N_{2,2}^{(2)}(2.5)$ is analogous.

Arguments

Type	Intent	Optional		Name
class(basis),	intent(in)		::	me
real(kind=wp),	intent(in)		::	u	Given $u$ value
integer,	intent(in)		::	n	Number of derivatives ( $n \leq p$ )
real(kind=wp),	intent(out)		::	ders(0:me%p,0:n)	$= \begin{bmatrix} {N_{i - p,p}^{(0)}(u)} & \cdots &{N_{i - p,p}^{(k)}(u)}\\ \vdots & \ddots & \vdots \\ {N_{i,p}^{(0)}(u)} & \cdots &{N_{i,p}^{(k)}(u)} \end{bmatrix}$
integer,	intent(in),	optional	::	span	Knot span where $u$ lies on

Calls

Help

Source Code

DersBasisFuns

Source Code

        pure subroutine DersBasisFuns(me,u,n,ders,span)
            class(basis), intent(in)            :: me
            real(wp),   intent(in)              :: u            !< Given \( u \) value
            integer,    intent(in)              :: n            !< Number of derivatives (\(n \leq p\))
            integer,    intent(in), optional    :: span         !< Knot span where \( u \) lies on
            real(wp),   intent(out)             :: ders(0:me%p, 0:n)    !<
                                                !< $$ = \begin{bmatrix} {N_{i - p,p}^{(0)}(u)} & \cdots &{N_{i - p,p}^{(k)}(u)}\\
                                                !< \vdots & \ddots & \vdots \\ {N_{i,p}^{(0)}(u)} & \cdots &{N_{i,p}^{(k)}(u)} \end{bmatrix} $$  
            !locals
            integer :: i    !< Knot span index
            integer :: j, k, r, s1, s2, rk, pk, j1, j2
            real(wp) :: saved, temp, d
            real(wp) :: left(me%p), right(me%p)
            real(wp) :: ndu(0:me%p,0:me%p), a(0:1,0:me%p)
            
            ders = 0.0_wp
            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
                
            ndu(0,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
                    !/* Lower triangle */
                    ndu(j,r) = right(r+1) + left(j-r) 
                    temp = ndu(r,j-1) / ndu(j,r)
                    !/* Upper triangle */
                    ndu(r,j) = saved + right(r+1) * temp
                    saved = left(j-r) * temp
                end do
                ndu(j,j) = saved
            end do
            
            !/* Load the basis functions */
            ders(:,0) = ndu(:,p)
            
            !/* This section computes the derivatives (Eq. [2.9]) */
            do r = 0, p !/* Loop over function index */
                s1 = 0; s2 = 1 !/* Alternate rows in array a */
                a(0,0) = 1.0_wp
                !/* Loop to compute kth derivative */
                do k = 1, n
                    d = 0.0_wp
                    rk = r-k; pk = p-k;
                    if (r >= k) then
                        a(s2,0) = a(s1,0) / ndu(pk+1,rk)
                        d =  a(s2,0) * ndu(rk,pk)
                    end if
                    if (rk > -1) then
                        j1 = 1
                    else
                        j1 = -rk
                    end if
                    if (r-1 <= pk) then
                        j2 = k-1
                    else
                        j2 = p-r
                    end if
                    do j = j1, j2
                        a(s2,j) = (a(s1,j) - a(s1,j-1)) / ndu(pk+1,rk+j)
                        d =  d + a(s2,j) * ndu(rk+j,pk)
                    end do
                    if (r <= pk) then
                        a(s2,k) = - a(s1,k-1) / ndu(pk+1,r)
                        d =  d + a(s2,k) * ndu(r,pk)
                    end if
                    ders(r,k) = d
                    j = s1; s1 = s2; s2 = j !/* Switch rows */
                end do
            end do
            !/* Multiply through by the correct factors */
            !/* (Eq. [2.9]) */
            r = p
            do k = 1, n
                ders(:,k) = ders(:,k) * r
                r = r * (p-k)
            end do
            end associate

        end subroutine DersBasisFuns

DersBasisFuns Subroutine

Contents

Source Code

private pure subroutine DersBasisFuns(me, u, n, ders, span)

Implementation Details of Algorithm A2.3

Ex2.4 Computing Derivatives of the Basis Functions

Arguments

Calls

Contents

Source Code

Source Code