diff --git a/src/forcad_nurbs_curve.f90 b/src/forcad_nurbs_curve.f90 index 6a99ca824..35d94dd28 100644 --- a/src/forcad_nurbs_curve.f90 +++ b/src/forcad_nurbs_curve.f90 @@ -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 @@ -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))) @@ -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) @@ -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) diff --git a/src/forcad_nurbs_surface.f90 b/src/forcad_nurbs_surface.f90 index 89da83809..7bf7b1c2c 100644 --- a/src/forcad_nurbs_surface.f90 +++ b/src/forcad_nurbs_surface.f90 @@ -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 @@ -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))) @@ -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) @@ -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) diff --git a/src/forcad_nurbs_volume.f90 b/src/forcad_nurbs_volume.f90 index c2d2303f8..4d16aaf1b 100644 --- a/src/forcad_nurbs_volume.f90 +++ b/src/forcad_nurbs_volume.f90 @@ -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 @@ -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))) @@ -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) @@ -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) @@ -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) @@ -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 !===============================================================================