Skip to content

Commit

Permalink
Fix NURBS derivatives.
Browse files Browse the repository at this point in the history 
Signed-off-by: Seyed Ali Ghasemi <info@gha3mi.com>
  • Loading branch information
@gha3mi
gha3mi committed Jun 7, 2024
1 parent 6434ba8 commit 6df1ce7
Show file tree
Hide file tree
Showing 4 changed files with 98 additions and 74 deletions.
19 changes: 12 additions & 7 deletions src/forcad_nurbs_curve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
34 changes: 23 additions & 11 deletions src/forcad_nurbs_surface.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
40 changes: 29 additions & 11 deletions src/forcad_nurbs_volume.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
79 changes: 34 additions & 45 deletions src/forcad_utils.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
!===============================================================================


Expand Down

0 comments on commit 6df1ce7

Please sign in to comment.