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 b58223e70..7bf7b1c2c 100644 --- 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 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 !=============================================================================== diff --git a/test/test_nurbs_curve.f90 b/test/test_nurbs_curve.f90 index b3fe3661b..8ee719c80 100644 --- a/test/test_nurbs_curve.f90 +++ b/test/test_nurbs_curve.f90 @@ -9,32 +9,55 @@ program test_nurbs_curve real(rk), allocatable :: Xg(:,:), Xgb(:,:) real(rk) :: knot(6) integer, allocatable :: elemConn(:,:) - real(rk), allocatable :: Tgc(:,:), dTgc(:,:) + real(rk), allocatable :: Tgc(:,:), dTgc(:,:), Tgcb(:,:), dTgcb(:,:), d2Tgc(:,:), d2Tgcb(:,:) + real(rk), allocatable :: Tgc1(:), dTgc1(:), Tgc1b(:), dTgc1b(:), d2Tgc1(:), d2Tgc1b(:) integer :: i, id real(rk), allocatable :: nearest_Xg(:) - real(rk) :: nearest_Xt + real(rk) :: nearest_Xt, length, lengthb type(unit_test) :: ut allocate(Xc(3, 3)) Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk] - Xc(2,:) = [0.0_rk, 5.0_rk, 0.0_rk] - Xc(3,:) = [5.0_rk, 5.0_rk, 0.0_rk] + Xc(2,:) = [1.0_rk, 0.0_rk, 0.0_rk] + Xc(3,:) = [2.0_rk, 0.0_rk, 0.0_rk] allocate(Wc(3)) - Wc = [1.0_rk, 2.0_rk, 0.3_rk] + Wc = [1.0_rk, 0.9_rk, 0.8_rk] knot = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk, 1.0_rk] call nurbs%set(knot, Xc, Wc) call bsp%set(knot, Xc) + call nurbs%set(degree=nurbs%get_degree(), nc=nurbs%get_nc(), Xc=nurbs%get_Xc(), Wc=nurbs%get_Wc()) + call bsp%set(degree=nurbs%get_degree(), nc=nurbs%get_nc(), Xc=nurbs%get_Xc()) + call nurbs%create(res = 23) call bsp%create(res = 23) + call nurbs%cmp_length(length) + call bsp%cmp_length(lengthb) + + call ut%check(res=length, expected=2.0_rk, tol=1e-5_rk, msg="test_nurbs_curve: 01") + call ut%check(res=lengthb, expected=2.0_rk, tol=1e-5_rk, msg="test_nurbs_curve: 02") + call nurbs%nearest_point([0.0_rk, 0.0_rk, 0.5_rk], nearest_Xg, nearest_Xt, id) - call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_curve: 01") - call ut%check(res=nearest_Xt, expected=0.0_rk, tol=1e-5_rk, msg="test_nurbs_curve: 02") - call ut%check(res=id, expected=1, msg="test_nurbs_curve: 03") + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_curve: 03") + call ut%check(res=nearest_Xt, expected=0.0_rk, tol=1e-5_rk, msg="test_nurbs_curve: 04") + call ut%check(res=id, expected=1, msg="test_nurbs_curve: 05") + + call bsp%nearest_point([0.0_rk, 0.0_rk, 0.5_rk], nearest_Xg, nearest_Xt, id) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_curve: 06") + call ut%check(res=nearest_Xt, expected=0.0_rk, tol=1e-5_rk, msg="test_nurbs_curve: 07") + call ut%check(res=id, expected=1, msg="test_nurbs_curve: 08") + + call nurbs%nearest_point2([0.0_rk, 0.0_rk, 0.5_rk], 1e-10_rk, 10, nearest_Xt, nearest_Xg) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_curve: 09") + call ut%check(res=nearest_Xt, expected=0.0_rk, tol=1e-5_rk, msg="test_nurbs_curve: 10") + + call bsp%nearest_point2([0.0_rk, 0.0_rk, 0.5_rk], 1e-10_rk, 10, nearest_Xt, nearest_Xg) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_curve: 11") + call ut%check(res=nearest_Xt, expected=0.0_rk, tol=1e-5_rk, msg="test_nurbs_curve: 12") Xg = nurbs%get_Xg() Xgb = bsp%get_Xg() @@ -42,112 +65,111 @@ program test_nurbs_curve call nurbs%set([0.0_rk, 1.0_rk], 2, [-1, -1], Xc, Wc) call bsp%set([0.0_rk, 1.0_rk], 2, [-1, -1], Xc) - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 04") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 05") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 13") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 14") call nurbs%set(Xc, Wc) call bsp%set(Xc) - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 06") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 07") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 15") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 16") call nurbs%create(Xt = nurbs%get_Xt()) call bsp%create(Xt = bsp%get_Xt()) - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 08") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 09") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 17") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 18") - call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 10") - call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 11") + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 19") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 20") - call ut%check(res=nurbs%get_Xc(1), expected=Xc(1,:), tol=1e-5_rk, msg="test_nurbs_curve: 12") - call ut%check(res=bsp%get_Xc(1), expected=Xc(1,:), tol=1e-5_rk, msg="test_nurbs_curve: 13") + call ut%check(res=nurbs%get_Xc(1), expected=Xc(1,:), tol=1e-5_rk, msg="test_nurbs_curve: 21") + call ut%check(res=bsp%get_Xc(1), expected=Xc(1,:), tol=1e-5_rk, msg="test_nurbs_curve: 22") - call ut%check(res=nurbs%get_Xc(1,1), expected=Xc(1,1), tol=1e-5_rk, msg="test_nurbs_curve: 14") - call ut%check(res=bsp%get_Xc(1,1), expected=Xc(1,1), tol=1e-5_rk, msg="test_nurbs_curve: 15") + call ut%check(res=nurbs%get_Xc(1,1), expected=Xc(1,1), tol=1e-5_rk, msg="test_nurbs_curve: 23") + call ut%check(res=bsp%get_Xc(1,1), expected=Xc(1,1), tol=1e-5_rk, msg="test_nurbs_curve: 24") - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 16") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 17") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 25") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 26") - call ut%check(res=nurbs%get_Xg(1), expected=Xg(1,:), tol=1e-5_rk, msg="test_nurbs_curve: 18") - call ut%check(res=bsp%get_Xg(1), expected=Xgb(1,:), tol=1e-5_rk, msg="test_nurbs_curve: 19") + call ut%check(res=nurbs%get_Xg(1), expected=Xg(1,:), tol=1e-5_rk, msg="test_nurbs_curve: 27") + call ut%check(res=bsp%get_Xg(1), expected=Xgb(1,:), tol=1e-5_rk, msg="test_nurbs_curve: 28") - call ut%check(res=nurbs%get_Xg(1,1), expected=Xg(1,1), tol=1e-5_rk, msg="test_nurbs_curve: 20") - call ut%check(res=bsp%get_Xg(1,1), expected=Xgb(1,1), tol=1e-5_rk, msg="test_nurbs_curve: 21") + call ut%check(res=nurbs%get_Xg(1,1), expected=Xg(1,1), tol=1e-5_rk, msg="test_nurbs_curve: 29") + call ut%check(res=bsp%get_Xg(1,1), expected=Xgb(1,1), tol=1e-5_rk, msg="test_nurbs_curve: 30") - call ut%check(res=nurbs%get_Wc(), expected=Wc, tol=1e-5_rk, msg="test_nurbs_curve: 22") + call ut%check(res=nurbs%get_Wc(), expected=Wc, tol=1e-5_rk, msg="test_nurbs_curve: 31") - call ut%check(res=nurbs%get_Wc(1), expected=Wc(1), tol=1e-5_rk, msg="test_nurbs_curve: 23") + call ut%check(res=nurbs%get_Wc(1), expected=Wc(1), tol=1e-5_rk, msg="test_nurbs_curve: 32") - call ut%check(res=nurbs%get_knot(), expected=knot, tol=1e-5_rk, msg="test_nurbs_curve: 24") - call ut%check(res=bsp%get_knot(), expected=knot, tol=1e-5_rk, msg="test_nurbs_curve: 25") + call ut%check(res=nurbs%get_knot(), expected=knot, tol=1e-5_rk, msg="test_nurbs_curve: 33") + call ut%check(res=bsp%get_knot(), expected=knot, tol=1e-5_rk, msg="test_nurbs_curve: 34") - call ut%check(res=nurbs%get_knot(1), expected=knot(1), tol=1e-5_rk, msg="test_nurbs_curve: 26") - call ut%check(res=bsp%get_knot(1), expected=knot(1), tol=1e-5_rk, msg="test_nurbs_curve: 27") + call ut%check(res=nurbs%get_knot(1), expected=knot(1), tol=1e-5_rk, msg="test_nurbs_curve: 35") + call ut%check(res=bsp%get_knot(1), expected=knot(1), tol=1e-5_rk, msg="test_nurbs_curve: 36") - call ut%check(res=nurbs%get_ng(), expected=size(Xg,1), msg="test_nurbs_curve: 28") - call ut%check(res=bsp%get_ng(), expected=size(Xgb,1), msg="test_nurbs_curve: 29") + call ut%check(res=nurbs%get_ng(), expected=size(Xg,1), msg="test_nurbs_curve: 37") + call ut%check(res=bsp%get_ng(), expected=size(Xgb,1), msg="test_nurbs_curve: 38") - call ut%check(res=nurbs%get_degree(), expected=2, msg="test_nurbs_curve: 30") - call ut%check(res=bsp%get_degree(), expected=2, msg="test_nurbs_curve: 31") + call ut%check(res=nurbs%get_degree(), expected=2, msg="test_nurbs_curve: 39") + call ut%check(res=bsp%get_degree(), expected=2, msg="test_nurbs_curve: 40") - call ut%check(res=nurbs%get_multiplicity(), expected=[3,3], msg="test_nurbs_curve: 32") - call ut%check(res=bsp%get_multiplicity(), expected=[3,3], msg="test_nurbs_curve: 33") + call ut%check(res=nurbs%get_multiplicity(), expected=[3,3], msg="test_nurbs_curve: 41") + call ut%check(res=bsp%get_multiplicity(), expected=[3,3], msg="test_nurbs_curve: 42") - call ut%check(res=nurbs%get_continuity(), expected=[-1,-1], msg="test_nurbs_curve: 34") - call ut%check(res=bsp%get_continuity(), expected=[-1,-1], msg="test_nurbs_curve: 35") + call ut%check(res=nurbs%get_continuity(), expected=[-1,-1], msg="test_nurbs_curve: 43") + call ut%check(res=bsp%get_continuity(), expected=[-1,-1], msg="test_nurbs_curve: 44") - call ut%check(res=nurbs%get_nc(), expected=size(Xc,1), msg="test_nurbs_curve: 36") - call ut%check(res=bsp%get_nc(), expected=size(Xc,1), msg="test_nurbs_curve: 37") + call ut%check(res=nurbs%get_nc(), expected=size(Xc,1), msg="test_nurbs_curve: 45") + call ut%check(res=bsp%get_nc(), expected=size(Xc,1), msg="test_nurbs_curve: 46") call nurbs%cmp_nc() call bsp%cmp_nc() - call ut%check(res=nurbs%get_nc(), expected=size(Xc,1), msg="test_nurbs_curve: 38") - call ut%check(res=bsp%get_nc(), expected=size(Xc,1), msg="test_nurbs_curve: 39") - + call ut%check(res=nurbs%get_nc(), expected=size(Xc,1), msg="test_nurbs_curve: 47") + call ut%check(res=bsp%get_nc(), expected=size(Xc,1), msg="test_nurbs_curve: 48") elemConn = nurbs%cmp_elem_Xc_vis(2) call nurbs%set_elem_Xc_vis(elemConn) - call ut%check(res=nurbs%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_curve: 40") + call ut%check(res=nurbs%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_curve: 49") deallocate(elemConn) elemConn = nurbs%cmp_elem_Xc_vis() call nurbs%set_elem_Xc_vis(elemConn) - call ut%check(res=nurbs%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_curve: 41") + call ut%check(res=nurbs%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_curve: 50") deallocate(elemConn) elemConn = bsp%cmp_elem_Xc_vis(2) call bsp%set_elem_Xc_vis(elemConn) - call ut%check(res=bsp%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_curve: 42") + call ut%check(res=bsp%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_curve: 51") deallocate(elemConn) elemConn = bsp%cmp_elem_Xc_vis() call bsp%set_elem_Xc_vis(elemConn) - call ut%check(res=bsp%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_curve: 43") + call ut%check(res=bsp%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_curve: 52") deallocate(elemConn) elemConn = nurbs%cmp_elem_Xg_vis(2) call nurbs%set_elem_Xg_vis(elemConn) - call ut%check(res=nurbs%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_curve: 44") + call ut%check(res=nurbs%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_curve: 53") deallocate(elemConn) elemConn = nurbs%cmp_elem_Xg_vis() call nurbs%set_elem_Xg_vis(elemConn) - call ut%check(res=nurbs%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_curve: 45") + call ut%check(res=nurbs%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_curve: 54") deallocate(elemConn) elemConn = bsp%cmp_elem_Xg_vis(2) call bsp%set_elem_Xg_vis(elemConn) - call ut%check(res=bsp%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_curve: 46") + call ut%check(res=bsp%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_curve: 55") deallocate(elemConn) elemConn = bsp%cmp_elem_Xg_vis() call bsp%set_elem_Xg_vis(elemConn) - call ut%check(res=bsp%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_curve: 47") + call ut%check(res=bsp%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_curve: 56") deallocate(elemConn) elemConn = nurbs%cmp_elem() call nurbs%set_elem(elemConn) - call ut%check(res=nurbs%get_elem(), expected=elemConn, msg="test_nurbs_curve: 48") + call ut%check(res=nurbs%get_elem(), expected=elemConn, msg="test_nurbs_curve: 57") deallocate(elemConn) elemConn = bsp%cmp_elem() call bsp%set_elem(elemConn) - call ut%check(res=bsp%get_elem(), expected=elemConn, msg="test_nurbs_curve: 49") + call ut%check(res=bsp%get_elem(), expected=elemConn, msg="test_nurbs_curve: 58") deallocate(elemConn) call nurbs%modify_Xc(Xc(1,1), 1,1) @@ -158,20 +180,41 @@ program test_nurbs_curve call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 50") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 51") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 59") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 60") call nurbs%basis(res=23, Tgc=Tgc) call bsp%basis(res=23, Tgc=Tgc) + call nurbs%basis(Xt=0.0_rk, Tgc=Tgc1) + call bsp%basis(Xt=0.0_rk, Tgc=Tgc1b) + call nurbs%basis(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], Tgc=Tgc) call bsp%basis(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], Tgc=Tgc) call nurbs%derivative(res=23, dTgc=dTgc, Tgc=Tgc) - call bsp%derivative(res=23, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative(res=23, dTgc=dTgcb, Tgc=Tgcb) call nurbs%derivative(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], dTgc=dTgc, Tgc=Tgc) - call bsp%derivative(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], dTgc=dTgc, Tgc=Tgc) + call bsp%derivative(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative(Xt=0.0_rk, dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative(Xt=0.0_rk, dTgc=dTgc1b, Tgc=Tgc1b) + + call nurbs%derivative(Xt=0.0_rk, dTgc=dTgc1b, Tgc=Tgc1b, elem=[1,2,3]) + call bsp%derivative(Xt=0.0_rk, dTgc=dTgc1b, Tgc=Tgc1b, elem=[1,2,3]) + + call nurbs%derivative2(res=23, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative2(res=23, d2Tgc=d2Tgcb, dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative2(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative2(Xt=[(real(i-1, rk) / real(23-1, rk), i=1, 23)], d2Tgc=d2Tgcb, dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative2(Xt=0.0_rk, d2Tgc=d2Tgc1, dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative2(Xt=0.0_rk, d2Tgc=d2Tgc1b, dTgc=dTgc1b, Tgc=Tgc1b) + + call nurbs%derivative2(Xt=0.0_rk, d2Tgc=d2Tgc1, dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative2(Xt=0.0_rk, d2Tgc=d2Tgc1b, dTgc=dTgc1b, Tgc=Tgc1b) call nurbs%rotate_Xc(45.0_rk, 0.0_rk, 0.0_rk) call nurbs%rotate_Xc(-45.0_rk, 0.0_rk, 0.0_rk) @@ -179,8 +222,8 @@ program test_nurbs_curve call bsp%rotate_Xc(45.0_rk, 0.0_rk, 0.0_rk) call bsp%rotate_Xc(-45.0_rk, 0.0_rk, 0.0_rk) - call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 52") - call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 53") + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 61") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 62") call nurbs%rotate_Xc(0.0_rk, 45.0_rk, 0.0_rk) call nurbs%rotate_Xc(0.0_rk, -45.0_rk, 0.0_rk) @@ -188,8 +231,8 @@ program test_nurbs_curve call bsp%rotate_Xc(0.0_rk, 45.0_rk, 0.0_rk) call bsp%rotate_Xc(0.0_rk, -45.0_rk, 0.0_rk) - call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 54") - call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 55") + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 63") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 64") call nurbs%rotate_Xc(0.0_rk, 0.0_rk, 45.0_rk) call nurbs%rotate_Xc(0.0_rk, 0.0_rk, -45.0_rk) @@ -197,8 +240,8 @@ program test_nurbs_curve call bsp%rotate_Xc(0.0_rk, 0.0_rk, 45.0_rk) call bsp%rotate_Xc(0.0_rk, 0.0_rk, -45.0_rk) - call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 56") - call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 57") + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 65") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 66") call nurbs%rotate_Xg(45.0_rk, 0.0_rk, 0.0_rk) call nurbs%rotate_Xg(-45.0_rk, 0.0_rk, 0.0_rk) @@ -206,8 +249,8 @@ program test_nurbs_curve call bsp%rotate_Xg(45.0_rk, 0.0_rk, 0.0_rk) call bsp%rotate_Xg(-45.0_rk, 0.0_rk, 0.0_rk) - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 58") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 59") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 67") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 68") call nurbs%rotate_Xg(0.0_rk, 45.0_rk, 0.0_rk) call nurbs%rotate_Xg(0.0_rk, -45.0_rk, 0.0_rk) @@ -215,8 +258,8 @@ program test_nurbs_curve call bsp%rotate_Xg(0.0_rk, 45.0_rk, 0.0_rk) call bsp%rotate_Xg(0.0_rk, -45.0_rk, 0.0_rk) - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 60") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 61") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 69") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 70") call nurbs%rotate_Xg(0.0_rk, 0.0_rk, 45.0_rk) call nurbs%rotate_Xg(0.0_rk, 0.0_rk, -45.0_rk) @@ -224,8 +267,8 @@ program test_nurbs_curve call bsp%rotate_Xg(0.0_rk, 0.0_rk, 45.0_rk) call bsp%rotate_Xg(0.0_rk, 0.0_rk, -45.0_rk) - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 62") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 63") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 71") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 72") call nurbs%translate_Xc([5.0_rk, 5.0_rk, 5.0_rk]) call nurbs%translate_Xc([-5.0_rk, -5.0_rk, -5.0_rk]) @@ -233,8 +276,8 @@ program test_nurbs_curve call bsp%translate_Xc([5.0_rk, 5.0_rk, 5.0_rk]) call bsp%translate_Xc([-5.0_rk, -5.0_rk, -5.0_rk]) - call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 64") - call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 65") + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 73") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_curve: 74") call nurbs%translate_Xg([5.0_rk, 5.0_rk, 5.0_rk]) call nurbs%translate_Xg([-5.0_rk, -5.0_rk, -5.0_rk]) @@ -242,8 +285,8 @@ program test_nurbs_curve call bsp%translate_Xg([5.0_rk, 5.0_rk, 5.0_rk]) call bsp%translate_Xg([-5.0_rk, -5.0_rk, -5.0_rk]) - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 66") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 67") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 75") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 76") call nurbs%export_Xc("vtk/test_nurbs_curve_Xc.vtk") call nurbs%export_Xg("vtk/test_nurbs_curve_Xg.vtk") @@ -257,8 +300,8 @@ program test_nurbs_curve call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 68") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 69") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 77") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 78") call nurbs%elevate_degree(2) call bsp%elevate_degree(2) @@ -266,8 +309,8 @@ program test_nurbs_curve call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 70") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 71") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 79") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 80") call nurbs%remove_knots([0.25_rk, 0.75_rk], [2,1]) call bsp%remove_knots([0.25_rk, 0.75_rk], [2,1]) @@ -275,8 +318,8 @@ program test_nurbs_curve call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 72") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 73") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_curve: 81") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_curve: 82") call nurbs%set_circle([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk) @@ -287,4 +330,10 @@ program test_nurbs_curve call bsp%finalize() deallocate(Xc, Wc, Xg, Xgb) + call nurbs%set(knot=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], Xc=[0.0_rk,2.0_rk], Wc=[1.0_rk, 0.9_rk]) + call bsp%set(knot=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], Xc=[0.0_rk,2.0_rk]) + + call nurbs%finalize() + call bsp%finalize() + end program diff --git a/test/test_nurbs_surface.f90 b/test/test_nurbs_surface.f90 index 19ce0398d..159738bd9 100644 --- a/test/test_nurbs_surface.f90 +++ b/test/test_nurbs_surface.f90 @@ -5,40 +5,315 @@ program test_nurbs_surface implicit none type(nurbs_surface) :: nurbs, bsp - real(rk) :: knot1(6), knot2(6) + real(rk) :: knot1(4), knot2(4), area, areab real(rk), allocatable :: Xc(:,:), Wc(:) real(rk), allocatable :: Xg(:,:), Xgb(:,:) + integer, allocatable :: elemConn(:,:) + real(rk), allocatable :: Tgc(:,:), dTgc(:,:,:), Tgcb(:,:), dTgcb(:,:,:), d2Tgc(:,:,:), d2Tgcb(:,:,:) + real(rk), allocatable :: Tgc1(:), dTgc1(:,:), Tgc1b(:), dTgc1b(:,:), d2Tgc1(:,:), d2Tgc1b(:,:) + integer :: i, id + real(rk), allocatable :: nearest_Xg(:), nearest_Xt(:) type(unit_test) :: ut - Xc = generate_Xc(3, 3, 1.0_rk) + allocate(Xc(4, 3)) + Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk] + Xc(2,:) = [5.0_rk, 0.0_rk, 0.0_rk] + Xc(3,:) = [0.0_rk, 5.0_rk, 0.0_rk] + Xc(4,:) = [5.0_rk, 5.0_rk, 0.0_rk] allocate(Wc(size(Xc, 1))) Wc = 1.0_rk - Wc(2) = 2.0_rk + Wc(2) = 0.9_rk - knot1 = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk, 1.0_rk] - knot2 = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk, 1.0_rk] + knot1 = [0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk] + knot2 = [0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk] call nurbs%set(knot1, knot2, Xc, Wc) call bsp%set(knot1, knot2, Xc) + call nurbs%set(degree=nurbs%get_degree(), nc=nurbs%get_nc(), Xc=nurbs%get_Xc(), Wc=nurbs%get_Wc()) + call bsp%set(degree=nurbs%get_degree(), nc=nurbs%get_nc(), Xc=nurbs%get_Xc()) + call nurbs%create(30, 30) call bsp%create(30, 30) + call nurbs%cmp_area(area) + call bsp%cmp_area(areab) + + call ut%check(res=area, expected=25.0_rk, tol=1e-5_rk, msg="test_nurbs_surface: 01") + call ut%check(res=areab, expected=25.0_rk, tol=1e-5_rk, msg="test_nurbs_surface: 02") + + call nurbs%nearest_point([0.0_rk, 0.0_rk, 0.5_rk], nearest_Xg, nearest_Xt, id) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_surface: 03") + call ut%check(res=nearest_Xt, expected=[0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_surface: 04") + call ut%check(res=id, expected=1, msg="test_nurbs_surface: 05") + + call bsp%nearest_point([0.0_rk, 0.0_rk, 0.5_rk], nearest_Xg, nearest_Xt, id) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_surface: 06") + call ut%check(res=nearest_Xt, expected=[0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_surface: 07") + call ut%check(res=id, expected=1, msg="test_nurbs_surface: 08") + + call nurbs%nearest_point2([0.0_rk, 0.0_rk, 0.5_rk], 1e-10_rk, 10, nearest_Xt, nearest_Xg) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_surface: 09") + call ut%check(res=nearest_Xt, expected=[0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_surface: 10") + + call bsp%nearest_point2([0.0_rk, 0.0_rk, 0.5_rk], 1e-10_rk, 10, nearest_Xt, nearest_Xg) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_surface: 11") + call ut%check(res=nearest_Xt, expected=[0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_surface: 12") + Xg = nurbs%get_Xg() Xgb = bsp%get_Xg() - call nurbs%insert_knots(1, [0.25_rk, 0.75_rk], [2,1]) - call nurbs%insert_knots(2, [0.25_rk, 0.75_rk], [2,1]) + call nurbs%set([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [1,1], [-1, -1], [-1, -1], Xc, Wc) + call bsp%set([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [1,1], [-1, -1], [-1, -1], Xc) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 13") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 14") + + call nurbs%set([2,2], Xc, Wc) + call bsp%set([2,2], Xc) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 15") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 16") + + + call nurbs%create(Xt1 = nurbs%get_Xt(1), Xt2 = nurbs%get_Xt(2)) + call bsp%create(Xt1 = bsp%get_Xt(1), Xt2 = nurbs%get_Xt(2)) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 17") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 18") + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 19") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 20") + + call ut%check(res=nurbs%get_Xc(1), expected=Xc(1,:), tol=1e-5_rk, msg="test_nurbs_surface: 21") + call ut%check(res=bsp%get_Xc(1), expected=Xc(1,:), tol=1e-5_rk, msg="test_nurbs_surface: 22") + + call ut%check(res=nurbs%get_Xc(1,1), expected=Xc(1,1), tol=1e-5_rk, msg="test_nurbs_surface: 23") + call ut%check(res=bsp%get_Xc(1,1), expected=Xc(1,1), tol=1e-5_rk, msg="test_nurbs_surface: 24") + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 25") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 26") + + call ut%check(res=nurbs%get_Xg(1), expected=Xg(1,:), tol=1e-5_rk, msg="test_nurbs_surface: 27") + call ut%check(res=bsp%get_Xg(1), expected=Xgb(1,:), tol=1e-5_rk, msg="test_nurbs_surface: 28") + + call ut%check(res=nurbs%get_Xg(1,1), expected=Xg(1,1), tol=1e-5_rk, msg="test_nurbs_surface: 29") + call ut%check(res=bsp%get_Xg(1,1), expected=Xgb(1,1), tol=1e-5_rk, msg="test_nurbs_surface: 30") + + call ut%check(res=nurbs%get_Wc(), expected=Wc, tol=1e-5_rk, msg="test_nurbs_surface: 31") + + call ut%check(res=nurbs%get_Wc(1), expected=Wc(1), tol=1e-5_rk, msg="test_nurbs_surface: 32") + + call ut%check(res=nurbs%get_knot(1), expected=knot1, tol=1e-5_rk, msg="test_nurbs_surface: 33") + call ut%check(res=bsp%get_knot(1), expected=knot1, tol=1e-5_rk, msg="test_nurbs_surface: 34") + + call ut%check(res=nurbs%get_knot(1,1), expected=knot1(1), tol=1e-5_rk, msg="test_nurbs_surface: 35") + call ut%check(res=bsp%get_knot(1,1), expected=knot1(1), tol=1e-5_rk, msg="test_nurbs_surface: 36") + + ! call ut%check(res=nurbs%get_ng(), expected=size(Xg,1), msg="test_nurbs_surface: 37") + ! call ut%check(res=bsp%get_ng(), expected=size(Xgb,1), msg="test_nurbs_surface: 38") + + call ut%check(res=nurbs%get_degree(1), expected=1, msg="test_nurbs_surface: 39") + call ut%check(res=bsp%get_degree(1), expected=1, msg="test_nurbs_surface: 40") + + call ut%check(res=nurbs%get_multiplicity(1), expected=[2,2], msg="test_nurbs_surface: 41") + call ut%check(res=bsp%get_multiplicity(1), expected=[2,2], msg="test_nurbs_surface: 42") + + call ut%check(res=nurbs%get_continuity(1), expected=[-1,-1], msg="test_nurbs_surface: 43") + call ut%check(res=bsp%get_continuity(1), expected=[-1,-1], msg="test_nurbs_surface: 44") + + ! call ut%check(res=nurbs%get_nc(), expected=size(Xc,1), msg="test_nurbs_surface: 45") + ! call ut%check(res=bsp%get_nc(), expected=size(Xc,1), msg="test_nurbs_surface: 46") + + call nurbs%cmp_nc() + call bsp%cmp_nc() + + ! call ut%check(res=nurbs%get_nc(), expected=size(Xc,1), msg="test_nurbs_surface: 47") + ! call ut%check(res=bsp%get_nc(), expected=size(Xc,1), msg="test_nurbs_surface: 48") + + elemConn = nurbs%cmp_elem_Xc_vis([1,1]) + call nurbs%set_elem_Xc_vis(elemConn) + call ut%check(res=nurbs%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_surface: 49") + deallocate(elemConn) + elemConn = nurbs%cmp_elem_Xc_vis() + call nurbs%set_elem_Xc_vis(elemConn) + call ut%check(res=nurbs%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_surface: 50") + deallocate(elemConn) + elemConn = bsp%cmp_elem_Xc_vis([1,1]) + call bsp%set_elem_Xc_vis(elemConn) + call ut%check(res=bsp%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_surface: 51") + deallocate(elemConn) + elemConn = bsp%cmp_elem_Xc_vis() + call bsp%set_elem_Xc_vis(elemConn) + call ut%check(res=bsp%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_surface: 52") + deallocate(elemConn) + + elemConn = nurbs%cmp_elem_Xg_vis([1,1]) + call nurbs%set_elem_Xg_vis(elemConn) + call ut%check(res=nurbs%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_surface: 53") + deallocate(elemConn) + elemConn = nurbs%cmp_elem_Xg_vis() + call nurbs%set_elem_Xg_vis(elemConn) + call ut%check(res=nurbs%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_surface: 54") + deallocate(elemConn) + elemConn = bsp%cmp_elem_Xg_vis([1,1]) + call bsp%set_elem_Xg_vis(elemConn) + call ut%check(res=bsp%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_surface: 55") + deallocate(elemConn) + elemConn = bsp%cmp_elem_Xg_vis() + call bsp%set_elem_Xg_vis(elemConn) + call ut%check(res=bsp%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_surface: 56") + deallocate(elemConn) + + elemConn = nurbs%cmp_elem() + call nurbs%set_elem(elemConn) + call ut%check(res=nurbs%get_elem(), expected=elemConn, msg="test_nurbs_surface: 57") + deallocate(elemConn) + elemConn = bsp%cmp_elem() + call bsp%set_elem(elemConn) + call ut%check(res=bsp%get_elem(), expected=elemConn, msg="test_nurbs_surface: 58") + deallocate(elemConn) + + call nurbs%modify_Xc(Xc(1,1), 1,1) + call bsp%modify_Xc(Xc(1,1), 1,1) + + call nurbs%modify_Wc(Wc(1),1) + + call nurbs%create() + call bsp%create() + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 59") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 60") + + call nurbs%basis(res1=30, res2=30, Tgc=Tgc) + call bsp%basis(res1=30, res2=30, Tgc=Tgc) - call bsp%insert_knots(1, [0.25_rk, 0.75_rk], [2,1]) - call bsp%insert_knots(2, [0.25_rk, 0.75_rk], [2,1]) + call nurbs%basis(Xt=[0.0_rk, 0.0_rk], Tgc=Tgc1) + call bsp%basis(Xt=[0.0_rk, 0.0_rk], Tgc=Tgc1b) + + call nurbs%basis(Xt1=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], Xt2=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], Tgc=Tgc) + call bsp%basis(Xt1=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], Xt2=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], Tgc=Tgc) + + call nurbs%derivative(res1=30, res2=30, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative(res1=30, res2=30, dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative(& + Xt1=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], Xt2=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], dTgc=dTgc, Tgc=Tgc) + call bsp%derivative(& + Xt1=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], Xt2=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative(Xt=[0.0_rk,0.0_rk], dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative(Xt=[0.0_rk,0.0_rk], dTgc=dTgc1b, Tgc=Tgc1b) + + call nurbs%derivative(Xt=[0.0_rk,0.0_rk], dTgc=dTgc1b, Tgc=Tgc1b, elem=[1,2,3]) + call bsp%derivative(Xt=[0.0_rk,0.0_rk], dTgc=dTgc1b, Tgc=Tgc1b, elem=[1,2,3]) + + call nurbs%derivative2(res1=30, res2=30, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative2(res1=30, res2=30, d2Tgc=d2Tgcb, dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative2(& + Xt1=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], Xt2=[(real(i-1, rk) / real(30-1, rk), i=1, 30)],& + d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative2(& + Xt1=[(real(i-1, rk) / real(30-1, rk), i=1, 30)], Xt2=[(real(i-1, rk) / real(30-1, rk), i=1, 30)],& + d2Tgc=d2Tgcb, dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative2(Xt=[0.0_rk,0.0_rk], d2Tgc=d2Tgc1, dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative2(Xt=[0.0_rk,0.0_rk], d2Tgc=d2Tgc1b, dTgc=dTgc1b, Tgc=Tgc1b) + + call nurbs%derivative2(Xt=[0.0_rk,0.0_rk], d2Tgc=d2Tgc1, dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative2(Xt=[0.0_rk,0.0_rk], d2Tgc=d2Tgc1b, dTgc=dTgc1b, Tgc=Tgc1b) + + call nurbs%rotate_Xc(45.0_rk, 0.0_rk, 0.0_rk) + call nurbs%rotate_Xc(-45.0_rk, 0.0_rk, 0.0_rk) + + call bsp%rotate_Xc(45.0_rk, 0.0_rk, 0.0_rk) + call bsp%rotate_Xc(-45.0_rk, 0.0_rk, 0.0_rk) + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 61") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 62") + + call nurbs%rotate_Xc(0.0_rk, 45.0_rk, 0.0_rk) + call nurbs%rotate_Xc(0.0_rk, -45.0_rk, 0.0_rk) + + call bsp%rotate_Xc(0.0_rk, 45.0_rk, 0.0_rk) + call bsp%rotate_Xc(0.0_rk, -45.0_rk, 0.0_rk) + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 63") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 64") + + call nurbs%rotate_Xc(0.0_rk, 0.0_rk, 45.0_rk) + call nurbs%rotate_Xc(0.0_rk, 0.0_rk, -45.0_rk) + + call bsp%rotate_Xc(0.0_rk, 0.0_rk, 45.0_rk) + call bsp%rotate_Xc(0.0_rk, 0.0_rk, -45.0_rk) + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 65") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 66") + + call nurbs%rotate_Xg(45.0_rk, 0.0_rk, 0.0_rk) + call nurbs%rotate_Xg(-45.0_rk, 0.0_rk, 0.0_rk) + + call bsp%rotate_Xg(45.0_rk, 0.0_rk, 0.0_rk) + call bsp%rotate_Xg(-45.0_rk, 0.0_rk, 0.0_rk) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 67") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 68") + + call nurbs%rotate_Xg(0.0_rk, 45.0_rk, 0.0_rk) + call nurbs%rotate_Xg(0.0_rk, -45.0_rk, 0.0_rk) + + call bsp%rotate_Xg(0.0_rk, 45.0_rk, 0.0_rk) + call bsp%rotate_Xg(0.0_rk, -45.0_rk, 0.0_rk) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 69") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 70") + + call nurbs%rotate_Xg(0.0_rk, 0.0_rk, 45.0_rk) + call nurbs%rotate_Xg(0.0_rk, 0.0_rk, -45.0_rk) + + call bsp%rotate_Xg(0.0_rk, 0.0_rk, 45.0_rk) + call bsp%rotate_Xg(0.0_rk, 0.0_rk, -45.0_rk) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 71") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 72") + + call nurbs%translate_Xc([5.0_rk, 5.0_rk, 5.0_rk]) + call nurbs%translate_Xc([-5.0_rk, -5.0_rk, -5.0_rk]) + + call bsp%translate_Xc([5.0_rk, 5.0_rk, 5.0_rk]) + call bsp%translate_Xc([-5.0_rk, -5.0_rk, -5.0_rk]) + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 73") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_surface: 74") + + call nurbs%translate_Xg([5.0_rk, 5.0_rk, 5.0_rk]) + call nurbs%translate_Xg([-5.0_rk, -5.0_rk, -5.0_rk]) + + call bsp%translate_Xg([5.0_rk, 5.0_rk, 5.0_rk]) + call bsp%translate_Xg([-5.0_rk, -5.0_rk, -5.0_rk]) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 75") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 76") + + call nurbs%export_Xc("vtk/test_nurbs_surface_Xc.vtk") + call nurbs%export_Xg("vtk/test_nurbs_surface_Xg.vtk") + + call bsp%export_Xc("vtk/test_nurbs_surface_Xc.vtk") + call bsp%export_Xg("vtk/test_nurbs_surface_Xg.vtk") + + call nurbs%insert_knots(1, [0.25_rk, 0.75_rk], [1,1]) + call nurbs%insert_knots(2, [0.25_rk, 0.75_rk], [1,1]) + + call bsp%insert_knots(1, [0.25_rk, 0.75_rk], [1,1]) + call bsp%insert_knots(2, [0.25_rk, 0.75_rk], [1,1]) call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 01") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 02") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 77") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 78") call nurbs%elevate_degree(1, 2) call nurbs%elevate_degree(2, 2) @@ -49,8 +324,8 @@ program test_nurbs_surface call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 03") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 04") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 79") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 80") call nurbs%remove_knots(1, [0.25_rk, 0.75_rk], [2,1]) call nurbs%remove_knots(2, [0.25_rk, 0.75_rk], [2,1]) @@ -61,37 +336,18 @@ program test_nurbs_surface call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 05") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 06") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_surface: 81") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_surface: 82") + + + call nurbs%set_tetragon([2.0_rk, 2.0_rk], [2,2]) + call bsp%set_tetragon([2.0_rk, 2.0_rk], [2,2], [1.0_rk,1.0_rk,0.9_rk,0.9_rk]) + call nurbs%set_ring([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk, 2.0_rk) + call nurbs%set_half_ring([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk, 2.0_rk) + call nurbs%set_C([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk, 2.0_rk) call nurbs%finalize() call bsp%finalize() deallocate(Xc, Wc, Xg, Xgb) -contains - - !----------------------------------------------------------------------------- - function generate_Xc(num_rows, num_cols, peak_height) result(control_points) - integer, intent(in) :: num_rows, num_cols - real(rk), intent(in) :: peak_height - real(rk), allocatable :: control_points(:,:) - integer :: i, j - real(rk) :: x_spacing, y_spacing, x_offset, y_offset - x_spacing = 1.0_rk / real(num_cols - 1) - y_spacing = 1.0_rk / real(num_rows - 1) - x_offset = -0.5_rk - y_offset = -0.5_rk - allocate(control_points(num_rows * num_cols, 3)) - do i = 1, num_rows - do j = 1, num_cols - control_points((i - 1) * num_cols + j, 1) = x_offset + real(j - 1) * x_spacing - control_points((i - 1) * num_cols + j, 2) = y_offset + real(i - 1) * y_spacing - control_points((i - 1) * num_cols + j, 3) = & - peak_height * exp(-((control_points((i - 1) * num_cols + j, 1) ** 2) & - + (control_points((i - 1) * num_cols + j, 2) ** 2))) + 0.5_rk * peak_height * 0.2_rk - end do - end do - end function - !----------------------------------------------------------------------------- - end program diff --git a/test/test_nurbs_volume.f90 b/test/test_nurbs_volume.f90 index eeefc18e4..a43e1822a 100644 --- a/test/test_nurbs_volume.f90 +++ b/test/test_nurbs_volume.f90 @@ -7,13 +7,27 @@ program test_nurbs_volume type(nurbs_volume) :: nurbs, bsp real(rk), allocatable :: Xc(:,:), Wc(:) real(rk), allocatable :: Xg(:,:), Xgb(:,:) - real(rk) :: knot1(4), knot2(4), knot3(4) + integer, allocatable :: elemConn(:,:) + real(rk), allocatable :: Tgc(:,:), dTgc(:,:,:), Tgcb(:,:), dTgcb(:,:,:), d2Tgc(:,:,:), d2Tgcb(:,:,:) + real(rk), allocatable :: Tgc1(:), dTgc1(:,:), Tgc1b(:), dTgc1b(:,:), d2Tgc1(:,:), d2Tgc1b(:,:) + real(rk) :: knot1(4), knot2(4), knot3(4), volume, volumeb + integer :: i, id + real(rk), allocatable :: nearest_Xg(:), nearest_Xt(:) type(unit_test) :: ut - Xc = generate_Xc(5.0_rk) + allocate(Xc(8, 3)) + Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk] + Xc(2,:) = [5.0_rk, 0.0_rk, 0.0_rk] + Xc(3,:) = [0.0_rk, 5.0_rk, 0.0_rk] + Xc(4,:) = [5.0_rk, 5.0_rk, 0.0_rk] + Xc(5,:) = [0.0_rk, 0.0_rk, 5.0_rk] + Xc(6,:) = [5.0_rk, 0.0_rk, 5.0_rk] + Xc(7,:) = [0.0_rk, 5.0_rk, 5.0_rk] + Xc(8,:) = [5.0_rk, 5.0_rk, 5.0_rk] - allocate(Wc(size(Xc,1)), source=1.0_rk) - Wc(2) = 5.0_rk + allocate(Wc(size(Xc, 1))) + Wc = 1.0_rk + Wc(2) = 0.9_rk knot1 = [0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk] knot2 = [0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk] @@ -22,12 +36,292 @@ program test_nurbs_volume call nurbs%set(knot1, knot2, knot3, Xc, Wc) call bsp%set(knot1, knot2, knot3, Xc) + call nurbs%set(degree=nurbs%get_degree(), nc=nurbs%get_nc(), Xc=nurbs%get_Xc(), Wc=nurbs%get_Wc()) + call bsp%set(degree=nurbs%get_degree(), nc=nurbs%get_nc(), Xc=nurbs%get_Xc()) + call nurbs%create(20, 20, 20) call bsp%create(20, 20, 20) + call nurbs%cmp_volume(volume) + call bsp%cmp_volume(volumeb) + + call ut%check(res=volume, expected=125.0_rk, tol=1e-5_rk, msg="test_nurbs_volume: 01") + call ut%check(res=volumeb, expected=125.0_rk, tol=1e-5_rk, msg="test_nurbs_volume: 02") + + call nurbs%nearest_point([0.0_rk, 0.0_rk, -0.5_rk], nearest_Xg, nearest_Xt, id) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_volume: 03") + call ut%check(res=nearest_Xt, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_volume: 04") + call ut%check(res=id, expected=1, msg="test_nurbs_volume: 05") + + call bsp%nearest_point([0.0_rk, 0.0_rk, -0.5_rk], nearest_Xg, nearest_Xt, id) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_volume: 06") + call ut%check(res=nearest_Xt, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_volume: 07") + call ut%check(res=id, expected=1, msg="test_nurbs_volume: 08") + + call nurbs%nearest_point2([0.0_rk, 0.0_rk, -0.5_rk], 1e-10_rk, 10, nearest_Xt, nearest_Xg) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_volume: 09") + call ut%check(res=nearest_Xt, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_volume: 10") + + call bsp%nearest_point2([0.0_rk, 0.0_rk, -0.5_rk], 1e-10_rk, 10, nearest_Xt, nearest_Xg) + call ut%check(res=nearest_Xg, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_volume: 11") + call ut%check(res=nearest_Xt, expected=[0.0_rk, 0.0_rk, 0.0_rk], tol=1e-5_rk, msg="test_nurbs_volume: 12") + Xg = nurbs%get_Xg() Xgb = bsp%get_Xg() + call nurbs%set([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [1,1,1], [-1, -1], [-1, -1], [-1, -1], Xc, Wc) + call bsp%set([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [1,1,1], [-1, -1], [-1, -1], [-1, -1], Xc) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 13") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 14") + + call nurbs%set([2,2,2], Xc, Wc) + call bsp%set([2,2,2], Xc) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 15") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 16") + + + call nurbs%create(Xt1 = nurbs%get_Xt(1), Xt2 = nurbs%get_Xt(2)) + call bsp%create(Xt1 = bsp%get_Xt(1), Xt2 = nurbs%get_Xt(2)) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 17") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 18") + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 19") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 20") + + call ut%check(res=nurbs%get_Xc(1), expected=Xc(1,:), tol=1e-5_rk, msg="test_nurbs_volume: 21") + call ut%check(res=bsp%get_Xc(1), expected=Xc(1,:), tol=1e-5_rk, msg="test_nurbs_volume: 22") + + call ut%check(res=nurbs%get_Xc(1,1), expected=Xc(1,1), tol=1e-5_rk, msg="test_nurbs_volume: 23") + call ut%check(res=bsp%get_Xc(1,1), expected=Xc(1,1), tol=1e-5_rk, msg="test_nurbs_volume: 24") + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 25") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 26") + + call ut%check(res=nurbs%get_Xg(1), expected=Xg(1,:), tol=1e-5_rk, msg="test_nurbs_volume: 27") + call ut%check(res=bsp%get_Xg(1), expected=Xgb(1,:), tol=1e-5_rk, msg="test_nurbs_volume: 28") + + call ut%check(res=nurbs%get_Xg(1,1), expected=Xg(1,1), tol=1e-5_rk, msg="test_nurbs_volume: 29") + call ut%check(res=bsp%get_Xg(1,1), expected=Xgb(1,1), tol=1e-5_rk, msg="test_nurbs_volume: 30") + + call ut%check(res=nurbs%get_Wc(), expected=Wc, tol=1e-5_rk, msg="test_nurbs_volume: 31") + + call ut%check(res=nurbs%get_Wc(1), expected=Wc(1), tol=1e-5_rk, msg="test_nurbs_volume: 32") + + call ut%check(res=nurbs%get_knot(1), expected=knot1, tol=1e-5_rk, msg="test_nurbs_volume: 33") + call ut%check(res=bsp%get_knot(1), expected=knot1, tol=1e-5_rk, msg="test_nurbs_volume: 34") + + call ut%check(res=nurbs%get_knot(1,1), expected=knot1(1), tol=1e-5_rk, msg="test_nurbs_volume: 35") + call ut%check(res=bsp%get_knot(1,1), expected=knot1(1), tol=1e-5_rk, msg="test_nurbs_volume: 36") + + ! call ut%check(res=nurbs%get_ng(), expected=size(Xg,1), msg="test_nurbs_volume: 37") + ! call ut%check(res=bsp%get_ng(), expected=size(Xgb,1), msg="test_nurbs_volume: 38") + + call ut%check(res=nurbs%get_degree(1), expected=1, msg="test_nurbs_volume: 39") + call ut%check(res=bsp%get_degree(1), expected=1, msg="test_nurbs_volume: 40") + + call ut%check(res=nurbs%get_multiplicity(1), expected=[2,2], msg="test_nurbs_volume: 41") + call ut%check(res=bsp%get_multiplicity(1), expected=[2,2], msg="test_nurbs_volume: 42") + + call ut%check(res=nurbs%get_continuity(1), expected=[-1,-1], msg="test_nurbs_volume: 43") + call ut%check(res=bsp%get_continuity(1), expected=[-1,-1], msg="test_nurbs_volume: 44") + + ! call ut%check(res=nurbs%get_nc(), expected=size(Xc,1), msg="test_nurbs_volume: 45") + ! call ut%check(res=bsp%get_nc(), expected=size(Xc,1), msg="test_nurbs_volume: 46") + + call nurbs%cmp_nc() + call bsp%cmp_nc() + + ! call ut%check(res=nurbs%get_nc(), expected=size(Xc,1), msg="test_nurbs_volume: 47") + ! call ut%check(res=bsp%get_nc(), expected=size(Xc,1), msg="test_nurbs_volume: 48") + + elemConn = nurbs%cmp_elem_Xc_vis([1,1,1]) + call nurbs%set_elem_Xc_vis(elemConn) + call ut%check(res=nurbs%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_volume: 49") + deallocate(elemConn) + elemConn = nurbs%cmp_elem_Xc_vis() + call nurbs%set_elem_Xc_vis(elemConn) + call ut%check(res=nurbs%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_volume: 50") + deallocate(elemConn) + elemConn = bsp%cmp_elem_Xc_vis([1,1,1]) + call bsp%set_elem_Xc_vis(elemConn) + call ut%check(res=bsp%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_volume: 51") + deallocate(elemConn) + elemConn = bsp%cmp_elem_Xc_vis() + call bsp%set_elem_Xc_vis(elemConn) + call ut%check(res=bsp%get_elem_Xc_vis(), expected=elemConn, msg="test_nurbs_volume: 52") + deallocate(elemConn) + + elemConn = nurbs%cmp_elem_Xg_vis([1,1,1]) + call nurbs%set_elem_Xg_vis(elemConn) + call ut%check(res=nurbs%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_volume: 53") + deallocate(elemConn) + elemConn = nurbs%cmp_elem_Xg_vis() + call nurbs%set_elem_Xg_vis(elemConn) + call ut%check(res=nurbs%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_volume: 54") + deallocate(elemConn) + elemConn = bsp%cmp_elem_Xg_vis([1,1,1]) + call bsp%set_elem_Xg_vis(elemConn) + call ut%check(res=bsp%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_volume: 55") + deallocate(elemConn) + elemConn = bsp%cmp_elem_Xg_vis() + call bsp%set_elem_Xg_vis(elemConn) + call ut%check(res=bsp%get_elem_Xg_vis(), expected=elemConn, msg="test_nurbs_volume: 56") + deallocate(elemConn) + + elemConn = nurbs%cmp_elem() + call nurbs%set_elem(elemConn) + call ut%check(res=nurbs%get_elem(), expected=elemConn, msg="test_nurbs_volume: 57") + deallocate(elemConn) + elemConn = bsp%cmp_elem() + call bsp%set_elem(elemConn) + call ut%check(res=bsp%get_elem(), expected=elemConn, msg="test_nurbs_volume: 58") + deallocate(elemConn) + + call nurbs%modify_Xc(Xc(1,1), 1,1) + call bsp%modify_Xc(Xc(1,1), 1,1) + + call nurbs%modify_Wc(Wc(1),1) + + call nurbs%create() + call bsp%create() + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 59") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 60") + + call nurbs%basis(res1=20, res2=20, res3=20, Tgc=Tgc) + call bsp%basis(res1=20, res2=20, res3=20, Tgc=Tgc) + + call nurbs%basis(Xt=[0.0_rk, 0.0_rk, 0.0_rk], Tgc=Tgc1) + call bsp%basis(Xt=[0.0_rk, 0.0_rk, 0.0_rk], Tgc=Tgc1b) + + call nurbs%basis(& + Xt1=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt2=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt3=[(real(i-1, rk) / real(20-1, rk), i=1, 20)], Tgc=Tgc) + call bsp%basis(& + Xt1=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt2=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt3=[(real(i-1, rk) / real(20-1, rk), i=1, 20)], Tgc=Tgc) + + call nurbs%derivative(res1=20, res2=20, res3=20, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative(res1=20, res2=20, res3=20, dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative(& + Xt1=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt2=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt3=[(real(i-1, rk) / real(20-1, rk), i=1, 20)], dTgc=dTgc, Tgc=Tgc) + call bsp%derivative(& + Xt1=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt2=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt3=[(real(i-1, rk) / real(20-1, rk), i=1, 20)], dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative(Xt=[0.0_rk,0.0_rk,0.0_rk], dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative(Xt=[0.0_rk,0.0_rk,0.0_rk], dTgc=dTgc1b, Tgc=Tgc1b) + + call nurbs%derivative(Xt=[0.0_rk,0.0_rk,0.0_rk], dTgc=dTgc1b, Tgc=Tgc1b, elem=[1,2,3]) + call bsp%derivative(Xt=[0.0_rk,0.0_rk,0.0_rk], dTgc=dTgc1b, Tgc=Tgc1b, elem=[1,2,3]) + + call nurbs%derivative2(res1=20, res2=20, res3=20, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative2(res1=20, res2=20, res3=20, d2Tgc=d2Tgcb, dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative2(& + Xt1=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt2=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt3=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc) + call bsp%derivative2(& + Xt1=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt2=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + Xt3=[(real(i-1, rk) / real(20-1, rk), i=1, 20)],& + d2Tgc=d2Tgcb, dTgc=dTgcb, Tgc=Tgcb) + + call nurbs%derivative2(Xt=[0.0_rk,0.0_rk,0.0_rk], d2Tgc=d2Tgc1, dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative2(Xt=[0.0_rk,0.0_rk,0.0_rk], d2Tgc=d2Tgc1b, dTgc=dTgc1b, Tgc=Tgc1b) + + call nurbs%derivative2(Xt=[0.0_rk,0.0_rk,0.0_rk], d2Tgc=d2Tgc1, dTgc=dTgc1b, Tgc=Tgc1b) + call bsp%derivative2(Xt=[0.0_rk,0.0_rk,0.0_rk], d2Tgc=d2Tgc1b, dTgc=dTgc1b, Tgc=Tgc1b) + + call nurbs%rotate_Xc(45.0_rk, 0.0_rk, 0.0_rk) + call nurbs%rotate_Xc(-45.0_rk, 0.0_rk, 0.0_rk) + + call bsp%rotate_Xc(45.0_rk, 0.0_rk, 0.0_rk) + call bsp%rotate_Xc(-45.0_rk, 0.0_rk, 0.0_rk) + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 61") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 62") + + call nurbs%rotate_Xc(0.0_rk, 45.0_rk, 0.0_rk) + call nurbs%rotate_Xc(0.0_rk, -45.0_rk, 0.0_rk) + + call bsp%rotate_Xc(0.0_rk, 45.0_rk, 0.0_rk) + call bsp%rotate_Xc(0.0_rk, -45.0_rk, 0.0_rk) + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 63") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 64") + + call nurbs%rotate_Xc(0.0_rk, 0.0_rk, 45.0_rk) + call nurbs%rotate_Xc(0.0_rk, 0.0_rk, -45.0_rk) + + call bsp%rotate_Xc(0.0_rk, 0.0_rk, 45.0_rk) + call bsp%rotate_Xc(0.0_rk, 0.0_rk, -45.0_rk) + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 65") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 66") + + call nurbs%rotate_Xg(45.0_rk, 0.0_rk, 0.0_rk) + call nurbs%rotate_Xg(-45.0_rk, 0.0_rk, 0.0_rk) + + call bsp%rotate_Xg(45.0_rk, 0.0_rk, 0.0_rk) + call bsp%rotate_Xg(-45.0_rk, 0.0_rk, 0.0_rk) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 67") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 68") + + call nurbs%rotate_Xg(0.0_rk, 45.0_rk, 0.0_rk) + call nurbs%rotate_Xg(0.0_rk, -45.0_rk, 0.0_rk) + + call bsp%rotate_Xg(0.0_rk, 45.0_rk, 0.0_rk) + call bsp%rotate_Xg(0.0_rk, -45.0_rk, 0.0_rk) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 69") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 70") + + call nurbs%rotate_Xg(0.0_rk, 0.0_rk, 45.0_rk) + call nurbs%rotate_Xg(0.0_rk, 0.0_rk, -45.0_rk) + + call bsp%rotate_Xg(0.0_rk, 0.0_rk, 45.0_rk) + call bsp%rotate_Xg(0.0_rk, 0.0_rk, -45.0_rk) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 71") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 72") + + call nurbs%translate_Xc([5.0_rk, 5.0_rk, 5.0_rk]) + call nurbs%translate_Xc([-5.0_rk, -5.0_rk, -5.0_rk]) + + call bsp%translate_Xc([5.0_rk, 5.0_rk, 5.0_rk]) + call bsp%translate_Xc([-5.0_rk, -5.0_rk, -5.0_rk]) + + call ut%check(res=nurbs%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 73") + call ut%check(res=bsp%get_Xc(), expected=Xc, tol=1e-5_rk, msg="test_nurbs_volume: 74") + + call nurbs%translate_Xg([5.0_rk, 5.0_rk, 5.0_rk]) + call nurbs%translate_Xg([-5.0_rk, -5.0_rk, -5.0_rk]) + + call bsp%translate_Xg([5.0_rk, 5.0_rk, 5.0_rk]) + call bsp%translate_Xg([-5.0_rk, -5.0_rk, -5.0_rk]) + + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 75") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 76") + + call nurbs%export_Xc("vtk/test_nurbs_volume_Xc.vtk") + call nurbs%export_Xg("vtk/test_nurbs_volume_Xg.vtk") + + call bsp%export_Xc("vtk/test_nurbs_volume_Xc.vtk") + call bsp%export_Xg("vtk/test_nurbs_volume_Xg.vtk") + call nurbs%insert_knots(1, [0.25_rk, 0.75_rk], [1,1]) call nurbs%insert_knots(2, [0.25_rk, 0.75_rk], [1,1]) call nurbs%insert_knots(3, [0.25_rk, 0.75_rk], [1,1]) @@ -39,8 +333,8 @@ program test_nurbs_volume call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 01") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 02") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 77") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 78") call nurbs%elevate_degree(1, 2) call nurbs%elevate_degree(2, 2) @@ -53,46 +347,32 @@ program test_nurbs_volume call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 03") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 04") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 79") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 80") - call nurbs%remove_knots(1, [0.25_rk, 0.75_rk], [1,1]) - call nurbs%remove_knots(2, [0.25_rk, 0.75_rk], [1,1]) - call nurbs%remove_knots(3, [0.25_rk, 0.75_rk], [1,1]) + call nurbs%remove_knots(1, [0.25_rk, 0.75_rk], [2,1]) + call nurbs%remove_knots(2, [0.25_rk, 0.75_rk], [2,1]) + call nurbs%remove_knots(3, [0.25_rk, 0.75_rk], [2,1]) - call bsp%remove_knots(1, [0.25_rk, 0.75_rk], [1,1]) - call bsp%remove_knots(2, [0.25_rk, 0.75_rk], [1,1]) - call bsp%remove_knots(3, [0.25_rk, 0.75_rk], [1,1]) + ! call bsp%remove_knots(1, [0.25_rk, 0.75_rk], [2,1]) + ! call bsp%remove_knots(2, [0.25_rk, 0.75_rk], [2,1]) + ! call bsp%remove_knots(3, [0.25_rk, 0.75_rk], [2,1]) call nurbs%create() call bsp%create() - call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 05") - call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 06") + call ut%check(res=nurbs%get_Xg(), expected=Xg, tol=1e-5_rk, msg="test_nurbs_volume: 81") + call ut%check(res=bsp%get_Xg(), expected=Xgb, tol=1e-5_rk, msg="test_nurbs_volume: 82") + + + call nurbs%set_hexahedron([2.0_rk, 2.0_rk, 2.0_rk], [2,2,2]) + call bsp%set_hexahedron([2.0_rk, 2.0_rk, 2.0_rk], [2,2,2], [1.0_rk,1.0_rk,0.9_rk,0.9_rk,1.0_rk,1.0_rk,1.0_rk,0.9_rk]) + call nurbs%set_ring([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk, 2.0_rk, 2.0_rk) + call nurbs%set_half_ring([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk, 2.0_rk, 2.0_rk) + call nurbs%set_C([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk, 2.0_rk, 2.0_rk) call nurbs%finalize() call bsp%finalize() deallocate(Xc, Wc, Xg, Xgb) -contains - - !----------------------------------------------------------------------------- - function generate_Xc(L) result(control_points) - implicit none - real(rk), intent(in) :: L - real(rk), allocatable :: control_points(:,:) - real(rk) :: L2 - L2 = L / 2.0_rk - allocate(control_points(8, 3)) - control_points(1,:) = [ L2, -L2, L2] - control_points(2,:) = [ L2, -L2, -L2] - control_points(3,:) = [-L2, -L2, L2] - control_points(4,:) = [-L2, -L2, -L2] - control_points(5,:) = [ L2, L2, L2] - control_points(6,:) = [ L2, L2, -L2] - control_points(7,:) = [-L2, L2, L2] - control_points(8,:) = [-L2, L2, -L2] - end function - !----------------------------------------------------------------------------- - end program