Add derivative and basis method

README.md

@@ -66,7 +66,8 @@ ford ford.yml
 - [x] Add `elevate_degree()` method for curves.
 - [ ] Add `elevate_degree()` method for surfaces and volumes.
 - [ ] Add `reduce_degree()` method for curves, surfaces and volumes.
-- [ ] Add `derivative()` method.
+- [x] Add `derivative()` method for curves.
+- [ ] Add `derivative()` method for surfaces and volumes.
 - [ ] Add support for multiple patches.
 - [ ] Add extraction of piecewise Bezier objects from NURBS.
 

src/NURBS/forcad_nurbs_curve.f90

@@ -1,6 +1,6 @@
 module forcad_nurbs_curve
 
-    use forcad_utils, only: rk, basis_bspline, elemConn_C0, compute_multiplicity, compute_knot_vector
+    use forcad_utils, only: rk, basis_bspline, elemConn_C0, compute_multiplicity, compute_knot_vector, basis_bspline_der
 
     implicit none
 
@@ -43,6 +43,8 @@ module forcad_nurbs_curve
         procedure :: get_nc              !!> Get number of required control points
         procedure :: insert_knot         !!> Insert a new knot
         procedure :: elevate_degree      !!> Elevate the degree of the curve
+        procedure :: derivative          !!>
+        procedure :: basis               !!>
     end type
     !===============================================================================
 
@@ -845,4 +847,96 @@ pure function factln(n) result(f)
     end function
     !===============================================================================
 
+
+    !===============================================================================
+    !> author: Seyed Ali Ghasemi
+    !> license: BSD 3-Clause
+    pure subroutine derivative(this, res, Xt, dTgc)
+        class(nurbs_curve), intent(inout) :: this
+        integer, intent(in), optional :: res
+        real(rk), intent(in), optional :: Xt(:)
+        real(rk), allocatable, intent(out) :: dTgc(:,:)
+        real(rk), allocatable :: dTgc_Xt(:)
+        integer :: i
+
+        ! check
+        if (.not.allocated(this%Xc)) then
+            error stop 'Control points are not set.'
+        end if
+
+        ! Set parameter values
+        if (present(Xt)) then
+            if (allocated(this%Xt)) deallocate(this%Xt)
+            this%Xt = Xt
+        elseif (present(res)) then
+            if (allocated(this%Xt)) deallocate(this%Xt)
+            allocate(this%Xt(res))
+            this%Xt = [(real(i-1, rk) / real(res-1, rk), i=1, res)]
+            ! else
+            ! this%Xt = this%Xt
+        end if
+
+        allocate(dTgc(size(this%Xt, 1), this%nc))
+
+        if (allocated(this%Wc)) then
+            do i = 1, size(this%Xt, 1)
+                dTgc_Xt = basis_bspline_der(this%Xt(i), this%knot, this%nc, this%order)
+                dTgc_Xt = dTgc_Xt*(this%Wc/(dot_product(dTgc_Xt,this%Wc)))
+                dTgc(i,:) = dTgc_Xt
+            end do
+        else
+            do i = 1, size(this%Xt, 1)
+                dTgc_Xt = basis_bspline_der(this%Xt(i), this%knot, this%nc, this%order)
+                dTgc(i,:) = dTgc_Xt
+            end do
+        end if
+    end subroutine
+    !===============================================================================
+
+
+    !===============================================================================
+    !> author: Seyed Ali Ghasemi
+    !> license: BSD 3-Clause
+    pure subroutine basis(this, res, Xt, Tgc)
+        class(nurbs_curve), intent(inout) :: this
+        integer, intent(in), optional :: res
+        real(rk), intent(in), optional :: Xt(:)
+        real(rk), allocatable, intent(out) :: Tgc(:,:)
+        real(rk), allocatable :: Tgci(:)
+        integer :: i
+
+        ! check
+        if (.not.allocated(this%Xc)) then
+            error stop 'Control points are not set.'
+        end if
+
+        ! Set parameter values
+        if (present(Xt)) then
+            if (allocated(this%Xt)) deallocate(this%Xt)
+            this%Xt = Xt
+        elseif (present(res)) then
+            if (allocated(this%Xt)) deallocate(this%Xt)
+            allocate(this%Xt(res))
+            this%Xt = [(real(i-1, rk) / real(res-1, rk), i=1, res)]
+            ! else
+            ! this%Xt = this%Xt
+        end if
+
+        allocate(Tgc(size(this%Xt, 1), this%nc))
+
+        if (allocated(this%Wc)) then
+            do i = 1, size(this%Xt, 1)
+                Tgci = basis_bspline(this%Xt(i), this%knot, this%nc, this%order)
+                Tgci = Tgci*(this%Wc/(dot_product(Tgci,this%Wc)))
+                Tgc(i,:) = Tgci
+            end do
+        else
+            do i = 1, size(this%Xt, 1)
+                Tgci = basis_bspline(this%Xt(i), this%knot, this%nc, this%order)
+                Tgc(i,:) = Tgci
+            end do
+        end if
+    end subroutine
+    !===============================================================================
+
 end module forcad_nurbs_curve

