diff --git a/src/forcad_nurbs_curve.f90 b/src/forcad_nurbs_curve.f90
index c32154b87..9d92248da 100644
--- a/src/forcad_nurbs_curve.f90
+++ b/src/forcad_nurbs_curve.f90
@@ -1592,15 +1592,17 @@ impure function compute_dTgc_nurbs_1d(Xt, knot, degree, nc, ng, Wc) result(dTgc)
     integer, intent(in) :: ng
     real(rk), intent(in), contiguous :: Wc(:)
     real(rk), allocatable :: dTgc(:,:)
-    real(rk), allocatable :: dTgci(:)
+    real(rk), allocatable :: dBi(:), Tgci(:)
     integer :: i
 
     allocate(dTgc(ng, nc))
-    allocate(dTgci(nc))
-    !$OMP PARALLEL DO PRIVATE(dTgci)
+    allocate(dBi(nc))
+    allocate(Tgci(nc))
+    !$OMP PARALLEL DO PRIVATE(Tgci, dBi)
     do i = 1, size(Xt)
-        dTgci = basis_bspline_der(Xt(i), knot, nc, degree)
-        dTgc(i,:) = dTgci*(Wc/(dot_product(dTgci,Wc)))
+        call basis_bspline_der(Xt(i), knot, nc, degree, dBi, Tgci)
+        Tgci = Tgci*(Wc/(dot_product(Tgci,Wc)))
+        dTgc(i,:) = ( dBi*Wc - Tgci*dot_product(Tgci,Wc) ) / dot_product(dBi,Wc)
     end do
     !$OMP END PARALLEL DO
 end function
@@ -1620,12 +1622,15 @@ impure function compute_dTgc_bspline_1d(Xt, knot, degree, nc, ng) result(dTgc)
     integer, intent(in) :: nc
     integer, intent(in) :: ng
     real(rk), allocatable :: dTgc(:,:)
+    real(rk), allocatable :: dTgci(:)
     integer :: i
 
     allocate(dTgc(ng, nc))
-    !$OMP PARALLEL DO
+    allocate(dTgci(nc))
+    !$OMP PARALLEL DO PRIVATE(dTgci)
     do i = 1, size(Xt)
-        dTgc(i,:) = basis_bspline_der(Xt(i), knot, nc, degree)
+        call basis_bspline_der(Xt(i), knot, nc, degree, dTgci)
+        dTgc(i,:) = dTgci
     end do
     !$OMP END PARALLEL DO
 end function
diff --git a/src/forcad_nurbs_surface.f90 b/src/forcad_nurbs_surface.f90
index 1a7d776cb..11690c7dd 100644
--- a/src/forcad_nurbs_surface.f90
+++ b/src/forcad_nurbs_surface.f90
@@ -2198,17 +2198,27 @@ impure function compute_dTgc_nurbs_2d(Xt, knot1, knot2, degree, nc, ng, Wc) resu
     integer, intent(in) :: ng(2)
     real(rk), intent(in), contiguous :: Wc(:)
     real(rk), allocatable :: dTgc(:,:)
-    real(rk), allocatable :: dTgci(:)
+    real(rk), allocatable :: dTgci(:), dTgc1(:), dTgc2(:)
+    real(rk), allocatable :: Tgci(:), Tgc1(:), Tgc2(:)
     integer :: i
 
     allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2)))
     allocate(dTgci(nc(1)*nc(2)))
-    !$OMP PARALLEL DO PRIVATE(dTgci)
+    allocate(dTgc1(nc(1)))
+    allocate(dTgc2(nc(2)))
+    allocate(Tgci(nc(1)*nc(2)))
+    allocate(Tgc1(nc(1)))
+    allocate(Tgc2(nc(2)))
+    !$OMP PARALLEL DO PRIVATE(dTgci, dTgc1, dTgc2, Tgci, Tgc1, Tgc2)
     do i = 1, size(Xt, 1)
-        dTgci = kron(&
-            basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2)),&
-            basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1)))
-        dTgc(i,:) = dTgci*(Wc/(dot_product(dTgci,Wc)))
+        call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dTgc1, Tgc1)
+        call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dTgc2, Tgc2)
+
+        dTgci = kron(dTgc2, dTgc1)
+        Tgci  = kron( Tgc2,  Tgc1)
+        Tgci = Tgci*(Wc/(dot_product(Tgci,Wc)))
+
+        dTgc(i,:) = ( dTgci*Wc - Tgci*dot_product(Tgci,Wc) ) / dot_product(dTgci,Wc)
     end do
     !$OMP END PARALLEL DO
 end function
