Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improved Unit Tests and Bug Fixes #18

Merged
merged 5 commits into from
Jul 1, 2024
Merged

Improved Unit Tests and Bug Fixes #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
f868f79
Update unit tests for nurbs curve and surface
gha3mi Jun 29, 2024
e514129
Fix bug in set4 and compute_dTgc_bspline_2d_vector
gha3mi Jul 1, 2024
8c11af3
Update unit tests for nurbs_surface
gha3mi Jul 1, 2024
a17dcf1
Improve nearest_point2()
gha3mi Jul 1, 2024
5438872
Update unit tests for surface and volume
gha3mi Jul 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 5 additions & 2 deletions src/forcad_nurbs_curve.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1840,7 +1840,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
integer, intent(in) :: maxit
real(rk), intent(out) :: nearest_Xt
real(rk), allocatable, intent(out), optional :: nearest_Xg(:)
real(rk):: xk, obj, grad, hess, dk, alphak, tau, beta, lower_bounds, upper_bounds
real(rk):: xk, xkn, obj, grad, hess, dk, alphak, tau, beta, lower_bounds, upper_bounds
real(rk), allocatable :: Xg(:), Tgc(:), dTgc(:), d2Tgc(:)
integer :: k, l
logical :: convergenz
Expand All @@ -1865,6 +1865,8 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
xk = maxval(this%knot)
end if

xkn = xk

convergenz = .false.

allocate(Xg(size(this%Xc,2)))
Expand All @@ -1885,7 +1887,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
! debug
print '(i3,1x,e20.10,1x,e20.10)', k, xk, abs(grad)

if (abs(grad) <= tol) then
if (abs(grad) <= tol .or. (k>0 .and. abs(xk-xkn) <= tol)) then
convergenz = .true.
nearest_Xt = xk
if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt)
Expand All @@ -1902,6 +1904,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
l = l + 1
end do

xkn = xk
xk = xk + alphak*dk
! Check if xk is within the knot vector range
xk = max(min(xk, upper_bounds), lower_bounds)
Expand Down
15 changes: 9 additions & 6 deletions src/forcad_nurbs_surface.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -411,7 +411,7 @@ pure subroutine set4(this, degree, nc, Xc, Wc)
integer, intent(in), contiguous :: nc(:)
real(rk), intent(in), contiguous :: Xc(:,:)
real(rk), intent(in), contiguous, optional :: Wc(:)
integer :: m(3), i
integer :: m(2), i

if (allocated(this%Xc)) deallocate(this%Xc)

Expand Down Expand Up @@ -2446,7 +2446,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
real(rk), intent(out) :: nearest_Xt(2)
real(rk), allocatable, intent(out), optional :: nearest_Xg(:)
real(rk):: obj, grad(2), hess(2,2), dk(2), alphak, tau, beta, lower_bounds(2), upper_bounds(2)
real(rk), allocatable :: Xg(:), xk(:), Tgc(:), dTgc(:,:), d2Tgc(:,:)
real(rk), allocatable :: Xg(:), xk(:), xkn(:), Tgc(:), dTgc(:,:), d2Tgc(:,:)
integer :: k, l
logical :: convergenz
type(nurbs_surface) :: copy_this
Expand Down Expand Up @@ -2476,6 +2476,8 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
xk(2) = maxval(this%knot2)
end if

xkn = xk

convergenz = .false.

allocate(Xg(size(this%Xc,2)))
Expand Down Expand Up @@ -2509,7 +2511,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
! debug
print '(i3,1x,2e20.10,1x,e20.10)', k, xk, norm2(grad)

if (norm2(grad) <= tol) then
if (norm2(grad) <= tol .or. (k>0 .and. norm2(xk-xkn) <= tol)) then
convergenz = .true.
nearest_Xt = xk
if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt)
Expand All @@ -2527,6 +2529,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
l = l + 1
end do

xkn = xk
xk = xk + alphak*dk
! Check if xk is within the knot vector range
xk = max(min(xk, upper_bounds), lower_bounds)
Expand Down Expand Up @@ -2837,8 +2840,9 @@ impure subroutine compute_dTgc_bspline_2d_vector(Xt, knot1, knot2, degree, nc, n
real(rk), allocatable :: B1(:), B2(:)
integer :: i

allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2), 2))
!$OMP PARALLEL DO PRIVATE(dB1, dB2)
allocate(dTgc(ng(1)*ng(2), nc(1)*nc(2), 2), Tgc(ng(1)*ng(2), nc(1)*nc(2)))
allocate(B1(nc(1)), B2(nc(2)), dB1(nc(1)), dB2(nc(2)))

do i = 1, size(Xt, 1)
call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dB1, B1)
call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dB2, B2)
Expand All @@ -2847,7 +2851,6 @@ impure subroutine compute_dTgc_bspline_2d_vector(Xt, knot1, knot2, degree, nc, n
dTgc(i,:,1) = kron(B2, dB1)
dTgc(i,:,2) = kron(dB2, B1)
end do
!$OMP END PARALLEL DO
end subroutine
!===============================================================================

Expand Down
13 changes: 9 additions & 4 deletions src/forcad_nurbs_volume.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2967,11 +2967,13 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
real(rk), intent(out) :: nearest_Xt(3)
real(rk), allocatable, intent(out), optional :: nearest_Xg(:)
real(rk):: obj, grad(3), hess(3,3), dk(3), alphak, tau, beta, lower_bounds(3), upper_bounds(3)
real(rk), allocatable :: Xg(:), xk(:), Tgc(:), dTgc(:,:), d2Tgc(:,:)
real(rk), allocatable :: Xg(:), xk(:), xkn(:), Tgc(:), dTgc(:,:), d2Tgc(:,:)
integer :: k, l
logical :: convergenz
type(nurbs_volume) :: copy_this

alphak = 0.0_rk
dk = 0.0_rk
k = 0

! lower and upper bounds
Expand Down Expand Up @@ -3003,6 +3005,8 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
xk(3) = maxval(this%knot3)
end if

xkn = xk

convergenz = .false.

allocate(Xg(size(this%Xc,2)))
Expand Down Expand Up @@ -3052,7 +3056,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
! debug
print '(i3,1x,3e20.10,1x,e20.10)', k, xk, norm2(grad)

if (norm2(grad) <= tol) then
if (norm2(grad) <= tol .or. (k>0 .and. norm2(xk-xkn) <= tol)) then
convergenz = .true.
nearest_Xt = xk
if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt)
Expand All @@ -3070,6 +3074,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest
l = l + 1
end do

xkn = xk
xk = xk + alphak*dk
! Check if xk is within the knot vector range
xk = max(min(xk, upper_bounds), lower_bounds)
Expand Down Expand Up @@ -3581,7 +3586,7 @@ impure subroutine compute_dTgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree

allocate(dTgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3), 3))
allocate(Tgc(ng(1)*ng(2)*ng(3), nc(1)*nc(2)*nc(3)))
!$OMP PARALLEL DO PRIVATE(dB1, dB2, dB3)

do i = 1, size(Xt, 1)
call basis_bspline_der(Xt(i,1), knot1, nc(1), degree(1), dB1, B1)
call basis_bspline_der(Xt(i,2), knot2, nc(2), degree(2), dB2, B2)
Expand All @@ -3593,7 +3598,7 @@ impure subroutine compute_dTgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree
dTgc(i,:,2) = kron(kron(B3,dB2),B1)
dTgc(i,:,3) = kron(kron(dB3,B2),B1)
end do
!$OMP END PARALLEL DO

end subroutine
!===============================================================================

Expand Down
Loading
Loading