src/forcad_utils.f90

@@ -3,7 +3,8 @@ module forcad_utils
     implicit none
 
     private
-    public :: rk, basis_bernstein, basis_bspline, elemConn_C0, kron, ndgrid, compute_multiplicity, compute_knot_vector
+    public :: rk, basis_bernstein, basis_bspline, elemConn_C0, kron, ndgrid, compute_multiplicity, compute_knot_vector, &
+              basis_bspline_der
 
     integer, parameter :: rk = kind(1.0d0)
 
@@ -58,6 +59,63 @@ pure function basis_bspline(Xt, knot, nc, order) result(B)
     !===============================================================================
 
 
+    !===============================================================================
+    !> author: Seyed Ali Ghasemi
+    pure function basis_bspline_der(Xt, knot, nc, order) result(dB)
+        integer, intent(in)   :: order
+        real(rk), intent(in)  :: knot(:)
+        integer, intent(in)   :: nc
+        real(rk), intent(in)  :: Xt
+        real(rk), allocatable :: dB(:)
+        real(rk), allocatable :: Nt(:,:), dNt_dXt(:,:)
+        real(rk)              :: R, L, Rp, Lp, knot_i, knot_ip, knot_jk, knot_jkm, knot_end, a, b, c, d
+        integer               :: i, p, k, n, m, jk
+
+        k = order + 1
+        n = nc - 1
+        allocate(Nt(nc+order, order+1))
+        Nt = 0.0_rk
+        do i = 1, n+k
+            knot_i   = knot(i)
+            knot_ip  = knot(i+1)
+            knot_end = knot(size(knot))
+            if ( abs(Xt - knot_end) > tiny(0.0_rk) ) then
+                if ( Xt >= knot_i .and. Xt < knot_ip ) Nt(i,1) = 1.0_rk
+            elseif ( abs(Xt - knot_end) < tiny(0.0_rk) ) then
+                if ( Xt >= knot_i .and. Xt <= knot_ip ) Nt(i,1) = 1.0_rk
+            end if
+        end do
+        allocate(dNt_dXt(nc+order, order+1))
+        dNt_dXt = 0.0_rk
+        m = 0
+        do jk = 2, k
+            m = m + 1
+            do i = 1, n + k - m
+                knot_i   = knot(i)
+                knot_ip  = knot(i+1)
+                knot_jk  = knot(i+jk)
+                knot_jkm = knot(i+jk-1)
+                a        = (knot_jkm - knot_i)
+                b        = (knot_jk - Xt)
+                c        = (knot_jk - knot_ip)
+                d        = (Xt - knot_i)
+                R = d/a
+                if ( isnan(R) .or. isinf(R) .or. abs(R) < 1.0e-14_rk ) R = 0.0_rk
+                L = b/c
+                if ( isnan(L) .or. isinf(L) .or. abs(L) < 1.0e-14_rk ) L = 0.0_rk
+                Nt(i,jk) = R*Nt(i,jk-1) + L*Nt(i+1,jk-1)
+                Rp = (Nt(i,jk-1) + d*dNt_dXt(i,jk-1)) / a
+                if ( isnan(Rp) .or. isinf(Rp) .or. abs(Rp) < 1.0e-14_rk ) Rp = 0.0_rk
+                Lp = (b*dNt_dXt(i+1,jk-1) - Nt(i+1,jk-1)) / c
+                if ( isnan(Lp) .or. isinf(Lp) .or. abs(Lp) < 1.0e-14_rk ) Lp = 0.0_rk
+                dNt_dXt(i,jk) = Rp + Lp
+            end do
+        end do
+        dB = dNt_dXt(1:nc,k)
+    end function
+    !===============================================================================
+
+
     !===============================================================================
     !> author: Seyed Ali Ghasemi
     !> license: BSD 3-Clause
@@ -286,4 +344,18 @@ pure function compute_knot_vector(Xth_dir, order, continuity) result(knot)
     end function
     !===============================================================================
 
+
+    !===============================================================================
+    !> author: Seyed Ali Ghasemi
+    !> license: BSD 3-Clause
+    elemental pure function isinf(x) result(output)
+        real(rk), intent(in) :: x
+        logical :: output
+
+        output=.false.
+        if (x >  huge(x)) output=.true.
+        if (x < -huge(x)) output=.true.
+    end function
+!===============================================================================
+
 end module forcad_utils