@@ -2227,15 +2237,17 @@ impure function compute_dTgc_bspline_2d(Xt, knot1, knot2, degree, nc, ng) result
     integer, intent(in) :: degree(2)
     integer, intent(in) :: nc(2)
     integer, intent(in) :: ng(2)
-    real(rk), allocatable :: dTgc(:,:)
+    real(rk), allocatable :: dTgc(:,:), dTgc1(:), dTgc2(:)
     integer :: i
 
     allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2)))
-    !$OMP PARALLEL DO
+    allocate(dTgc1(nc(1)))
+    allocate(dTgc2(nc(2)))
+    !$OMP PARALLEL DO PRIVATE(dTgc1, dTgc2)
     do i = 1, size(Xt, 1)
-        dTgc(i,:) = kron(&
-            basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2)),&
-            basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1)))
+        call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dTgc1)
+        call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dTgc2)
+        dTgc(i,:) = kron(dTgc2, dTgc1)
     end do
     !$OMP END PARALLEL DO
 end function
diff --git a/src/forcad_nurbs_volume.f90 b/src/forcad_nurbs_volume.f90
index 174b2bb47..919302621 100644
--- a/src/forcad_nurbs_volume.f90
+++ b/src/forcad_nurbs_volume.f90
@@ -2883,17 +2883,30 @@ impure function compute_dTgc_nurbs_3d(Xt, knot1, knot2, knot3, degree, nc, ng, W
     integer, intent(in) :: ng(3)
     real(rk), intent(in), contiguous :: Wc(:)
     real(rk), allocatable :: dTgc(:,:)
-    real(rk), allocatable :: dTgci(:)
+    real(rk), allocatable :: dTgci(:), dTgc1(:), dTgc2(:), dTgc3(:)
+    real(rk), allocatable :: Tgci(:), Tgc1(:), Tgc2(:), Tgc3(:)
     integer :: i
 
     allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3)))
     allocate(dTgci(nc(1)*nc(2)*nc(3)))
-    !$OMP PARALLEL DO PRIVATE(dTgci)
+    allocate(dTgc1(nc(1)))
+    allocate(dTgc2(nc(2)))
+    allocate(dTgc2(nc(3)))
+    allocate(Tgci(nc(1)*nc(2)*nc(3)))
+    allocate(Tgc1(nc(1)))
+    allocate(Tgc2(nc(2)))
+    allocate(Tgc3(nc(3)))
+    !$OMP PARALLEL DO PRIVATE(dTgci, dTgc1, dTgc2, dTgc3, Tgci, Tgc1, Tgc2, Tgc3)
     do i = 1, size(Xt, 1)
-        dTgci = kron(basis_bspline_der(Xt(i,3), knot3, nc(3), degree(3)), kron(&
-            basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2)),&
-            basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1))))
-        dTgc(i,:) = dTgci*(Wc/(dot_product(dTgci,Wc)))
+        call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dTgc1, Tgc1)
+        call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dTgc2, Tgc2)
+        call basis_bspline_der(Xt(i,3), knot3, nc(3), degree(3), dTgc3, Tgc3)
+
+        dTgci = kron(dTgc3, kron(dTgc2, dTgc1))
+        Tgci  = kron( Tgc3, kron( Tgc2,  Tgc1))
+        Tgci = Tgci*(Wc/(dot_product(Tgci,Wc)))
+
+        dTgc(i,:) = ( dTgci*Wc - Tgci*dot_product(Tgci,Wc) ) / dot_product(dTgci,Wc)
     end do
     !$OMP END PARALLEL DO
 end function
@@ -2912,15 +2925,20 @@ impure function compute_dTgc_bspline_3d(Xt, knot1, knot2, knot3, degree, nc, ng)
     integer, intent(in) :: degree(3)
     integer, intent(in) :: nc(3)
     integer, intent(in) :: ng(3)
