Merge pull request #18 from gha3mi/dev

Improved Unit Tests and Bug Fixes
gha3mi · Jul 1, 2024 · 6809304 · 6809304
2 parents f59ff9c + 5438872
commit 6809304
Show file tree

Hide file tree

Showing 6 changed files with 767 additions and 171 deletions.
diff --git 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
@@ -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)
 
@@ -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)
@@ -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)
@@ -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
 !===============================================================================
 

diff --git 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
 !===============================================================================