-    real(rk), allocatable :: dTgc(:,:)
+    real(rk), allocatable :: dTgc(:,:), dTgc1(:), dTgc2(:), dTgc3(:)
     integer :: i
 
     allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3)))
-    !$OMP PARALLEL DO
+    allocate(dTgc1(nc(1)))
+    allocate(dTgc2(nc(2)))
+    allocate(dTgc3(nc(3)))
+    !$OMP PARALLEL DO PRIVATE(dTgc1, dTgc2, dTgc3)
     do i = 1, size(Xt, 1)
-        dTgc(i,:) = kron(basis_bspline_der(Xt(i,3), knot3, nc(3), degree(3)), kron(&
-            basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2)),&
-            basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1))))
+        call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dTgc1)
+        call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dTgc2)
+        call basis_bspline_der(Xt(i,3), knot3, nc(3), degree(3), dTgc3)
+
+        dTgc(i,:) = kron(dTgc3, kron(dTgc2, dTgc1))
     end do
     !$OMP END PARALLEL DO
 end function
diff --git a/src/forcad_utils.f90 b/src/forcad_utils.f90
index ac2bc25ff..2e8f6fd3b 100644
--- a/src/forcad_utils.f90
+++ b/src/forcad_utils.f90
@@ -91,58 +91,47 @@ pure function basis_bspline(Xt, knot, nc, degree) result(B)
     !===============================================================================
     !> author: Seyed Ali Ghasemi
     !> license: BSD 3-Clause
-    pure function basis_bspline_der(Xt, knot, nc, degree) result(dB)
+    pure subroutine basis_bspline_der(Xt, knot, nc, degree, dB, B)
         integer, intent(in)   :: degree
         real(rk), intent(in), contiguous :: 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, k, n, m, jk
+        real(rk), allocatable, intent(out) :: dB(:)
+        real(rk), allocatable, intent(out), optional :: B(:)
+        real(rk), allocatable :: N(:,:), dN_dXt(:,:)
+        real(rk)              :: temp, Xth_i, Xth_i1, Xth_ip, Xth_ip1
+        integer               :: i, p
+
+        temp = abs(Xt - knot(size(knot)))
+        allocate(N(nc, 0:degree), source=0.0_rk)
+        allocate(dN_dXt(nc, 0:degree), source=0.0_rk)
+
+        do p = 0, degree
+            do concurrent (i = 1:nc)
+                Xth_i   = knot(i)
+                Xth_i1  = knot(i+1)
+                Xth_ip  = knot(i+p)
+                Xth_ip1 = knot(i+p+1)
+
+                if ( temp /= tiny(0.0_rk) .and. Xth_i <= Xt  .and. Xt <= Xth_i1 ) then
+                    N(i,0) = 1.0_rk
+                    dN_dXt(i,0) = 0.0_rk
+                end if
+                if ( Xth_ip /= Xth_i ) then
+                    N(i,p) = (Xt - Xth_i)/(Xth_ip - Xth_i) * N(i,p-1)
+                    dN_dXt(i,p) = ( N(i,p-1) + (Xt - Xth_i)*dN_dXt(i,p-1) ) / (Xth_ip - Xth_i) 
+                end if
+                if ( Xth_ip1 /= Xth_i1 ) then
+                    N(i,p) = N(i,p) + (Xth_ip1 - Xt)/(Xth_ip1 - Xth_i1) * N(i+1,p-1)
+                    dN_dXt(i,p) = dN_dXt(i,p) - ( N(i+1,p-1) - (Xth_ip1 - Xt)*dN_dXt(i+1,p-1) ) / (Xth_ip1  - Xth_i1)
+                end if
 
-        k = degree + 1
-        n = nc - 1
-        allocate(Nt(nc+degree, degree+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+degree, degree+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) < tiny(0.0_rk) ) R = 0.0_rk
-                L = b/c
-                if ( isnan(L) .or. isinf(L) .or. abs(L) < tiny(0.0_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) < tiny(0.0_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) < tiny(0.0_rk) ) Lp = 0.0_rk
-                dNt_dXt(i,jk) = Rp + Lp
             end do
         end do
-        dB = dNt_dXt(1:nc,k)
-    end function
+
+        dB = dN_dXt(:,degree)
+        if (present(B)) B = N(:,degree)
+    end subroutine
     !===============================================================================